Changeset 5df7d0 in git for factory/algext.cc


Ignore:
Timestamp:
Jul 1, 2011, 12:02:02 PM (13 years ago)
Author:
Martin Lee <martinlee84@…>
Branches:
(u'spielwiese', 'fe61d9c35bf7c61f2b6cbf1b56e25e2f08d536cc')
Children:
33e850a665f1c7a320edcde4faacdd7602f49b76
Parents:
6f08f33cc3e3c2a69b00521f71d91ff88066e3c0
Message:
replaced tryCRA by tryNewtonInterp and deleted obsolete tryCRA
reduce result in tryBrownGCD modulo M
make modular image really monic (i.e. invert firstLC (g_image) instead of lc (g_image))
put reduce into header


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

Legend:

Unmodified
Added
Removed

  • factory/algext.cc

    r6f08f3 r5df7d0  
    314314      return;
    315315    result = inv*P; // monify result (not reduced, yet)
     316    result= reduce (result, M);
    316317    return;
    317318  }
    … …  
    327328    {
    328329      result *= inv;
     330      result= reduce (result, M);
    329331      return;
    330332    }
    … …  
    359361bool isEqual(int *a, int *b, int lower, int upper);
    360362CanonicalForm firstLC(const CanonicalForm & f);
    361 void tryCRA( const CanonicalForm & x1, const CanonicalForm & q1, const CanonicalForm & x2, const CanonicalForm & q2, const CanonicalForm & M, CanonicalForm & xnew, CanonicalForm & qnew, bool & fail );
    362 void tryExtgcd( const CanonicalForm & F, const CanonicalForm & G, const CanonicalForm & M, CanonicalForm & result, CanonicalForm & s, CanonicalForm & t, bool & fail );
    363363static CanonicalForm trycontent ( const CanonicalForm & f, const Variable & x, const CanonicalForm & M, bool & fail );
    364364static CanonicalForm tryvcontent ( const CanonicalForm & f, const Variable & x, const CanonicalForm & M, bool & fail );
    … …  
    366366static void tryDivide( const CanonicalForm & f, const CanonicalForm & g, const CanonicalForm & M, CanonicalForm & result, bool & divides, bool & fail );
    367367
     368static inline CanonicalForm
     369tryNewtonInterp (const CanonicalForm alpha, const CanonicalForm u,
     370              const CanonicalForm newtonPoly, const CanonicalForm oldInterPoly,
     371              const Variable & x, const CanonicalForm& M, bool& fail)
     372{
     373  CanonicalForm interPoly;
     374
     375  CanonicalForm inv;
     376  tryInvert (newtonPoly (alpha, x), M, inv, fail);
     377  if (fail)
     378    return 0;
     379
     380  interPoly= oldInterPoly+reduce ((u - oldInterPoly (alpha, x))*inv*newtonPoly, M);
     381  return interPoly;
     382}
    368383
    369384void tryBrownGCD( const CanonicalForm & F, const CanonicalForm & G, const CanonicalForm & M, CanonicalForm & result, bool & fail, bool topLevel )
    … …  
    391406      return;
    392407    result = inv*G;
     408    result= reduce (result, M);
    393409    return;
    394410  }
    … …  
    409425      return;
    410426    result = inv*F;
     427    result= reduce (result, M);
    411428    return;
    412429  }
    … …  
    448465    if(fail)
    449466      return;
    450     result = NN(result); // do not forget to map back
     467    result= NN (reduce (result, M)); // do not forget to map back
    451468    return;
    452469  }
    … …  
    509526  CanonicalForm gamma_image, m=1;
    510527  CanonicalForm gm=0;
    511   CanonicalForm g_image, alpha, gnew, mnew;
     528  CanonicalForm g_image, alpha, gnew;
    512529  FFGenerator gen = FFGenerator();
    513530  Variable x= Variable (1);
    … …  
    535552    if(isEqual(dg_im, L, 2, mv))
    536553    {
    537       g_image /= lc(g_image); // make g_image monic over Z/p
     554      CanonicalForm inv;
     555      tryInvert (firstLC (g_image), M, inv, fail);
     556      if (fail)
     557        return;
     558      g_image *= inv;
    538559      g_image *= gamma_image; // multiply by multiple of image lc(gcd)
    539       tryCRA( g_image, x-alpha, gm, m, M, gnew, mnew, fail );
     560      g_image= reduce (g_image, M);
     561      gnew= tryNewtonInterp (alpha, g_image, m, gm, x, M, fail);
    540562      // gnew = gm mod m
    541563      // gnew = g_image mod var(1)-alpha
    … …  
    543565      if(fail)
    544566        return;
    545       m = mnew;
     567      m *= (x - alpha);
    546568      if(gnew == gm) // gnew did not change
    547569      {
    … …  
    563585          if(divides)
    564586          {
    565             result = NN(c*g_image);
     587            result = NN(reduce (c*g_image, M));
    566588            return;
    567589          }
    … …  
    672694    // here: mipo is monic
    673695    tryBrownGCD( mapinto(f), mapinto(g), mipo, Dp, fail );
    674     Dp = reduce( Dp, mipo );
    675696    setCharacteristic(0);
    676697    if( fail ) // mipo splits in char p
    … …  
    781802}
    782803
    783 
    784 void tryCRA( const CanonicalForm & x1, const CanonicalForm & q1, const CanonicalForm & x2, const CanonicalForm & q2, const CanonicalForm & M, CanonicalForm & xnew, CanonicalForm & qnew, bool & fail )
    785 { // as CRA, but takes care of zero divisors
    786   CanonicalForm tmp;
    787   if(x1.level() <= 1 && x2.level() <= 1) // base case
    788   {
    789     tryExtgcd(q1,q2,M,tmp,xnew,qnew,fail);
    790     if(fail)
    791       return;
    792     xnew = x1 + (x2-x1) * xnew * q1;
    793     qnew = q1*q2;
    794     xnew = mod(xnew,qnew);
    795     return;
    796   }
    797   CanonicalForm tmp2;
    798   xnew = 0;
    799   qnew = q1 * q2;
    800   // here: x1.level() > 1 || x2.level() > 1
    801   if(x1.level() > x2.level())
    802   {
    803     for(CFIterator i=x1; i.hasTerms(); i++)
    804     {
    805       if(i.exp() == 0) // const. term
    806       {
    807         tryCRA(i.coeff(),q1,x2,q2,M,tmp,tmp2,fail);
    808         if(fail)
    809           return;
    810         xnew += tmp;
    811       }
    812       else
    813       {
    814         tryCRA(i.coeff(),q1,0,q2,M,tmp,tmp2,fail);
    815         if(fail)
    816           return;
    817         xnew += tmp * power(x1.mvar(),i.exp());
    818       }
    819     }
    820     return;
    821   }
    822   // here: x1.level() <= x2.level() && ( x1.level() > 1 || x2.level() > 1 )
    823   if(x2.level() > x1.level())
    824   {
    825     for(CFIterator j=x2; j.hasTerms(); j++)
    826     {
    827       if(j.exp() == 0) // const. term
    828       {
    829         tryCRA(x1,q1,j.coeff(),q2,M,tmp,tmp2,fail);
    830         if(fail)
    831           return;
    832         xnew += tmp;
    833       }
    834       else
    835       {
    836         tryCRA(0,q1,j.coeff(),q2,M,tmp,tmp2,fail);
    837         if(fail)
    838           return;
    839         xnew += tmp * power(x2.mvar(),j.exp());
    840       }
    841     }
    842     return;
    843   }
    844   // here: x1.level() == x2.level() && x1.level() > 1 && x2.level() > 1
    845   CFIterator i = x1;
    846   CFIterator j = x2;
    847   while(i.hasTerms() || j.hasTerms())
    848   {
    849     if(i.hasTerms())
    850     {
    851       if(j.hasTerms())
    852       {
    853         if(i.exp() == j.exp())
    854         {
    855           tryCRA(i.coeff(),q1,j.coeff(),q2,M,tmp,tmp2,fail);
    856           if(fail)
    857             return;
    858           xnew += tmp * power(x1.mvar(),i.exp());
    859           i++; j++;
    860         }
    861         else
    862         {
    863           if(i.exp() < j.exp())
    864           {
    865             tryCRA(i.coeff(),q1,0,q2,M,tmp,tmp2,fail);
    866             if(fail)
    867               return;
    868             xnew += tmp * power(x1.mvar(),i.exp());
    869             i++;
    870           }
    871           else // i.exp() > j.exp()
    872           {
    873             tryCRA(0,q1,j.coeff(),q2,M,tmp,tmp2,fail);
    874             if(fail)
    875               return;
    876             xnew += tmp * power(x1.mvar(),j.exp());
    877             j++;
    878           }
    879         }
    880       }
    881       else // j is out of terms
    882       {
    883         tryCRA(i.coeff(),q1,0,q2,M,tmp,tmp2,fail);
    884         if(fail)
    885           return;
    886         xnew += tmp * power(x1.mvar(),i.exp());
    887         i++;
    888       }
    889     }
    890     else // i is out of terms
    891     {
    892       tryCRA(0,q1,j.coeff(),q2,M,tmp,tmp2,fail);
    893       if(fail)
    894         return;
    895       xnew += tmp * power(x1.mvar(),j.exp());
    896       j++;
    897     }
    898   }
    899 }
    900 
    901 
    902804void tryExtgcd( const CanonicalForm & F, const CanonicalForm & G, const CanonicalForm & M, CanonicalForm & result, CanonicalForm & s, CanonicalForm & t, bool & fail )
    903805{ // F, G are univariate polynomials (i.e. they have exactly one polynomial variable)
    … …  
    1067969  divides = true;
    1068970}
    1069 
    1070 
    1071971
    1072972void tryCRA( const CanonicalForm & x1, const CanonicalForm & q1, const CanonicalForm & x2, const CanonicalForm & q2, CanonicalForm & xnew, CanonicalForm & qnew, bool & fail )
Note: See TracChangeset for help on using the changeset viewer.