Changeset ad8e1b in git for factory/algext.cc


Ignore:
Timestamp:
Aug 6, 2008, 3:53:15 PM (16 years ago)
Author:
Hans Schönemann <hannes@…>
Branches:
(u'spielwiese', 'fe61d9c35bf7c61f2b6cbf1b56e25e2f08d536cc')
Children:
e8a0aaf34722c3d93f882019f7622c3419ed764d
Parents:
69b66ef25238d635e8046b0a178b0552b57b23c3
Message:
*hannes: QGCD, fieldGCD, etc


git-svn-id: file:///usr/local/Singular/svn/trunk@10949 2c84dea3-7e68-4137-9b89-c4e89433aadc
File:
1 edited

Legend:

Unmodified
Added
Removed

  • factory/algext.cc

    r69b66ef rad8e1b  
    1313#include "cf_generator.h"
    1414
    15 void tryEuclid( const CanonicalForm & A, const CanonicalForm & B, const CanonicalForm M, CanonicalForm & result, bool & fail )
    16 {
    17   CanonicalForm P;
    18   if( degree(A) > degree(B) )
    19   {
    20     P = A; result = B;
    21   }
    22   else
    23   {
    24     P = B; result = A;
    25   }
    26   if( P.isZero() ) // then result is zero, too
    27     return;
    28   CanonicalForm inv;
    29   if( result.isZero() )
    30   {
    31     tryInvert( Lc(P), M, inv, fail );
    32     if(fail)
    33       return;
    34     result = inv*P; // monify result
    35     return;
    36   }
    37   CanonicalForm rem;
    38   // here: degree(P) >= degree(result)
    39   while(true)
    40   {
    41     tryInvert( Lc(result), M, inv, fail );
    42     if(fail)
    43       return;
    44     // here: Lc(result) is invertible modulo M
    45     rem = P - Lc(P)*inv*result * power( P.mvar(), degree(P)-degree(result) );
    46     if( rem.isZero() )
    47     {
    48       result *= inv; // monify result
    49       return;
    50     }
    51     P = result;
    52     result = rem;
    53   }
    54 }
    55 
     15CanonicalForm reduce(const CanonicalForm & f, const CanonicalForm & M)
     16{ // polynomials in M.mvar() are considered coefficients
     17  // M univariate monic polynomial
     18  // the coefficients of f are reduced modulo M
     19  if(f.inBaseDomain() || f.level() < M.level())
     20    return f;
     21  if(f.level() == M.level())
     22  {
     23    if(f.degree() < M.degree())
     24      return f;
     25    CanonicalForm tmp = f;
     26    do
     27      tmp -= lc(tmp)*M*power(M.mvar(),tmp.degree()-M.degree());
     28    while(tmp.degree() >= M.degree());
     29    return tmp;
     30  }
     31  // here: f.level() > M.level()
     32  CanonicalForm result = 0;
     33  for(CFIterator i=f; i.hasTerms(); i++)
     34    result += reduce(i.coeff(),M) * power(f.mvar(),i.exp());
     35  return result;
     36}
    5637
    5738void tryInvert( const CanonicalForm & F, const CanonicalForm & M, CanonicalForm & inv, bool & fail )
    58 {
    59   // F, M are required to be "univariate" polynomials in an algebraic variable
     39{ // F, M are required to be "univariate" polynomials in an algebraic variable
    6040  // we try to invert F modulo M
     41  if(F.inBaseDomain())
     42  {
     43    if(F.isZero())
     44    {
     45      fail = true;
     46      return;
     47    }
     48    inv = 1/F;
     49    return;
     50  }
    6151  CanonicalForm b;
    6252  Variable a = M.mvar();
    6353  Variable x = Variable(1);
    6454  if(!extgcd( replacevar( F, a, x ), replacevar( M, a, x ), inv, b ).isOne())
    65   {
    66     printf("failed ");
    6755    fail = true;
    68   }
    6956  else
    70     inv = replacevar( inv, a, x); // change back to alg var
    71 }
    72 
     57    inv = replacevar( inv, x, a ); // change back to alg var
     58}
     59
     60void tryEuclid( const CanonicalForm & A, const CanonicalForm & B, const CanonicalForm & M, CanonicalForm & result, bool & fail )
     61{
     62  CanonicalForm P;
     63  if(A.inCoeffDomain())
     64  {
     65    tryInvert( A, M, P, fail );
     66    if(fail)
     67      return;
     68    result = 1;
     69    return;
     70  }
     71  if(B.inCoeffDomain())
     72  {
     73    tryInvert( B, M, P, fail );
     74    if(fail)
     75      return;
     76    result = 1;
     77    return;
     78  }
     79  // here: both not inCoeffDomain
     80  if( A.degree() > B.degree() )
     81  {
     82    P = A; result = B;
     83  }
     84  else
     85  {
     86    P = B; result = A;
     87  }
     88  CanonicalForm inv;
     89  if( result.isZero() )
     90  {
     91    tryInvert( Lc(P), M, inv, fail );
     92    if(fail)
     93      return;
     94    result = inv*P; // monify result (not reduced, yet)
     95    return;
     96  }
     97  Variable x = P.mvar();
     98  CanonicalForm rem;
     99  // here: degree(P) >= degree(result)
     100  while(true)
     101  {
     102    tryInvert( Lc(result), M, inv, fail );
     103    if(fail)
     104      return;
     105    // here: Lc(result) is invertible modulo M
     106    rem = reduce( P - Lc(P)*inv*result*power( x, P.degree(x)-result.degree(x) ), M );
     107    if( rem.isZero() )
     108    {
     109      result *= inv;
     110      return;
     111    }
     112    if(result.degree(x) >= rem.degree(x))
     113    {
     114      P = result;
     115      result = rem;
     116    }
     117    else
     118      P = rem;
     119  }
     120}
    73121
    74122bool hasFirstAlgVar( const CanonicalForm & f, Variable & a )
    … …  
    76124  if( f.inBaseDomain() ) // f has NO alg. variable
    77125    return false;
    78 
    79126  if( f.level()<0 ) // f has only alg. vars, so take the first one
    80127  {
    … …  
    85132    if( hasFirstAlgVar( i.coeff(), a ))
    86133      return true; // 'a' is already set
    87 
    88134  return false;
    89135}
    … …  
    183229
    184230
    185 CanonicalForm QGCD( const CanonicalForm & F, const CanonicalForm & G )
    186 { // F,G in Q(a)[x1,...,xn]
     231
     232CanonicalForm QGCD( const CanonicalForm & F, const CanonicalForm & G );
     233int * leadDeg(const CanonicalForm & f, int *degs);
     234bool isLess(int *a, int *b, int lower, int upper);
     235bool isEqual(int *a, int *b, int lower, int upper);
     236CanonicalForm firstLC(const CanonicalForm & f);
     237void tryCRA( const CanonicalForm & x1, const CanonicalForm & q1, const CanonicalForm & x2, const CanonicalForm & q2, const CanonicalForm & M, CanonicalForm & xnew, CanonicalForm & qnew, bool & fail );
     238void tryExtgcd( const CanonicalForm & F, const CanonicalForm & G, const CanonicalForm & M, CanonicalForm & result, CanonicalForm & s, CanonicalForm & t, bool & fail );
     239static CanonicalForm trycontent ( const CanonicalForm & f, const Variable & x, const CanonicalForm & M, bool & fail );
     240static CanonicalForm tryvcontent ( const CanonicalForm & f, const Variable & x, const CanonicalForm & M, bool & fail );
     241static CanonicalForm trycf_content ( const CanonicalForm & f, const CanonicalForm & g, const CanonicalForm & M, bool & fail );
     242static void tryDivide( const CanonicalForm & f, const CanonicalForm & g, const CanonicalForm & M, CanonicalForm & result, bool & divides, bool & fail );
     243
     244
     245void tryBrownGCD( const CanonicalForm & F, const CanonicalForm & G, const CanonicalForm & M, CanonicalForm & result, bool & fail )
     246{ // assume F,G are multivariate polys over Z/p(a) for big prime p, M "univariate" polynomial in an algebraic variable
     247  // M is assumed to be monic
    187248  if(F.isZero())
    188249  {
    189250    if(G.isZero())
    190       return G; // G is zero
     251    {
     252      result = G; // G is zero
     253      return;
     254    }
    191255    if(G.inCoeffDomain())
    192       return CanonicalForm(1);
    193     return G/Lc(G); // return monic G
    194   }
    195   if(G.isZero()) // F is non-zero
    196   {
    197     if(F.inCoeffDomain())
    198       return CanonicalForm(1);
    199     return F/Lc(F); // return monic F
    200   }
    201   if(F.inCoeffDomain() || G.inCoeffDomain())
    202     return CanonicalForm(1);
    203 
    204   CanonicalForm D;
    205   if (getCharacteristic()==0)
    206   {
    207     CanonicalForm f,g,tmp, M, q, Dp, cl, newD, newq;
    208     int p, i;
    209     int *bound, *other; // degree vectors
    210     bool fail;
    211     On(SW_RATIONAL);
    212     f = F * bCommonDen(F);
    213     g = G * bCommonDen(G);
    214     Variable a,b;
    215     if( !hasFirstAlgVar( f, a ) && !hasFirstAlgVar( g, b ))
    216     {
    217       // F and G are in Q[x1,...,xn], call...
    218       return gcd_poly_0( f, g );
    219     }
    220     if( b.level() > a.level() )
    221       a = b;
    222     // here: a is the biggest alg. var in f and g
    223     tmp = getMipo(a);
    224     M = tmp * bCommonDen(tmp);
    225     Off( SW_RATIONAL );
    226     // here: f, g in Z[y][x1,...,xn], M in Z[y] not necessarily monic
    227     // calculate upper bound for degree of gcd
    228     int mv = f.level();
    229     if(g.level() > mv)
    230       mv = g.level();
    231     // here: mv is level of the largest variable in f, g
    232     bound = new int[mv+1];
    233     other = new int[mv+1];
    234     for(int i=1; i<=mv; i++)
    235       bound[i] = other[i] = 0;
    236     bound = leadDeg(f,bound);
    237     other = leadDeg(g,other);
    238     for(int i=1; i<=mv; i++) // entry at i=0 is not visited
    239       if(other[i] < bound[i])
    240         bound[i] = other[i];
    241     // now bound points on the smaller vector
    242     cl = lc(M) * lc(f) * lc(g);
    243     q = 1;
    244     D = 0;
    245     for(int i=cf_getNumBigPrimes()-1; i>-1; i--)
    246     {
    247       p = cf_getBigPrime(i);
    248       if( mod( cl, p ).isZero() ) // skip lc-bad primes
    249         continue;
    250 
    251       fail = false;
    252       setCharacteristic(p);
    253       tryBrownGCD( mapinto(f), mapinto(g), mapinto(M), Dp, fail );
    254       setCharacteristic(0);
    255       if( fail ) // M splits in char p
    256         continue;
    257 
    258       for(int i=1; i<=mv; i++)
    259         other[i] = 0; // reset (this is necessary, because some entries may not be updated by call to leadDeg)
    260       other = leadDeg(Dp,other);
    261 
    262       if( other==0 ) // early termination
    263       {
    264         // Dp is in coeff domain
    265         CanonicalForm inv;
    266         tryInvert(Dp,M,inv,fail); // check if zero-divisor
    267         if(fail)
    268           continue;
    269         return CanonicalForm(1);
    270       }
    271 
    272       if( isLess(bound, other, 1, mv) ) // current prime unlucky
    273         continue;
    274 
    275       if( isLess(other, bound, 1, mv) ) // all previous primes unlucky
    276       {
    277         q = p;
    278         D = mapinto(Dp); // shortcut CRA
    279         for(int i=1; i<=mv; i++) // tighten bound
    280           bound[i] = other[i];
    281         continue;
    282       }
    283       chineseRemainder( D, q, mapinto(Dp), p, newD, newq );
    284       // newD = Dp mod p
    285       // newD = D mod q
    286       // newq = p*q
    287       q = newq;
    288       if( D != newD )
    289       {
    290         D = newD;
    291         continue;
    292       }
    293       On( SW_RATIONAL );
    294       tmp = Farey( D, q );
    295       if( fdivides( tmp, f ) && fdivides( tmp, g )) // trial division
    296       {
    297         Off( SW_RATIONAL );
    298         return tmp;
    299       }
    300       Off( SW_RATIONAL );
    301     }
    302   }
    303   // hopefully, we never reach this point
    304   Off( SW_USE_QGCD );
    305   D = gcd( F, G );
    306   On( SW_USE_QGCD );
    307   return D;
    308 }
    309 
    310 void tryBrownGCD( const CanonicalForm & F, const CanonicalForm & G, const CanonicalForm & M, CanonicalForm & result, bool & fail )
    311 {// assume F,G are multivariate polys over Z/p(a) for big prime p, M "univariate" polynomial in an algebraic variable
    312   printf("Brown ");
    313   if(F.isZero())
    314   {
    315     if(G.isZero())
    316     {
    317       result = G; // G is zero
    318       return;
    319     }
    320     if(G.inCoeffDomain())
    321256    {
    322257      tryInvert(G,M,result,fail);
     258      if(fail)
     259        return;
     260      result = 1;
    323261      return;
    324262    }
    … …  
    336274    {
    337275      tryInvert(F,M,result,fail);
     276      if(fail)
     277        return;
     278      result = 1;
    338279      return;
    339280    }
    … …  
    346287    return;
    347288  }
     289  // here: F,G both nonzero
    348290  if(F.inCoeffDomain())
    349291  {
    350292    tryInvert(F,M,result,fail);
     293    if(fail)
     294      return;
     295    result = 1;
    351296    return;
    352297  }
    … …  
    354299  {
    355300    tryInvert(G,M,result,fail);
     301    if(fail)
     302      return;
     303    result = 1;
    356304    return;
    357305  }
    … …  
    363311  CanonicalForm f=MM(F);
    364312  CanonicalForm g=MM(G);
    365   // here: f, g are compressed
     313  // here: f,g are compressed
    366314  // compute largest variable in f or g (least one is Variable(1))
    367315  int mv = f.level();
    … …  
    378326  }
    379327  // here: mv > 1
    380   CanonicalForm cf = vcontent(f, Variable(2)); // cf is univariate poly in var(1)
    381   CanonicalForm cg = vcontent(g, Variable(2));
     328  CanonicalForm cf = tryvcontent(f, Variable(2), M, fail); // cf is univariate poly in var(1)
     329  if(fail)
     330    return;
     331  CanonicalForm cg = tryvcontent(g, Variable(2), M, fail);
     332  if(fail)
     333    return;
    382334  CanonicalForm c;
    383335  tryEuclid(cf,cg,M,c,fail);
    384336  if(fail)
    385337    return;
    386   f/=cf;
    387   g/=cg;
     338  bool divides = true;
     339  CanonicalForm tmp;
     340  // f /= cf
     341  tryDivide(f,cf,M,tmp,divides,fail); // 'divides' is not checked
     342  if(fail)
     343    return;
     344  f = tmp;
     345  // g /= cg
     346  tryDivide(g,cg,M,tmp,divides,fail); // 'divides' is not checked
     347  if(fail)
     348    return;
     349  g = tmp;
    388350  if(f.inCoeffDomain())
    389351  {
    … …  
    391353    if(fail)
    392354      return;
    393     result = NN(result);
     355    result = NN(c);
    394356    return;
    395357  }
    … …  
    399361    if(fail)
    400362      return;
    401     result = NN(result);
     363    result = NN(c);
    402364    return;
    403365  }
    … …  
    412374  if(fail)
    413375    return;
    414   for(int i=2; i<=mv; i++) // entry at i=1 is not visited
     376  for(int i=2; i<=mv; i++) // entries at i=0,1 not visited
    415377    if(N[i] < L[i])
    416378      L[i] = N[i];
    … …  
    426388  {
    427389    alpha = gen.item();
    428     gamma_image = gamma(alpha, Variable(1)); // plug in alpha for var(1)
     390    gamma_image = reduce(gamma(alpha, Variable(1)),M); // plug in alpha for var(1)
    429391    if(gamma_image.isZero()) // skip lc-bad points var(1)-alpha
    430392      continue;
    … …  
    432394    if(fail)
    433395      return;
     396    g_image = reduce(g_image, M);
    434397    if(g_image.inCoeffDomain()) // early termination
    435398    {
    … …  
    447410      g_image /= lc(g_image); // make g_image monic over Z/p
    448411      g_image *= gamma_image; // multiply by multiple of image lc(gcd)
    449       tryCRA( g_image, Variable(1)-alpha, gm, m, gnew, mnew, fail );
     412      tryCRA( g_image, Variable(1)-alpha, gm, m, M, gnew, mnew, fail );
    450413      // gnew = gm mod m
    451414      // gnew = g_image mod var(1)-alpha
    … …  
    456419      if(gnew == gm) // gnew did not change
    457420      {
    458         g_image = gm / vcontent(gm, Variable(2));
    459         if(fdivides(g_image,f) && fdivides(g_image,g)) // trial division
     421        cf = tryvcontent(gm, Variable(2), M, fail);
     422        if(fail)
     423          return;
     424        divides = true;
     425        tryDivide(gm,cf,M,g_image,divides,fail); // 'divides' is ignored
     426        if(fail)
     427          return;
     428        tryDivide(f,g_image,M,tmp,divides,fail); // trial division (f)
     429        if(fail)
     430          return;
     431        if(divides)
    460432        {
    461           result = NN(c*g_image);
    462           return;
     433          tryDivide(g,g_image,M,tmp,divides,fail); // trial division (g)
     434          if(fail)
     435            return;
     436          if(divides)
     437          {
     438            result = NN(c*g_image);
     439            return;
     440          }
    463441        }
    464442      }
    … …  
    470448      continue;
    471449
    472     if(isLess(dg_im, L, 2, mv)) // dg_im < L --> all previous points were unlucky
    473     {
    474       m = CanonicalForm(1); // reset
    475       gm = 0; // reset
    476       for(int i=2; i<=mv; i++) // tighten bound
    477         L[i] = dg_im[i];
    478     }
     450    // here: isLess(dg_im, L, 2, mv) --> all previous points were unlucky
     451    m = CanonicalForm(1); // reset
     452    gm = 0; // reset
     453    for(int i=2; i<=mv; i++) // tighten bound
     454      L[i] = dg_im[i];
    479455  }
    480456  // we are out of evaluation points
    481457  fail = true;
    482458}
     459
     460CanonicalForm QGCD( const CanonicalForm & F, const CanonicalForm & G )
     461{ // f,g in Q(a)[x1,...,xn]
     462  if(F.isZero())
     463  {
     464    if(G.isZero())
     465      return G; // G is zero
     466    if(G.inCoeffDomain())
     467      return CanonicalForm(1);
     468    return G/Lc(G); // return monic G
     469  }
     470  if(G.isZero()) // F is non-zero
     471  {
     472    if(F.inCoeffDomain())
     473      return CanonicalForm(1);
     474    return F/Lc(F); // return monic F
     475  }
     476  if(F.inCoeffDomain() || G.inCoeffDomain())
     477    return CanonicalForm(1);
     478  // here: both NOT inCoeffDomain
     479  CanonicalForm f, g, tmp, M, q, D, Dp, cl, newq, mipo;
     480  int p, i;
     481  int *bound, *other; // degree vectors
     482  bool fail;
     483  On( SW_RATIONAL ); // needed by bCommonDen
     484  f = F * bCommonDen(F);
     485  g = G * bCommonDen(G);
     486  Variable a, b;
     487  if(hasFirstAlgVar(f,a))
     488  {
     489    if(hasFirstAlgVar(g,b))
     490    {
     491      if(b.level() > a.level())
     492        a = b;
     493    }
     494  }
     495  else
     496  {
     497    if(!hasFirstAlgVar(g,a))// both not in extension
     498    {
     499      Off( SW_RATIONAL );
     500      Off( SW_USE_QGCD );
     501      tmp = gcd( F, G );
     502      On( SW_USE_QGCD );
     503      return tmp;
     504    }
     505  }
     506  // here: a is the biggest alg. var in f and g AND some of f,g is in extension
     507  // (in the sequel b is used to swap alg/poly vars)
     508  setReduce(a,false); // do not reduce expressions modulo mipo
     509  tmp = getMipo(a);
     510  M = tmp * bCommonDen(tmp);
     511  // here: f, g in Z[a][x1,...,xn], M in Z[a] not necessarily monic
     512  Off( SW_RATIONAL ); // needed by mod
     513  // calculate upper bound for degree vector of gcd
     514  int mv = f.level(); i = g.level();
     515  if(i > mv)
     516    mv = i;
     517  // here: mv is level of the largest variable in f, g
     518  b = Variable(mv+1);
     519  bound = new int[mv+1]; // 'bound' could be indexed from 0 to mv, but we will only use from 1 to mv
     520  other = new int[mv+1];
     521  for(int i=1; i<=mv; i++) // initialize 'bound', 'other' with zeros
     522    bound[i] = other[i] = 0;
     523  bound = leadDeg(f,bound); // 'bound' is set the leading degree vector of f
     524  other = leadDeg(g,other);
     525  for(int i=1; i<=mv; i++) // entry at i=0 not visited
     526    if(other[i] < bound[i])
     527      bound[i] = other[i];
     528  // now 'bound' is the smaller vector
     529  cl = lc(M) * lc(f) * lc(g);
     530  q = 1;
     531  D = 0;
     532  for( i=cf_getNumBigPrimes()-1; i>-1; i-- )
     533  {
     534    p = cf_getBigPrime(i);
     535    if( mod( cl, p ).isZero() ) // skip lc-bad primes
     536      continue;
     537    fail = false;
     538    setCharacteristic(p);
     539    mipo = mapinto(M);
     540    mipo /= mipo.lc();
     541    // here: mipo is monic
     542    tryBrownGCD( mapinto(f), mapinto(g), mipo, Dp, fail );
     543    Dp = reduce( Dp, mipo );
     544    setCharacteristic(0);
     545    if( fail ) // mipo splits in char p
     546      continue;
     547    if( Dp.inCoeffDomain() ) // early termination
     548    {
     549      tryInvert(Dp,mipo,tmp,fail); // check if zero divisor
     550      if(fail)
     551        continue;
     552      setReduce(a,true);
     553      return CanonicalForm(1);
     554    }
     555    // here: Dp NOT inCoeffDomain
     556    for(int i=1; i<=mv; i++)
     557      other[i] = 0; // reset (this is necessary, because some entries may not be updated by call to leadDeg)
     558    other = leadDeg(Dp,other);
     559   
     560    if(isEqual(bound, other, 1, mv)) // equal
     561    {
     562      chineseRemainder( D, q, replacevar( mapinto(Dp), a, b ), p, tmp, newq );
     563      // tmp = Dp mod p
     564      // tmp = D mod q
     565      // newq = p*q
     566      q = newq;
     567      if( D != tmp )
     568        D = tmp;
     569      On( SW_RATIONAL );
     570      tmp = replacevar( Farey( D, q ), b, a ); // Farey and switch back to alg var
     571      setReduce(a,true); // reduce expressions modulo mipo
     572      On( SW_RATIONAL ); // needed by fdivides
     573      if( fdivides( tmp, f ) && fdivides( tmp, g )) // trial division
     574      {
     575        Off( SW_RATIONAL );
     576        setReduce(a,true);
     577        return tmp;
     578      }
     579      Off( SW_RATIONAL );
     580      setReduce(a,false); // do not reduce expressions modulo mipo
     581      continue;
     582    }
     583    if( isLess(bound, other, 1, mv) ) // current prime unlucky
     584      continue;
     585    // here: isLess(other, bound, 1, mv) ) ==> all previous primes unlucky
     586    q = p;
     587    D = replacevar( mapinto(Dp), a, b ); // shortcut CRA // shortcut CRA
     588    for(int i=1; i<=mv; i++) // tighten bound
     589      bound[i] = other[i];
     590  }
     591  // hopefully, we never reach this point
     592  setReduce(a,true);
     593  Off( SW_USE_QGCD );
     594  D = gcd( f, g );
     595  On( SW_USE_QGCD );
     596  return D;
     597}
     598
     599
     600int * leadDeg(const CanonicalForm & f, int *degs)
     601{ // leading degree vector w.r.t. lex. monomial order x(i+1) > x(i)
     602  // if f is in a coeff domain, the zero pointer is returned
     603  // 'a' should point to an array of sufficient size level(f)+1
     604  if(f.inCoeffDomain())
     605    return 0;
     606  CanonicalForm tmp = f;
     607  do
     608  {
     609    degs[tmp.level()] = tmp.degree();
     610    tmp = LC(tmp);
     611  }
     612  while(!tmp.inCoeffDomain());
     613  return degs;
     614}
     615
     616
     617bool isLess(int *a, int *b, int lower, int upper)
     618{ // compares the degree vectors a,b on the specified part. Note: x(i+1) > x(i)
     619  for(int i=upper; i>=lower; i--)
     620    if(a[i] == b[i])
     621      continue;
     622    else
     623      return a[i] < b[i];
     624  return true;
     625}
     626
     627
     628bool isEqual(int *a, int *b, int lower, int upper)
     629{ // compares the degree vectors a,b on the specified part. Note: x(i+1) > x(i)
     630  for(int i=lower; i<=upper; i++)
     631    if(a[i] != b[i])
     632      return false;
     633  return true;
     634}
     635
     636
     637CanonicalForm firstLC(const CanonicalForm & f)
     638{ // returns the leading coefficient (LC) of level <= 1
     639  CanonicalForm ret = f;
     640  while(ret.level() > 1)
     641    ret = LC(ret);
     642  return ret;
     643}
     644
     645
     646void tryCRA( const CanonicalForm & x1, const CanonicalForm & q1, const CanonicalForm & x2, const CanonicalForm & q2, const CanonicalForm & M, CanonicalForm & xnew, CanonicalForm & qnew, bool & fail )
     647{ // as CRA, but takes care of zero divisors
     648  CanonicalForm tmp;
     649  if(x1.level() <= 1 && x2.level() <= 1) // base case
     650  {
     651    tryExtgcd(q1,q2,M,tmp,xnew,qnew,fail);
     652    if(fail)
     653      return;
     654    xnew = x1 + (x2-x1) * xnew * q1;
     655    qnew = q1*q2;
     656    xnew = mod(xnew,qnew);
     657    return;
     658  }
     659  CanonicalForm tmp2;
     660  xnew = 0;
     661  qnew = q1 * q2;
     662  // here: x1.level() > 1 || x2.level() > 1
     663  if(x1.level() > x2.level())
     664  {
     665    for(CFIterator i=x1; i.hasTerms(); i++)
     666    {
     667      if(i.exp() == 0) // const. term
     668      {
     669        tryCRA(i.coeff(),q1,x2,q2,M,tmp,tmp2,fail);
     670        if(fail)
     671          return;
     672        xnew += tmp;
     673      }
     674      else
     675      {
     676        tryCRA(i.coeff(),q1,0,q2,M,tmp,tmp2,fail);
     677        if(fail)
     678          return;
     679        xnew += tmp * power(x1.mvar(),i.exp());
     680      }
     681    }
     682    return;
     683  }
     684  // here: x1.level() <= x2.level() && ( x1.level() > 1 || x2.level() > 1 )
     685  if(x2.level() > x1.level())
     686  {
     687    for(CFIterator j=x2; j.hasTerms(); j++)
     688    {
     689      if(j.exp() == 0) // const. term
     690      {
     691        tryCRA(x1,q1,j.coeff(),q2,M,tmp,tmp2,fail);
     692        if(fail)
     693          return;
     694        xnew += tmp;
     695      }
     696      else
     697      {
     698        tryCRA(0,q1,j.coeff(),q2,M,tmp,tmp2,fail);
     699        if(fail)
     700          return;
     701        xnew += tmp * power(x2.mvar(),j.exp());
     702      }
     703    }
     704    return;
     705  }
     706  // here: x1.level() == x2.level() && x1.level() > 1 && x2.level() > 1
     707  CFIterator i = x1;
     708  CFIterator j = x2;
     709  while(i.hasTerms() || j.hasTerms())
     710  {
     711    if(i.hasTerms())
     712    {
     713      if(j.hasTerms())
     714      {
     715        if(i.exp() == j.exp())
     716        {
     717          tryCRA(i.coeff(),q1,j.coeff(),q2,M,tmp,tmp2,fail);
     718          if(fail)
     719            return;
     720          xnew += tmp * power(x1.mvar(),i.exp());
     721          i++; j++;
     722        }
     723        else
     724        {
     725          if(i.exp() < j.exp())
     726          {
     727            tryCRA(i.coeff(),q1,0,q2,M,tmp,tmp2,fail);
     728            if(fail)
     729              return;
     730            xnew += tmp * power(x1.mvar(),i.exp());
     731            i++;
     732          }
     733          else // i.exp() > j.exp()
     734          {
     735            tryCRA(0,q1,j.coeff(),q2,M,tmp,tmp2,fail);
     736            if(fail)
     737              return;
     738            xnew += tmp * power(x1.mvar(),j.exp());
     739            j++;
     740          }
     741        }
     742      }
     743      else // j is out of terms
     744      {
     745        tryCRA(i.coeff(),q1,0,q2,M,tmp,tmp2,fail);
     746        if(fail)
     747          return;
     748        xnew += tmp * power(x1.mvar(),i.exp());
     749        i++;
     750      }
     751    }
     752    else // i is out of terms
     753    {
     754      tryCRA(0,q1,j.coeff(),q2,M,tmp,tmp2,fail);
     755      if(fail)
     756        return;
     757      xnew += tmp * power(x1.mvar(),j.exp());
     758      j++;
     759    }
     760  }
     761}
     762
     763
     764void tryExtgcd( const CanonicalForm & F, const CanonicalForm & G, const CanonicalForm & M, CanonicalForm & result, CanonicalForm & s, CanonicalForm & t, bool & fail )
     765{ // F, G are univariate polynomials (i.e. they have exactly one polynomial variable)
     766  // F and G must have the same level AND level > 0
     767  // we try to calculate gcd(F,G) = s*F + t*G
     768  // if a zero divisor is encontered, 'fail' is set to one
     769  // M is assumed to be monic
     770  CanonicalForm P;
     771  if(F.inCoeffDomain())
     772  {
     773    tryInvert( F, M, P, fail );
     774    if(fail)
     775      return;
     776    result = 1;
     777    s = P; t = 0;
     778    return;
     779  }
     780  if(G.inCoeffDomain())
     781  {
     782    tryInvert( G, M, P, fail );
     783    if(fail)
     784      return;
     785    result = 1;
     786    s = 0; t = P;
     787    return;
     788  }
     789  // here: both not inCoeffDomain
     790  CanonicalForm inv, rem, tmp, u, v, q, sum=0;
     791  if( F.degree() > G.degree() )
     792  {
     793    P = F; result = G;  s=v=0; t=u=1;
     794  }
     795  else
     796  {
     797    P = G; result = F; s=v=1; t=u=0;
     798  }
     799  Variable x = P.mvar();
     800  // here: degree(P) >= degree(result)
     801  while(true)
     802  {
     803    tryInvert( Lc(result), M, inv, fail );
     804    if(fail)
     805      return;
     806    // here: Lc(result) is invertible modulo M
     807    q = Lc(P)*inv*power( x, P.degree(x)-result.degree(x) );
     808    rem = reduce( P - q*result, M );
     809    if( rem.isZero() )
     810    {
     811      s*=inv;
     812      t*=inv;
     813      result *= inv; // monify result
     814      return;
     815    }
     816    sum += q;
     817    if(result.degree(x) >= rem.degree(x))
     818    {
     819      P=result;
     820      result=rem;
     821      tmp=u-sum*s;
     822      u=s;
     823      s=tmp;
     824      tmp=v-sum*t;
     825      v=t;
     826      t=tmp;
     827      sum = 0; // reset
     828    }
     829    else
     830      P = rem;
     831  }
     832}
     833
     834
     835static CanonicalForm trycontent ( const CanonicalForm & f, const Variable & x, const CanonicalForm & M, bool & fail )
     836{ // as 'content', but takes care of zero divisors
     837  ASSERT( x.level() > 0, "cannot calculate content with respect to algebraic variable" );
     838  Variable y = f.mvar();
     839  if ( y == x )
     840    return trycf_content( f, 0, M, fail );
     841  if ( y < x )
     842     return f;
     843  return swapvar( trycontent( swapvar( f, y, x ), y, M, fail ), y, x );
     844}
     845
     846
     847static CanonicalForm tryvcontent ( const CanonicalForm & f, const Variable & x, const CanonicalForm & M, bool & fail )
     848{ // as vcontent, but takes care of zero divisors
     849  ASSERT( x.level() > 0, "cannot calculate vcontent with respect to algebraic variable" );
     850  if ( f.mvar() <= x )
     851    return trycontent( f, x, M, fail );
     852  CFIterator i;
     853  CanonicalForm d = 0, e, ret;
     854  for ( i = f; i.hasTerms() && ! d.isOne() && ! fail; i++ )
     855  {
     856    e = tryvcontent( i.coeff(), x, M, fail );
     857    if(fail)
     858      break;
     859    tryBrownGCD( d, e, M, ret, fail );
     860    d = ret;
     861  }
     862  return d;
     863}
     864
     865
     866static CanonicalForm trycf_content ( const CanonicalForm & f, const CanonicalForm & g, const CanonicalForm & M, bool & fail )
     867{ // as cf_content, but takes care of zero divisors
     868  if ( f.inPolyDomain() || ( f.inExtension() && ! getReduce( f.mvar() ) ) )
     869  {
     870    CFIterator i = f;
     871    CanonicalForm tmp = g, result;
     872    while ( i.hasTerms() && ! tmp.isOne() && ! fail )
     873    {
     874      tryBrownGCD( i.coeff(), tmp, M, result, fail );
     875      tmp = result;
     876      i++;
     877    }
     878    return result;
     879  }
     880  return abs( f );
     881}
     882
     883
     884static void tryDivide( const CanonicalForm & f, const CanonicalForm & g, const CanonicalForm & M, CanonicalForm & result, bool & divides, bool & fail )
     885{ // M "univariate" monic polynomial
     886  // f, g polynomials with coeffs modulo M.
     887  // if f is divisible by g, 'divides' is set to 1 and 'result' == f/g mod M coefficientwise.
     888  // 'fail' is set to 1, iff a zero divisor is encountered.
     889  // divides==1 implies fail==0
     890  // required: getReduce(M.mvar())==0
     891  if(g.inBaseDomain())
     892  {
     893    result = f/g;
     894    divides = true;
     895    return;
     896  }
     897  if(g.inCoeffDomain())
     898  {
     899    tryInvert(g,M,result,fail);
     900    if(fail)
     901      return;
     902    result = reduce(f*result, M);
     903    divides = true;
     904    return;
     905  }
     906  // here: g NOT inCoeffDomain
     907  Variable x = g.mvar();
     908  if(f.degree(x) < g.degree(x))
     909  {
     910    divides = false;
     911    return;
     912  }
     913  // here: f.degree(x) > 0 and f.degree(x) >= g.degree(x)
     914  CanonicalForm F = f;
     915  CanonicalForm q, leadG = LC(g);
     916  result = 0;
     917  while(!F.isZero())
     918  {
     919    tryDivide(F.LC(x),leadG,M,q,divides,fail);
     920    if(fail || !divides)
     921      return;
     922    if(F.degree(x)<g.degree(x))
     923    {
     924      divides = false;
     925      return;
     926    }
     927    q *= power(x,F.degree(x)-g.degree(x));
     928    result += q;
     929    F = reduce(F-q*g, M);
     930  }
     931  result = reduce(result, M);
     932  divides = true;
     933}
     934
    483935
    484936
Note: See TracChangeset for help on using the changeset viewer.