Changeset f1e3bf5 in git for factory


Ignore:
Timestamp:
Jun 20, 2014, 12:25:12 PM (10 years ago)
Author:
Martin Lee <martinlee84@…>
Branches:
(u'spielwiese', 'fe61d9c35bf7c61f2b6cbf1b56e25e2f08d536cc')
Children:
5b4953b0435d257dfff1b9b6848b8c2a28e0f347
Parents:
baa2c3a8e2d77600e598268aca866716d7070466
git-author:
Martin Lee <martinlee84@web.de>2014-06-20 12:25:12+02:00
git-committer:
Martin Lee <martinlee84@web.de>2014-06-24 12:06:05+02:00
Message:
chg: moved EZGCD_P from cf_gcd_smallp to cfEzgcd
Location:
factory
Files:
1 added
4 edited

Legend:

Unmodified
Added
Removed

  • factory/cfEzgcd.cc

    rbaa2c3a rf1e3bf5  
    22 * Computer Algebra System SINGULAR
    33\*****************************************************************************/
    4 /** @file fac_ezgcd.cc
     4/** @file cfEzgcd.cc
    55 *
    6  * This file implements the GCD of two multivariate polynomials over Q using
    7  * EZ-GCD as described in "Algorithms for Computer Algebra" by Geddes, Czapor,
    8  * Labahnn
     6 * This file implements the GCD of two multivariate polynomials over Q or F_q
     7 * using EZ-GCD as described in "Algorithms for Computer Algebra" by Geddes,
     8 * Czapor, Labahnn
    99 *
    1010 * @author Martin Lee
    … …  
    2323#include "cf_defs.h"
    2424#include "canonicalform.h"
    25 #include "fac_util.h" // HEADER for this file
     25#include "cfEzgcd.h"
     26#include "cf_gcd_smallp.h"
     27#include "cf_util.h"
     28#include "cf_map_ext.h"
    2629#include "cf_algorithm.h"
    2730#include "cf_reval.h"
    … …  
    3336
    3437#ifdef HAVE_NTL
     38#include "NTLconvert.h"
     39
     40static const double log2exp= 1.442695041;
    3541
    3642TIMING_DEFINE_PRINT(ez_eval)
    … …  
    783789}
    784790
     791static inline
     792int Hensel_P (const CanonicalForm & UU, CFArray & G, const Evaluation & AA,
     793              const CFArray& LeadCoeffs )
     794{
     795  CFList factors;
     796  factors.append (G[1]);
     797  factors.append (G[2]);
     798
     799  CFMap NN, MM;
     800  Evaluation A= optimize4Lift (UU, MM, NN, AA);
     801
     802  CanonicalForm U= MM (UU);
     803  CFArray LCs= CFArray (1,2);
     804  LCs [1]= MM (LeadCoeffs [1]);
     805  LCs [2]= MM (LeadCoeffs [2]);
     806
     807  CFList evaluation;
     808  long termEstimate= size (U);
     809  for (int i= A.min(); i <= A.max(); i++)
     810  {
     811    if (!A[i].isZero() && (getCharacteristic() > degree (U,i))) //TODO find a good estimate for getCharacteristic() <= degree (U,i)
     812    {
     813      termEstimate *= degree (U,i)*2;
     814      termEstimate /= 3;
     815    }
     816    evaluation.append (A [i]);
     817  }
     818  if (termEstimate/getNumVars(U) > 500)
     819    return -1;
     820  CFList UEval;
     821  CanonicalForm shiftedU= myShift2Zero (U, UEval, evaluation);
     822
     823  if (size (shiftedU)/getNumVars (U) > 500)
     824    return -1;
     825
     826  CFArray shiftedLCs= CFArray (2);
     827  CFList shiftedLCsEval1, shiftedLCsEval2;
     828  shiftedLCs[0]= myShift2Zero (LCs[1], shiftedLCsEval1, evaluation);
     829  shiftedLCs[1]= myShift2Zero (LCs[2], shiftedLCsEval2, evaluation);
     830  factors.insert (1);
     831  int liftBound= degree (UEval.getLast(), 2) + 1;
     832  CFArray Pi;
     833  CFMatrix M= CFMatrix (liftBound, factors.length() - 1);
     834  CFList diophant;
     835  CFArray lcs= CFArray (2);
     836  lcs [0]= shiftedLCsEval1.getFirst();
     837  lcs [1]= shiftedLCsEval2.getFirst();
     838  nonMonicHenselLift12 (UEval.getFirst(), factors, liftBound, Pi, diophant, M,
     839                        lcs, false);
     840
     841  for (CFListIterator i= factors; i.hasItem(); i++)
     842  {
     843    if (!fdivides (i.getItem(), UEval.getFirst()))
     844      return 0;
     845  }
     846
     847  int * liftBounds;
     848  bool noOneToOne= false;
     849  if (U.level() > 2)
     850  {
     851    liftBounds= new int [U.level() - 1]; /* index: 0.. U.level()-2 */
     852    liftBounds[0]= liftBound;
     853    for (int i= 1; i < U.level() - 1; i++)
     854      liftBounds[i]= degree (shiftedU, Variable (i + 2)) + 1;
     855    factors= nonMonicHenselLift2 (UEval, factors, liftBounds, U.level() - 1,
     856                                  false, shiftedLCsEval1, shiftedLCsEval2, Pi,
     857                                  diophant, noOneToOne);
     858    delete [] liftBounds;
     859    if (noOneToOne)
     860      return 0;
     861  }
     862  G[1]= factors.getFirst();
     863  G[2]= factors.getLast();
     864  G[1]= myReverseShift (G[1], evaluation);
     865  G[2]= myReverseShift (G[2], evaluation);
     866  G[1]= NN (G[1]);
     867  G[2]= NN (G[2]);
     868  return 1;
     869}
     870
     871static inline
     872bool findeval_P (const CanonicalForm & F, const CanonicalForm & G,
     873                 CanonicalForm & Fb, CanonicalForm & Gb, CanonicalForm & Db,
     874                 REvaluation & b, int delta, int degF, int degG, int maxeval,
     875                 int & count, int& k, int bound, int& l)
     876{
     877  if( count == 0 && delta != 0)
     878  {
     879    if( count++ > maxeval )
     880      return false;
     881  }
     882  if (count > 0)
     883  {
     884    b.nextpoint(k);
     885    if (k == 0)
     886      k++;
     887    l++;
     888    if (l > bound)
     889    {
     890      l= 1;
     891      k++;
     892      if (k > tmax (F.level(), G.level()) - 1)
     893        return false;
     894      b.nextpoint (k);
     895    }
     896    if (count++ > maxeval)
     897      return false;
     898  }
     899  while( true )
     900  {
     901    Fb = b( F );
     902    if( degree( Fb, 1 ) == degF )
     903    {
     904      Gb = b( G );
     905      if( degree( Gb, 1 ) == degG )
     906      {
     907        Db = gcd( Fb, Gb );
     908        if( delta > 0 )
     909        {
     910          if( degree( Db, 1 ) <= delta )
     911            return true;
     912        }
     913        else
     914          return true;
     915      }
     916    }
     917    if (k == 0)
     918      k++;
     919    b.nextpoint(k);
     920    l++;
     921    if (l > bound)
     922    {
     923      l= 1;
     924      k++;
     925      if (k > tmax (F.level(), G.level()) - 1)
     926        return false;
     927      b.nextpoint (k);
     928    }
     929    if( count++ > maxeval )
     930      return false;
     931  }
     932}
     933
     934// parameters for heuristic
     935static int maxNumEval= 200;
     936static int sizePerVars1= 500; //try dense gcd if size/#variables is bigger
     937
     938/// Extended Zassenhaus GCD for finite fields
     939CanonicalForm EZGCD_P( const CanonicalForm & FF, const CanonicalForm & GG )
     940{
     941  if (FF.isZero() && degree(GG) > 0) return GG/Lc(GG);
     942  else if (GG.isZero() && degree (FF) > 0) return FF/Lc(FF);
     943  else if (FF.isZero() && GG.isZero()) return FF.genOne();
     944  if (FF.inBaseDomain() || GG.inBaseDomain()) return FF.genOne();
     945  if (FF.isUnivariate() && fdivides(FF, GG)) return FF/Lc(FF);
     946  if (GG.isUnivariate() && fdivides(GG, FF)) return GG/Lc(GG);
     947  if (FF == GG) return FF/Lc(FF);
     948
     949  int maxNumVars= tmax (getNumVars (FF), getNumVars (GG));
     950  Variable a, oldA;
     951  int sizeF= size (FF);
     952  int sizeG= size (GG);
     953
     954  if (sizeF/maxNumVars > sizePerVars1 && sizeG/maxNumVars > sizePerVars1)
     955  {
     956    if (hasFirstAlgVar (FF, a) || hasFirstAlgVar (GG, a))
     957      return GCD_Fp_extension (FF, GG, a);
     958    else if (CFFactory::gettype() == GaloisFieldDomain)
     959      return GCD_GF (FF, GG);
     960    else
     961      return GCD_small_p (FF, GG);
     962  }
     963
     964  CanonicalForm F, G, f, g, d, Fb, Gb, Db, Fbt, Gbt, Dbt, B0, B, D0, lcF, lcG,
     965                lcD;
     966  CFArray DD( 1, 2 ), lcDD( 1, 2 );
     967  int degF, degG, delta, count;
     968  int maxeval;
     969  maxeval= tmin((getCharacteristic()/
     970                (int)(ilog2(getCharacteristic())*log2exp))*2, maxNumEval);
     971  count= 0; // number of eval. used
     972  REvaluation b, bt;
     973  int gcdfound = 0;
     974  Variable x = Variable(1);
     975
     976  F= FF;
     977  G= GG;
     978
     979  CFMap M,N;
     980  int smallestDegLev;
     981  TIMING_START (ez_p_compress)
     982  int best_level= compress4EZGCD (F, G, M, N, smallestDegLev);
     983
     984  if (best_level == 0) return G.genOne();
     985
     986  F= M (F);
     987  G= M (G);
     988  TIMING_END_AND_PRINT (ez_p_compress, "time for compression in EZ_P: ")
     989
     990  TIMING_START (ez_p_content)
     991  f = content( F, x ); g = content( G, x );
     992  d = gcd( f, g );
     993  F /= f; G /= g;
     994  TIMING_END_AND_PRINT (ez_p_content, "time to extract content in EZ_P: ")
     995
     996  if( F.isUnivariate() && G.isUnivariate() )
     997  {
     998    if( F.mvar() == G.mvar() )
     999      d *= gcd( F, G );
     1000    else
     1001      return N (d);
     1002    return N (d);
     1003  }
     1004  if ( F.isUnivariate())
     1005  {
     1006    g= content (G,G.mvar());
     1007    return N(d*gcd(F,g));
     1008  }
     1009  if ( G.isUnivariate())
     1010  {
     1011    f= content (F,F.mvar());
     1012    return N(d*gcd(G,f));
     1013  }
     1014
     1015  maxNumVars= tmax (getNumVars (F), getNumVars (G));
     1016  sizeF= size (F);
     1017  sizeG= size (G);
     1018
     1019  if (sizeF/maxNumVars > sizePerVars1 && sizeG/maxNumVars > sizePerVars1)
     1020  {
     1021    if (hasFirstAlgVar (F, a) || hasFirstAlgVar (G, a))
     1022      return N (d*GCD_Fp_extension (F, G, a));
     1023    else if (CFFactory::gettype() == GaloisFieldDomain)
     1024      return N (d*GCD_GF (F, G));
     1025    else
     1026      return N (d*GCD_small_p (F, G));
     1027  }
     1028
     1029  int dummy= 0;
     1030  if( gcd_test_one( F, G, false, dummy ) )
     1031  {
     1032    return N (d);
     1033  }
     1034
     1035  bool passToGF= false;
     1036  bool extOfExt= false;
     1037  int p= getCharacteristic();
     1038  bool algExtension= (hasFirstAlgVar(F,a) || hasFirstAlgVar(G,a));
     1039  int k= 1;
     1040  CanonicalForm primElem, imPrimElem;
     1041  CFList source, dest;
     1042  if (p < 50 && CFFactory::gettype() != GaloisFieldDomain && !algExtension)
     1043  {
     1044    if (p == 2)
     1045      setCharacteristic (2, 12, 'Z');
     1046    else if (p == 3)
     1047      setCharacteristic (3, 4, 'Z');
     1048    else if (p == 5 || p == 7)
     1049      setCharacteristic (p, 3, 'Z');
     1050    else
     1051      setCharacteristic (p, 2, 'Z');
     1052    passToGF= true;
     1053    F= F.mapinto();
     1054    G= G.mapinto();
     1055    maxeval= 2*ipower (p, getGFDegree());
     1056  }
     1057  else if (CFFactory::gettype() == GaloisFieldDomain &&
     1058           ipower (p , getGFDegree()) < 50)
     1059  {
     1060    k= getGFDegree();
     1061    if (ipower (p, 2*k) > 50)
     1062      setCharacteristic (p, 2*k, gf_name);
     1063    else
     1064      setCharacteristic (p, 3*k, gf_name);
     1065    F= GFMapUp (F, k);
     1066    G= GFMapUp (G, k);
     1067    maxeval= tmin (2*ipower (p, getGFDegree()), maxNumEval);
     1068  }
     1069  else if (p < 50 && algExtension && CFFactory::gettype() != GaloisFieldDomain)
     1070  {
     1071    int d= degree (getMipo (a));
     1072    oldA= a;
     1073    Variable v2;
     1074    if (p == 2 && d < 6)
     1075    {
     1076      if (fac_NTL_char != p)
     1077      {
     1078        fac_NTL_char= p;
     1079        zz_p::init (p);
     1080      }
     1081      bool primFail= false;
     1082      Variable vBuf;
     1083      primElem= primitiveElement (a, vBuf, primFail);
     1084      ASSERT (!primFail, "failure in integer factorizer");
     1085      if (d < 3)
     1086      {
     1087        zz_pX NTLIrredpoly;
     1088        BuildIrred (NTLIrredpoly, d*3);
     1089        CanonicalForm newMipo= convertNTLzzpX2CF (NTLIrredpoly, Variable (1));
     1090        v2= rootOf (newMipo);
     1091      }
     1092      else
     1093      {
     1094        zz_pX NTLIrredpoly;
     1095        BuildIrred (NTLIrredpoly, d*2);
     1096        CanonicalForm newMipo= convertNTLzzpX2CF (NTLIrredpoly, Variable (1));
     1097        v2= rootOf (newMipo);
     1098      }
     1099      imPrimElem= mapPrimElem (primElem, a, v2);
     1100      extOfExt= true;
     1101    }
     1102    else if ((p == 3 && d < 4) || ((p == 5 || p == 7) && d < 3))
     1103    {
     1104      if (fac_NTL_char != p)
     1105      {
     1106        fac_NTL_char= p;
     1107        zz_p::init (p);
     1108      }
     1109      bool primFail= false;
     1110      Variable vBuf;
     1111      primElem= primitiveElement (a, vBuf, primFail);
     1112      ASSERT (!primFail, "failure in integer factorizer");
     1113      zz_pX NTLIrredpoly;
     1114      BuildIrred (NTLIrredpoly, d*2);
     1115      CanonicalForm newMipo= convertNTLzzpX2CF (NTLIrredpoly, Variable (1));
     1116      v2= rootOf (newMipo);
     1117      imPrimElem= mapPrimElem (primElem, a, v2);
     1118      extOfExt= true;
     1119    }
     1120    if (extOfExt)
     1121    {
     1122      maxeval= tmin (2*ipower (p, degree (getMipo (v2))), maxNumEval);
     1123      F= mapUp (F, a, v2, primElem, imPrimElem, source, dest);
     1124      G= mapUp (G, a, v2, primElem, imPrimElem, source, dest);
     1125      a= v2;
     1126    }
     1127  }
     1128
     1129  lcF = LC( F, x ); lcG = LC( G, x );
     1130  lcD = gcd( lcF, lcG );
     1131
     1132  delta = 0;
     1133  degF = degree( F, x ); degG = degree( G, x );
     1134
     1135  if(hasFirstAlgVar(G,a))
     1136    b = REvaluation( 2, tmax(F.level(), G.level()), AlgExtRandomF( a ) );
     1137  else
     1138  { // both not in extension given by algebraic variable
     1139    if (CFFactory::gettype() != GaloisFieldDomain)
     1140      b = REvaluation( 2, tmax(F.level(), G.level()), FFRandom() );
     1141    else
     1142      b = REvaluation( 2, tmax(F.level(), G.level()), GFRandom() );
     1143  }
     1144
     1145  CanonicalForm cand, contcand;
     1146  CanonicalForm result;
     1147  int o, t;
     1148  o= 0;
     1149  t= 1;
     1150  int goodPointCount= 0;
     1151  while( !gcdfound )
     1152  {
     1153    TIMING_START (ez_p_eval);
     1154    if( !findeval_P( F, G, Fb, Gb, Db, b, delta, degF, degG, maxeval, count, o,
     1155         maxeval/maxNumVars, t ))
     1156    { // too many eval. used --> try another method
     1157      Off (SW_USE_EZGCD_P);
     1158      result= gcd (F,G);
     1159      On (SW_USE_EZGCD_P);
     1160      if (passToGF)
     1161      {
     1162        CanonicalForm mipo= gf_mipo;
     1163        setCharacteristic (p);
     1164        Variable alpha= rootOf (mipo.mapinto());
     1165        result= GF2FalphaRep (result, alpha);
     1166      }
     1167      if (k > 1)
     1168      {
     1169        result= GFMapDown (result, k);
     1170        setCharacteristic (p, k, gf_name);
     1171      }
     1172      if (extOfExt)
     1173        result= mapDown (result, primElem, imPrimElem, oldA, dest, source);
     1174      return N (d*result);
     1175    }
     1176    TIMING_END_AND_PRINT (ez_p_eval, "time for eval point search in EZ_P1: ");
     1177    delta = degree( Db );
     1178    if (delta == degF)
     1179    {
     1180      if (degF <= degG && fdivides (F, G))
     1181      {
     1182        if (passToGF)
     1183        {
     1184          CanonicalForm mipo= gf_mipo;
     1185          setCharacteristic (p);
     1186          Variable alpha= rootOf (mipo.mapinto());
     1187          F= GF2FalphaRep (F, alpha);
     1188        }
     1189        if (k > 1)
     1190        {
     1191          F= GFMapDown (F, k);
     1192          setCharacteristic (p, k, gf_name);
     1193        }
     1194        if (extOfExt)
     1195          F= mapDown (F, primElem, imPrimElem, oldA, dest, source);
     1196        return N (d*F);
     1197      }
     1198      else
     1199        delta--;
     1200    }
     1201    else if (delta == degG)
     1202    {
     1203      if (degG <= degF && fdivides (G, F))
     1204      {
     1205        if (passToGF)
     1206        {
     1207          CanonicalForm mipo= gf_mipo;
     1208          setCharacteristic (p);
     1209          Variable alpha= rootOf (mipo.mapinto());
     1210          G= GF2FalphaRep (G, alpha);
     1211        }
     1212        if (k > 1)
     1213        {
     1214          G= GFMapDown (G, k);
     1215          setCharacteristic (p, k, gf_name);
     1216        }
     1217        if (extOfExt)
     1218          G= mapDown (G, primElem, imPrimElem, oldA, dest, source);
     1219        return N (d*G);
     1220      }
     1221      else
     1222        delta--;
     1223    }
     1224    if( delta == 0 )
     1225    {
     1226      if (passToGF)
     1227        setCharacteristic (p);
     1228      if (k > 1)
     1229        setCharacteristic (p, k, gf_name);
     1230      return N (d);
     1231    }
     1232    while( true )
     1233    {
     1234      bt = b;
     1235      TIMING_START (ez_p_eval);
     1236      if( !findeval_P(F,G,Fbt,Gbt,Dbt, bt, delta, degF, degG, maxeval, count, o,
     1237           maxeval/maxNumVars, t ))
     1238      { // too many eval. used --> try another method
     1239        Off (SW_USE_EZGCD_P);
     1240        result= gcd (F,G);
     1241        On (SW_USE_EZGCD_P);
     1242        if (passToGF)
     1243        {
     1244          CanonicalForm mipo= gf_mipo;
     1245          setCharacteristic (p);
     1246          Variable alpha= rootOf (mipo.mapinto());
     1247          result= GF2FalphaRep (result, alpha);
     1248        }
     1249        if (k > 1)
     1250        {
     1251          result= GFMapDown (result, k);
     1252          setCharacteristic (p, k, gf_name);
     1253        }
     1254        if (extOfExt)
     1255          result= mapDown (result, primElem, imPrimElem, oldA, dest, source);
     1256        return N (d*result);
     1257      }
     1258      TIMING_END_AND_PRINT (ez_p_eval, "time for eval point search in EZ_P2: ");
     1259      int dd = degree( Dbt );
     1260      if( dd == 0 )
     1261      {
     1262        if (passToGF)
     1263          setCharacteristic (p);
     1264        if (k > 1)
     1265          setCharacteristic (p, k, gf_name);
     1266        return N (d);
     1267      }
     1268      if( dd == delta )
     1269      {
     1270        goodPointCount++;
     1271        if (goodPointCount == 5)
     1272          break;
     1273      }
     1274      if( dd < delta )
     1275      {
     1276        goodPointCount= 0;
     1277        delta = dd;
     1278        b = bt;
     1279        Db = Dbt; Fb = Fbt; Gb = Gbt;
     1280      }
     1281      if (delta == degF)
     1282      {
     1283        if (degF <= degG && fdivides (F, G))
     1284        {
     1285          if (passToGF)
     1286          {
     1287            CanonicalForm mipo= gf_mipo;
     1288            setCharacteristic (p);
     1289            Variable alpha= rootOf (mipo.mapinto());
     1290            F= GF2FalphaRep (F, alpha);
     1291          }
     1292          if (k > 1)
     1293          {
     1294            F= GFMapDown (F, k);
     1295            setCharacteristic (p, k, gf_name);
     1296          }
     1297          if (extOfExt)
     1298            F= mapDown (F, primElem, imPrimElem, oldA, dest, source);
     1299          return N (d*F);
     1300        }
     1301        else
     1302          delta--;
     1303      }
     1304      else if (delta == degG)
     1305      {
     1306        if (degG <= degF && fdivides (G, F))
     1307        {
     1308          if (passToGF)
     1309          {
     1310            CanonicalForm mipo= gf_mipo;
     1311            setCharacteristic (p);
     1312            Variable alpha= rootOf (mipo.mapinto());
     1313            G= GF2FalphaRep (G, alpha);
     1314          }
     1315          if (k > 1)
     1316          {
     1317            G= GFMapDown (G, k);
     1318            setCharacteristic (p, k, gf_name);
     1319          }
     1320          if (extOfExt)
     1321            G= mapDown (G, primElem, imPrimElem, oldA, dest, source);
     1322          return N (d*G);
     1323        }
     1324        else
     1325          delta--;
     1326      }
     1327      if( delta == 0 )
     1328      {
     1329        if (passToGF)
     1330          setCharacteristic (p);
     1331        if (k > 1)
     1332          setCharacteristic (p, k, gf_name);
     1333        return N (d);
     1334      }
     1335    }
     1336    if( delta != degF && delta != degG )
     1337    {
     1338      bool B_is_F;
     1339      CanonicalForm xxx1, xxx2;
     1340      CanonicalForm buf;
     1341      DD[1] = Fb / Db;
     1342      buf= Gb/Db;
     1343      xxx1 = gcd( DD[1], Db );
     1344      xxx2 = gcd( buf, Db );
     1345      if (((xxx1.inCoeffDomain() && xxx2.inCoeffDomain()) &&
     1346          (size (F) <= size (G)))
     1347          || (xxx1.inCoeffDomain() && !xxx2.inCoeffDomain()))
     1348      {
     1349        B = F;
     1350        DD[2] = Db;
     1351        lcDD[1] = lcF;
     1352        lcDD[2] = lcD;
     1353        B_is_F = true;
     1354      }
     1355      else if (((xxx1.inCoeffDomain() && xxx2.inCoeffDomain()) &&
     1356               (size (G) < size (F)))
     1357               || (!xxx1.inCoeffDomain() && xxx2.inCoeffDomain()))
     1358      {
     1359        DD[1] = buf;
     1360        B = G;
     1361        DD[2] = Db;
     1362        lcDD[1] = lcG;
     1363        lcDD[2] = lcD;
     1364        B_is_F = false;
     1365      }
     1366      else // special case handling
     1367      {
     1368        Off (SW_USE_EZGCD_P);
     1369        result= gcd (F,G);
     1370        On (SW_USE_EZGCD_P);
     1371        if (passToGF)
     1372        {
     1373          CanonicalForm mipo= gf_mipo;
     1374          setCharacteristic (p);
     1375          Variable alpha= rootOf (mipo.mapinto());
     1376          result= GF2FalphaRep (result, alpha);
     1377        }
     1378        if (k > 1)
     1379        {
     1380          result= GFMapDown (result, k);
     1381          setCharacteristic (p, k, gf_name);
     1382        }
     1383        if (extOfExt)
     1384          result= mapDown (result, primElem, imPrimElem, oldA, dest, source);
     1385        return N (d*result);
     1386      }
     1387      DD[2] = DD[2] * ( b( lcDD[2] ) / lc( DD[2] ) );
     1388      DD[1] = DD[1] * ( b( lcDD[1] ) / lc( DD[1] ) );
     1389
     1390      if (size (B*lcDD[2])/maxNumVars > sizePerVars1)
     1391      {
     1392        if (algExtension)
     1393        {
     1394          result= GCD_Fp_extension (F, G, a);
     1395          if (extOfExt)
     1396            result= mapDown (result, primElem, imPrimElem, oldA, dest, source);
     1397          return N (d*result);
     1398        }
     1399        if (CFFactory::gettype() == GaloisFieldDomain)
     1400        {
     1401          result= GCD_GF (F, G);
     1402          if (passToGF)
     1403          {
     1404            CanonicalForm mipo= gf_mipo;
     1405            setCharacteristic (p);
     1406            Variable alpha= rootOf (mipo.mapinto());
     1407            result= GF2FalphaRep (result, alpha);
     1408          }
     1409          if (k > 1)
     1410          {
     1411            result= GFMapDown (result, k);
     1412            setCharacteristic (p, k, gf_name);
     1413          }
     1414          return N (d*result);
     1415        }
     1416        else
     1417          return N (d*GCD_small_p (F,G));
     1418      }
     1419
     1420      TIMING_START (ez_p_hensel_lift);
     1421      gcdfound= Hensel_P (B*lcD, DD, b, lcDD);
     1422      TIMING_END_AND_PRINT (ez_p_hensel_lift, "time for Hensel lift in EZ_P: ");
     1423
     1424      if (gcdfound == -1) //things became dense
     1425      {
     1426        if (algExtension)
     1427        {
     1428          result= GCD_Fp_extension (F, G, a);
     1429          if (extOfExt)
     1430            result= mapDown (result, primElem, imPrimElem, oldA, dest, source);
     1431          return N (d*result);
     1432        }
     1433        if (CFFactory::gettype() == GaloisFieldDomain)
     1434        {
     1435          result= GCD_GF (F, G);
     1436          if (passToGF)
     1437          {
     1438            CanonicalForm mipo= gf_mipo;
     1439            setCharacteristic (p);
     1440            Variable alpha= rootOf (mipo.mapinto());
     1441            result= GF2FalphaRep (result, alpha);
     1442          }
     1443          if (k > 1)
     1444          {
     1445            result= GFMapDown (result, k);
     1446            setCharacteristic (p, k, gf_name);
     1447          }
     1448          return N (d*result);
     1449        }
     1450        else
     1451        {
     1452          if (p >= cf_getBigPrime(0))
     1453            return N (d*sparseGCDFp (F,G));
     1454          else
     1455            return N (d*GCD_small_p (F,G));
     1456        }
     1457      }
     1458
     1459      if (gcdfound == 1)
     1460      {
     1461        TIMING_START (termination_test);
     1462        contcand= content (DD[2], Variable (1));
     1463        cand = DD[2] / contcand;
     1464        if (B_is_F)
     1465          gcdfound = fdivides( cand, G ) && cand*(DD[1]/(lcD/contcand)) == F;
     1466        else
     1467          gcdfound = fdivides( cand, F ) && cand*(DD[1]/(lcD/contcand)) == G;
     1468        TIMING_END_AND_PRINT (termination_test,
     1469                              "time for termination test EZ_P: ");
     1470
     1471        if (passToGF && gcdfound)
     1472        {
     1473          CanonicalForm mipo= gf_mipo;
     1474          setCharacteristic (p);
     1475          Variable alpha= rootOf (mipo.mapinto());
     1476          cand= GF2FalphaRep (cand, alpha);
     1477        }
     1478        if (k > 1 && gcdfound)
     1479        {
     1480          cand= GFMapDown (cand, k);
     1481          setCharacteristic (p, k, gf_name);
     1482        }
     1483        if (extOfExt && gcdfound)
     1484          cand= mapDown (cand, primElem, imPrimElem, oldA, dest, source);
     1485      }
     1486    }
     1487    delta--;
     1488    goodPointCount= 0;
     1489  }
     1490  return N (d*cand);
     1491}
     1492

  • factory/cf_gcd_smallp.cc

    rbaa2c3a rf1e3bf5  
    38563856}
    38573857
    3858 static inline
    3859 int compress4EZGCD (const CanonicalForm& F, const CanonicalForm& G, CFMap & M,
    3860                     CFMap & N, int& both_non_zero)
    3861 {
    3862   int n= tmax (F.level(), G.level());
    3863   int * degsf= new int [n + 1];
    3864   int * degsg= new int [n + 1];
    3865 
    3866   for (int i = 0; i <= n; i++)
    3867     degsf[i]= degsg[i]= 0;
    3868 
    3869   degsf= degrees (F, degsf);
    3870   degsg= degrees (G, degsg);
    3871 
    3872   both_non_zero= 0;
    3873   int f_zero= 0;
    3874   int g_zero= 0;
    3875 
    3876   for (int i= 1; i <= n; i++)
    3877   {
    3878     if (degsf[i] != 0 && degsg[i] != 0)
    3879     {
    3880       both_non_zero++;
    3881       continue;
    3882     }
    3883     if (degsf[i] == 0 && degsg[i] != 0 && i <= G.level())
    3884     {
    3885       f_zero++;
    3886       continue;
    3887     }
    3888     if (degsg[i] == 0 && degsf[i] && i <= F.level())
    3889     {
    3890       g_zero++;
    3891       continue;
    3892     }
    3893   }
    3894 
    3895   if (both_non_zero == 0)
    3896   {
    3897     delete [] degsf;
    3898     delete [] degsg;
    3899     return 0;
    3900   }
    3901 
    3902   // map Variables which do not occur in both polynomials to higher levels
    3903   int k= 1;
    3904   int l= 1;
    3905   for (int i= 1; i <= n; i++)
    3906   {
    3907     if (degsf[i] != 0 && degsg[i] == 0 && i <= F.level())
    3908     {
    3909       if (k + both_non_zero != i)
    3910       {
    3911         M.newpair (Variable (i), Variable (k + both_non_zero));
    3912         N.newpair (Variable (k + both_non_zero), Variable (i));
    3913       }
    3914       k++;
    3915     }
    3916     if (degsf[i] == 0 && degsg[i] != 0 && i <= G.level())
    3917     {
    3918       if (l + g_zero + both_non_zero != i)
    3919       {
    3920         M.newpair (Variable (i), Variable (l + g_zero + both_non_zero));
    3921         N.newpair (Variable (l + g_zero + both_non_zero), Variable (i));
    3922       }
    3923       l++;
    3924     }
    3925   }
    3926 
    3927   // sort Variables x_{i} in decreasing order of
    3928   // min(deg_{x_{i}}(f),deg_{x_{i}}(g))
    3929   int m= tmin (F.level(), G.level());
    3930   int max_min_deg;
    3931   k= both_non_zero;
    3932   l= 0;
    3933   int i= 1;
    3934   while (k > 0)
    3935   {
    3936     max_min_deg= tmin (degsf[i], degsg[i]);
    3937     while (max_min_deg == 0)
    3938     {
    3939       i++;
    3940       max_min_deg= tmin (degsf[i], degsg[i]);
    3941     }
    3942     for (int j= i + 1; j <= m; j++)
    3943     {
    3944       if ((tmin (degsf[j],degsg[j]) < max_min_deg) &&
    3945           (tmin (degsf[j], degsg[j]) != 0))
    3946       {
    3947         max_min_deg= tmin (degsf[j], degsg[j]);
    3948         l= j;
    3949       }
    3950     }
    3951 
    3952     if (l != 0)
    3953     {
    3954       if (l != k)
    3955       {
    3956         M.newpair (Variable (l), Variable(k));
    3957         N.newpair (Variable (k), Variable(l));
    3958         degsf[l]= 0;
    3959         degsg[l]= 0;
    3960         l= 0;
    3961       }
    3962       else
    3963       {
    3964         degsf[l]= 0;
    3965         degsg[l]= 0;
    3966         l= 0;
    3967       }
    3968     }
    3969     else if (l == 0)
    3970     {
    3971       if (i != k)
    3972       {
    3973         M.newpair (Variable (i), Variable (k));
    3974         N.newpair (Variable (k), Variable (i));
    3975         degsf[i]= 0;
    3976         degsg[i]= 0;
    3977       }
    3978       else
    3979       {
    3980         degsf[i]= 0;
    3981         degsg[i]= 0;
    3982       }
    3983       i++;
    3984     }
    3985     k--;
    3986   }
    3987 
    3988   delete [] degsf;
    3989   delete [] degsg;
    3990 
    3991   return both_non_zero;
    3992 }
    3993 
    3994 static inline
    3995 CanonicalForm myShift2Zero (const CanonicalForm& F, CFList& Feval,
    3996                             const CFList& evaluation)
    3997 {
    3998   CanonicalForm A= F;
    3999   int k= 2;
    4000   for (CFListIterator i= evaluation; i.hasItem(); i++, k++)
    4001     A= A (Variable (k) + i.getItem(), k);
    4002 
    4003   CanonicalForm buf= A;
    4004   Feval= CFList();
    4005   Feval.append (buf);
    4006   for (k= evaluation.length() + 1; k > 2; k--)
    4007   {
    4008     buf= mod (buf, Variable (k));
    4009     Feval.insert (buf);
    4010   }
    4011   return A;
    4012 }
    4013 
    4014 static inline
    4015 CanonicalForm myReverseShift (const CanonicalForm& F, const CFList& evaluation)
    4016 {
    4017   int l= evaluation.length() + 1;
    4018   CanonicalForm result= F;
    4019   CFListIterator j= evaluation;
    4020   for (int i= 2; i < l + 1; i++, j++)
    4021   {
    4022     if (F.level() < i)
    4023       continue;
    4024     result= result (Variable (i) - j.getItem(), i);
    4025   }
    4026   return result;
    4027 }
    4028 
    4029 static inline
    4030 Evaluation optimize4Lift (const CanonicalForm& F, CFMap & M,
    4031                     CFMap & N, const Evaluation& A)
    4032 {
    4033   int n= F.level();
    4034   int * degsf= new int [n + 1];
    4035 
    4036   for (int i = 0; i <= n; i++)
    4037     degsf[i]= 0;
    4038 
    4039   degsf= degrees (F, degsf);
    4040 
    4041   Evaluation result= Evaluation (A.min(), A.max());
    4042   ASSERT (A.min() == 2, "expected A.min() == 2");
    4043   int max_deg;
    4044   int k= n;
    4045   int l= 1;
    4046   int i= 2;
    4047   int pos= 2;
    4048   while (k > 1)
    4049   {
    4050     max_deg= degsf [i];
    4051     while (max_deg == 0)
    4052     {
    4053       i++;
    4054       max_deg= degsf [i];
    4055     }
    4056     l= i;
    4057     for (int j= i + 1; j <=  n; j++)
    4058     {
    4059       if (degsf[j] > max_deg)
    4060       {
    4061         max_deg= degsf[j];
    4062         l= j;
    4063       }
    4064     }
    4065 
    4066     if (l <= n)
    4067     {
    4068       if (l != pos)
    4069       {
    4070         result.setValue (pos, A [l]);
    4071         M.newpair (Variable (l), Variable (pos));
    4072         N.newpair (Variable (pos), Variable (l));
    4073         degsf[l]= 0;
    4074         l= 2;
    4075         if (k == 2 && n == 3)
    4076         {
    4077           result.setValue (l, A [pos]);
    4078           M.newpair (Variable (pos), Variable (l));
    4079           N.newpair (Variable (l), Variable (pos));
    4080           degsf[pos]= 0;
    4081         }
    4082       }
    4083       else
    4084       {
    4085         result.setValue (l, A [l]);
    4086         degsf [l]= 0;
    4087       }
    4088     }
    4089     pos++;
    4090     k--;
    4091     l= 2;
    4092   }
    4093 
    4094   delete [] degsf;
    4095 
    4096   return result;
    4097 }
    4098 
    4099 static inline
    4100 int Hensel_P (const CanonicalForm & UU, CFArray & G, const Evaluation & AA,
    4101               const CFArray& LeadCoeffs )
    4102 {
    4103   CFList factors;
    4104   factors.append (G[1]);
    4105   factors.append (G[2]);
    4106 
    4107   CFMap NN, MM;
    4108   Evaluation A= optimize4Lift (UU, MM, NN, AA);
    4109 
    4110   CanonicalForm U= MM (UU);
    4111   CFArray LCs= CFArray (1,2);
    4112   LCs [1]= MM (LeadCoeffs [1]);
    4113   LCs [2]= MM (LeadCoeffs [2]);
    4114 
    4115   CFList evaluation;
    4116   long termEstimate= size (U);
    4117   for (int i= A.min(); i <= A.max(); i++)
    4118   {
    4119     if (!A[i].isZero() && (getCharacteristic() > degree (U,i))) //TODO find a good estimate for getCharacteristic() <= degree (U,i)
    4120     {
    4121       termEstimate *= degree (U,i)*2;
    4122       termEstimate /= 3;
    4123     }
    4124     evaluation.append (A [i]);
    4125   }
    4126   if (termEstimate/getNumVars(U) > 500)
    4127     return -1;
    4128   CFList UEval;
    4129   CanonicalForm shiftedU= myShift2Zero (U, UEval, evaluation);
    4130 
    4131   if (size (shiftedU)/getNumVars (U) > 500)
    4132     return -1;
    4133 
    4134   CFArray shiftedLCs= CFArray (2);
    4135   CFList shiftedLCsEval1, shiftedLCsEval2;
    4136   shiftedLCs[0]= myShift2Zero (LCs[1], shiftedLCsEval1, evaluation);
    4137   shiftedLCs[1]= myShift2Zero (LCs[2], shiftedLCsEval2, evaluation);
    4138   factors.insert (1);
    4139   int liftBound= degree (UEval.getLast(), 2) + 1;
    4140   CFArray Pi;
    4141   CFMatrix M= CFMatrix (liftBound, factors.length() - 1);
    4142   CFList diophant;
    4143   CFArray lcs= CFArray (2);
    4144   lcs [0]= shiftedLCsEval1.getFirst();
    4145   lcs [1]= shiftedLCsEval2.getFirst();
    4146   nonMonicHenselLift12 (UEval.getFirst(), factors, liftBound, Pi, diophant, M,
    4147                         lcs, false);
    4148 
    4149   for (CFListIterator i= factors; i.hasItem(); i++)
    4150   {
    4151     if (!fdivides (i.getItem(), UEval.getFirst()))
    4152       return 0;
    4153   }
    4154 
    4155   int * liftBounds;
    4156   bool noOneToOne= false;
    4157   if (U.level() > 2)
    4158   {
    4159     liftBounds= new int [U.level() - 1]; /* index: 0.. U.level()-2 */
    4160     liftBounds[0]= liftBound;
    4161     for (int i= 1; i < U.level() - 1; i++)
    4162       liftBounds[i]= degree (shiftedU, Variable (i + 2)) + 1;
    4163     factors= nonMonicHenselLift2 (UEval, factors, liftBounds, U.level() - 1,
    4164                                   false, shiftedLCsEval1, shiftedLCsEval2, Pi,
    4165                                   diophant, noOneToOne);
    4166     delete [] liftBounds;
    4167     if (noOneToOne)
    4168       return 0;
    4169   }
    4170   G[1]= factors.getFirst();
    4171   G[2]= factors.getLast();
    4172   G[1]= myReverseShift (G[1], evaluation);
    4173   G[2]= myReverseShift (G[2], evaluation);
    4174   G[1]= NN (G[1]);
    4175   G[2]= NN (G[2]);
    4176   return 1;
    4177 }
    4178 
    4179 static inline
    4180 bool findeval_P (const CanonicalForm & F, const CanonicalForm & G,
    4181                  CanonicalForm & Fb, CanonicalForm & Gb, CanonicalForm & Db,
    4182                  REvaluation & b, int delta, int degF, int degG, int maxeval,
    4183                  int & count, int& k, int bound, int& l)
    4184 {
    4185   if( count == 0 && delta != 0)
    4186   {
    4187     if( count++ > maxeval )
    4188       return false;
    4189   }
    4190   if (count > 0)
    4191   {
    4192     b.nextpoint(k);
    4193     if (k == 0)
    4194       k++;
    4195     l++;
    4196     if (l > bound)
    4197     {
    4198       l= 1;
    4199       k++;
    4200       if (k > tmax (F.level(), G.level()) - 1)
    4201         return false;
    4202       b.nextpoint (k);
    4203     }
    4204     if (count++ > maxeval)
    4205       return false;
    4206   }
    4207   while( true )
    4208   {
    4209     Fb = b( F );
    4210     if( degree( Fb, 1 ) == degF )
    4211     {
    4212       Gb = b( G );
    4213       if( degree( Gb, 1 ) == degG )
    4214       {
    4215         Db = gcd( Fb, Gb );
    4216         if( delta > 0 )
    4217         {
    4218           if( degree( Db, 1 ) <= delta )
    4219             return true;
    4220         }
    4221         else
    4222           return true;
    4223       }
    4224     }
    4225     if (k == 0)
    4226       k++;
    4227     b.nextpoint(k);
    4228     l++;
    4229     if (l > bound)
    4230     {
    4231       l= 1;
    4232       k++;
    4233       if (k > tmax (F.level(), G.level()) - 1)
    4234         return false;
    4235       b.nextpoint (k);
    4236     }
    4237     if( count++ > maxeval )
    4238       return false;
    4239   }
    4240 }
    4241 
    4242 // parameters for heuristic
    4243 static int maxNumEval= 200;
    4244 static int sizePerVars1= 500; //try dense gcd if size/#variables is bigger
    4245 
    4246 /// Extended Zassenhaus GCD for finite fields
    4247 CanonicalForm EZGCD_P( const CanonicalForm & FF, const CanonicalForm & GG )
    4248 {
    4249   if (FF.isZero() && degree(GG) > 0) return GG/Lc(GG);
    4250   else if (GG.isZero() && degree (FF) > 0) return FF/Lc(FF);
    4251   else if (FF.isZero() && GG.isZero()) return FF.genOne();
    4252   if (FF.inBaseDomain() || GG.inBaseDomain()) return FF.genOne();
    4253   if (FF.isUnivariate() && fdivides(FF, GG)) return FF/Lc(FF);
    4254   if (GG.isUnivariate() && fdivides(GG, FF)) return GG/Lc(GG);
    4255   if (FF == GG) return FF/Lc(FF);
    4256 
    4257   int maxNumVars= tmax (getNumVars (FF), getNumVars (GG));
    4258   Variable a, oldA;
    4259   int sizeF= size (FF);
    4260   int sizeG= size (GG);
    4261 
    4262   if (sizeF/maxNumVars > sizePerVars1 && sizeG/maxNumVars > sizePerVars1)
    4263   {
    4264     if (hasFirstAlgVar (FF, a) || hasFirstAlgVar (GG, a))
    4265       return GCD_Fp_extension (FF, GG, a);
    4266     else if (CFFactory::gettype() == GaloisFieldDomain)
    4267       return GCD_GF (FF, GG);
    4268     else
    4269       return GCD_small_p (FF, GG);
    4270   }
    4271 
    4272   CanonicalForm F, G, f, g, d, Fb, Gb, Db, Fbt, Gbt, Dbt, B0, B, D0, lcF, lcG,
    4273                 lcD;
    4274   CFArray DD( 1, 2 ), lcDD( 1, 2 );
    4275   int degF, degG, delta, count;
    4276   int maxeval;
    4277   maxeval= tmin((getCharacteristic()/
    4278                 (int)(ilog2(getCharacteristic())*log2exp))*2, maxNumEval);
    4279   count= 0; // number of eval. used
    4280   REvaluation b, bt;
    4281   int gcdfound = 0;
    4282   Variable x = Variable(1);
    4283 
    4284   F= FF;
    4285   G= GG;
    4286 
    4287   CFMap M,N;
    4288   int smallestDegLev;
    4289   TIMING_START (ez_p_compress)
    4290   int best_level= compress4EZGCD (F, G, M, N, smallestDegLev);
    4291 
    4292   if (best_level == 0) return G.genOne();
    4293 
    4294   F= M (F);
    4295   G= M (G);
    4296   TIMING_END_AND_PRINT (ez_p_compress, "time for compression in EZ_P: ")
    4297 
    4298   TIMING_START (ez_p_content)
    4299   f = content( F, x ); g = content( G, x );
    4300   d = gcd( f, g );
    4301   F /= f; G /= g;
    4302   TIMING_END_AND_PRINT (ez_p_content, "time to extract content in EZ_P: ")
    4303 
    4304   if( F.isUnivariate() && G.isUnivariate() )
    4305   {
    4306     if( F.mvar() == G.mvar() )
    4307       d *= gcd( F, G );
    4308     else
    4309       return N (d);
    4310     return N (d);
    4311   }
    4312   if ( F.isUnivariate())
    4313   {
    4314     g= content (G,G.mvar());
    4315     return N(d*gcd(F,g));
    4316   }
    4317   if ( G.isUnivariate())
    4318   {
    4319     f= content (F,F.mvar());
    4320     return N(d*gcd(G,f));
    4321   }
    4322 
    4323   maxNumVars= tmax (getNumVars (F), getNumVars (G));
    4324   sizeF= size (F);
    4325   sizeG= size (G);
    4326 
    4327   if (sizeF/maxNumVars > sizePerVars1 && sizeG/maxNumVars > sizePerVars1)
    4328   {
    4329     if (hasFirstAlgVar (F, a) || hasFirstAlgVar (G, a))
    4330       return N (d*GCD_Fp_extension (F, G, a));
    4331     else if (CFFactory::gettype() == GaloisFieldDomain)
    4332       return N (d*GCD_GF (F, G));
    4333     else
    4334       return N (d*GCD_small_p (F, G));
    4335   }
    4336 
    4337   int dummy= 0;
    4338   if( gcd_test_one( F, G, false, dummy ) )
    4339   {
    4340     return N (d);
    4341   }
    4342 
    4343   bool passToGF= false;
    4344   bool extOfExt= false;
    4345   int p= getCharacteristic();
    4346   bool algExtension= (hasFirstAlgVar(F,a) || hasFirstAlgVar(G,a));
    4347   int k= 1;
    4348   CanonicalForm primElem, imPrimElem;
    4349   CFList source, dest;
    4350   if (p < 50 && CFFactory::gettype() != GaloisFieldDomain && !algExtension)
    4351   {
    4352     if (p == 2)
    4353       setCharacteristic (2, 12, 'Z');
    4354     else if (p == 3)
    4355       setCharacteristic (3, 4, 'Z');
    4356     else if (p == 5 || p == 7)
    4357       setCharacteristic (p, 3, 'Z');
    4358     else
    4359       setCharacteristic (p, 2, 'Z');
    4360     passToGF= true;
    4361     F= F.mapinto();
    4362     G= G.mapinto();
    4363     maxeval= 2*ipower (p, getGFDegree());
    4364   }
    4365   else if (CFFactory::gettype() == GaloisFieldDomain &&
    4366            ipower (p , getGFDegree()) < 50)
    4367   {
    4368     k= getGFDegree();
    4369     if (ipower (p, 2*k) > 50)
    4370       setCharacteristic (p, 2*k, gf_name);
    4371     else
    4372       setCharacteristic (p, 3*k, gf_name);
    4373     F= GFMapUp (F, k);
    4374     G= GFMapUp (G, k);
    4375     maxeval= tmin (2*ipower (p, getGFDegree()), maxNumEval);
    4376   }
    4377   else if (p < 50 && algExtension && CFFactory::gettype() != GaloisFieldDomain)
    4378   {
    4379     int d= degree (getMipo (a));
    4380     oldA= a;
    4381     Variable v2;
    4382     if (p == 2 && d < 6)
    4383     {
    4384       if (fac_NTL_char != p)
    4385       {
    4386         fac_NTL_char= p;
    4387         zz_p::init (p);
    4388       }
    4389       bool primFail= false;
    4390       Variable vBuf;
    4391       primElem= primitiveElement (a, vBuf, primFail);
    4392       ASSERT (!primFail, "failure in integer factorizer");
    4393       if (d < 3)
    4394       {
    4395         zz_pX NTLIrredpoly;
    4396         BuildIrred (NTLIrredpoly, d*3);
    4397         CanonicalForm newMipo= convertNTLzzpX2CF (NTLIrredpoly, Variable (1));
    4398         v2= rootOf (newMipo);
    4399       }
    4400       else
    4401       {
    4402         zz_pX NTLIrredpoly;
    4403         BuildIrred (NTLIrredpoly, d*2);
    4404         CanonicalForm newMipo= convertNTLzzpX2CF (NTLIrredpoly, Variable (1));
    4405         v2= rootOf (newMipo);
    4406       }
    4407       imPrimElem= mapPrimElem (primElem, a, v2);
    4408       extOfExt= true;
    4409     }
    4410     else if ((p == 3 && d < 4) || ((p == 5 || p == 7) && d < 3))
    4411     {
    4412       if (fac_NTL_char != p)
    4413       {
    4414         fac_NTL_char= p;
    4415         zz_p::init (p);
    4416       }
    4417       bool primFail= false;
    4418       Variable vBuf;
    4419       primElem= primitiveElement (a, vBuf, primFail);
    4420       ASSERT (!primFail, "failure in integer factorizer");
    4421       zz_pX NTLIrredpoly;
    4422       BuildIrred (NTLIrredpoly, d*2);
    4423       CanonicalForm newMipo= convertNTLzzpX2CF (NTLIrredpoly, Variable (1));
    4424       v2= rootOf (newMipo);
    4425       imPrimElem= mapPrimElem (primElem, a, v2);
    4426       extOfExt= true;
    4427     }
    4428     if (extOfExt)
    4429     {
    4430       maxeval= tmin (2*ipower (p, degree (getMipo (v2))), maxNumEval);
    4431       F= mapUp (F, a, v2, primElem, imPrimElem, source, dest);
    4432       G= mapUp (G, a, v2, primElem, imPrimElem, source, dest);
    4433       a= v2;
    4434     }
    4435   }
    4436 
    4437   lcF = LC( F, x ); lcG = LC( G, x );
    4438   lcD = gcd( lcF, lcG );
    4439 
    4440   delta = 0;
    4441   degF = degree( F, x ); degG = degree( G, x );
    4442 
    4443   if(hasFirstAlgVar(G,a))
    4444     b = REvaluation( 2, tmax(F.level(), G.level()), AlgExtRandomF( a ) );
    4445   else
    4446   { // both not in extension given by algebraic variable
    4447     if (CFFactory::gettype() != GaloisFieldDomain)
    4448       b = REvaluation( 2, tmax(F.level(), G.level()), FFRandom() );
    4449     else
    4450       b = REvaluation( 2, tmax(F.level(), G.level()), GFRandom() );
    4451   }
    4452 
    4453   CanonicalForm cand, contcand;
    4454   CanonicalForm result;
    4455   int o, t;
    4456   o= 0;
    4457   t= 1;
    4458   int goodPointCount= 0;
    4459   while( !gcdfound )
    4460   {
    4461     TIMING_START (ez_p_eval);
    4462     if( !findeval_P( F, G, Fb, Gb, Db, b, delta, degF, degG, maxeval, count, o,
    4463          maxeval/maxNumVars, t ))
    4464     { // too many eval. used --> try another method
    4465       Off (SW_USE_EZGCD_P);
    4466       result= gcd (F,G);
    4467       On (SW_USE_EZGCD_P);
    4468       if (passToGF)
    4469       {
    4470         CanonicalForm mipo= gf_mipo;
    4471         setCharacteristic (p);
    4472         Variable alpha= rootOf (mipo.mapinto());
    4473         result= GF2FalphaRep (result, alpha);
    4474       }
    4475       if (k > 1)
    4476       {
    4477         result= GFMapDown (result, k);
    4478         setCharacteristic (p, k, gf_name);
    4479       }
    4480       if (extOfExt)
    4481         result= mapDown (result, primElem, imPrimElem, oldA, dest, source);
    4482       return N (d*result);
    4483     }
    4484     TIMING_END_AND_PRINT (ez_p_eval, "time for eval point search in EZ_P1: ");
    4485     delta = degree( Db );
    4486     if (delta == degF)
    4487     {
    4488       if (degF <= degG && fdivides (F, G))
    4489       {
    4490         if (passToGF)
    4491         {
    4492           CanonicalForm mipo= gf_mipo;
    4493           setCharacteristic (p);
    4494           Variable alpha= rootOf (mipo.mapinto());
    4495           F= GF2FalphaRep (F, alpha);
    4496         }
    4497         if (k > 1)
    4498         {
    4499           F= GFMapDown (F, k);
    4500           setCharacteristic (p, k, gf_name);
    4501         }
    4502         if (extOfExt)
    4503           F= mapDown (F, primElem, imPrimElem, oldA, dest, source);
    4504         return N (d*F);
    4505       }
    4506       else
    4507         delta--;
    4508     }
    4509     else if (delta == degG)
    4510     {
    4511       if (degG <= degF && fdivides (G, F))
    4512       {
    4513         if (passToGF)
    4514         {
    4515           CanonicalForm mipo= gf_mipo;
    4516           setCharacteristic (p);
    4517           Variable alpha= rootOf (mipo.mapinto());
    4518           G= GF2FalphaRep (G, alpha);
    4519         }
    4520         if (k > 1)
    4521         {
    4522           G= GFMapDown (G, k);
    4523           setCharacteristic (p, k, gf_name);
    4524         }
    4525         if (extOfExt)
    4526           G= mapDown (G, primElem, imPrimElem, oldA, dest, source);
    4527         return N (d*G);
    4528       }
    4529       else
    4530         delta--;
    4531     }
    4532     if( delta == 0 )
    4533     {
    4534       if (passToGF)
    4535         setCharacteristic (p);
    4536       if (k > 1)
    4537         setCharacteristic (p, k, gf_name);
    4538       return N (d);
    4539     }
    4540     while( true )
    4541     {
    4542       bt = b;
    4543       TIMING_START (ez_p_eval);
    4544       if( !findeval_P(F,G,Fbt,Gbt,Dbt, bt, delta, degF, degG, maxeval, count, o,
    4545            maxeval/maxNumVars, t ))
    4546       { // too many eval. used --> try another method
    4547         Off (SW_USE_EZGCD_P);
    4548         result= gcd (F,G);
    4549         On (SW_USE_EZGCD_P);
    4550         if (passToGF)
    4551         {
    4552           CanonicalForm mipo= gf_mipo;
    4553           setCharacteristic (p);
    4554           Variable alpha= rootOf (mipo.mapinto());
    4555           result= GF2FalphaRep (result, alpha);
    4556         }
    4557         if (k > 1)
    4558         {
    4559           result= GFMapDown (result, k);
    4560           setCharacteristic (p, k, gf_name);
    4561         }
    4562         if (extOfExt)
    4563           result= mapDown (result, primElem, imPrimElem, oldA, dest, source);
    4564         return N (d*result);
    4565       }
    4566       TIMING_END_AND_PRINT (ez_p_eval, "time for eval point search in EZ_P2: ");
    4567       int dd = degree( Dbt );
    4568       if( dd == 0 )
    4569       {
    4570         if (passToGF)
    4571           setCharacteristic (p);
    4572         if (k > 1)
    4573           setCharacteristic (p, k, gf_name);
    4574         return N (d);
    4575       }
    4576       if( dd == delta )
    4577       {
    4578         goodPointCount++;
    4579         if (goodPointCount == 5)
    4580           break;
    4581       }
    4582       if( dd < delta )
    4583       {
    4584         goodPointCount= 0;
    4585         delta = dd;
    4586         b = bt;
    4587         Db = Dbt; Fb = Fbt; Gb = Gbt;
    4588       }
    4589       if (delta == degF)
    4590       {
    4591         if (degF <= degG && fdivides (F, G))
    4592         {
    4593           if (passToGF)
    4594           {
    4595             CanonicalForm mipo= gf_mipo;
    4596             setCharacteristic (p);
    4597             Variable alpha= rootOf (mipo.mapinto());
    4598             F= GF2FalphaRep (F, alpha);
    4599           }
    4600           if (k > 1)
    4601           {
    4602             F= GFMapDown (F, k);
    4603             setCharacteristic (p, k, gf_name);
    4604           }
    4605           if (extOfExt)
    4606             F= mapDown (F, primElem, imPrimElem, oldA, dest, source);
    4607           return N (d*F);
    4608         }
    4609         else
    4610           delta--;
    4611       }
    4612       else if (delta == degG)
    4613       {
    4614         if (degG <= degF && fdivides (G, F))
    4615         {
    4616           if (passToGF)
    4617           {
    4618             CanonicalForm mipo= gf_mipo;
    4619             setCharacteristic (p);
    4620             Variable alpha= rootOf (mipo.mapinto());
    4621             G= GF2FalphaRep (G, alpha);
    4622           }
    4623           if (k > 1)
    4624           {
    4625             G= GFMapDown (G, k);
    4626             setCharacteristic (p, k, gf_name);
    4627           }
    4628           if (extOfExt)
    4629             G= mapDown (G, primElem, imPrimElem, oldA, dest, source);
    4630           return N (d*G);
    4631         }
    4632         else
    4633           delta--;
    4634       }
    4635       if( delta == 0 )
    4636       {
    4637         if (passToGF)
    4638           setCharacteristic (p);
    4639         if (k > 1)
    4640           setCharacteristic (p, k, gf_name);
    4641         return N (d);
    4642       }
    4643     }
    4644     if( delta != degF && delta != degG )
    4645     {
    4646       bool B_is_F;
    4647       CanonicalForm xxx1, xxx2;
    4648       CanonicalForm buf;
    4649       DD[1] = Fb / Db;
    4650       buf= Gb/Db;
    4651       xxx1 = gcd( DD[1], Db );
    4652       xxx2 = gcd( buf, Db );
    4653       if (((xxx1.inCoeffDomain() && xxx2.inCoeffDomain()) &&
    4654           (size (F) <= size (G)))
    4655           || (xxx1.inCoeffDomain() && !xxx2.inCoeffDomain()))
    4656       {
    4657         B = F;
    4658         DD[2] = Db;
    4659         lcDD[1] = lcF;
    4660         lcDD[2] = lcD;
    4661         B_is_F = true;
    4662       }
    4663       else if (((xxx1.inCoeffDomain() && xxx2.inCoeffDomain()) &&
    4664                (size (G) < size (F)))
    4665                || (!xxx1.inCoeffDomain() && xxx2.inCoeffDomain()))
    4666       {
    4667         DD[1] = buf;
    4668         B = G;
    4669         DD[2] = Db;
    4670         lcDD[1] = lcG;
    4671         lcDD[2] = lcD;
    4672         B_is_F = false;
    4673       }
    4674       else // special case handling
    4675       {
    4676         Off (SW_USE_EZGCD_P);
    4677         result= gcd (F,G);
    4678         On (SW_USE_EZGCD_P);
    4679         if (passToGF)
    4680         {
    4681           CanonicalForm mipo= gf_mipo;
    4682           setCharacteristic (p);
    4683           Variable alpha= rootOf (mipo.mapinto());
    4684           result= GF2FalphaRep (result, alpha);
    4685         }
    4686         if (k > 1)
    4687         {
    4688           result= GFMapDown (result, k);
    4689           setCharacteristic (p, k, gf_name);
    4690         }
    4691         if (extOfExt)
    4692           result= mapDown (result, primElem, imPrimElem, oldA, dest, source);
    4693         return N (d*result);
    4694       }
    4695       DD[2] = DD[2] * ( b( lcDD[2] ) / lc( DD[2] ) );
    4696       DD[1] = DD[1] * ( b( lcDD[1] ) / lc( DD[1] ) );
    4697 
    4698       if (size (B*lcDD[2])/maxNumVars > sizePerVars1)
    4699       {
    4700         if (algExtension)
    4701         {
    4702           result= GCD_Fp_extension (F, G, a);
    4703           if (extOfExt)
    4704             result= mapDown (result, primElem, imPrimElem, oldA, dest, source);
    4705           return N (d*result);
    4706         }
    4707         if (CFFactory::gettype() == GaloisFieldDomain)
    4708         {
    4709           result= GCD_GF (F, G);
    4710           if (passToGF)
    4711           {
    4712             CanonicalForm mipo= gf_mipo;
    4713             setCharacteristic (p);
    4714             Variable alpha= rootOf (mipo.mapinto());
    4715             result= GF2FalphaRep (result, alpha);
    4716           }
    4717           if (k > 1)
    4718           {
    4719             result= GFMapDown (result, k);
    4720             setCharacteristic (p, k, gf_name);
    4721           }
    4722           return N (d*result);
    4723         }
    4724         else
    4725           return N (d*GCD_small_p (F,G));
    4726       }
    4727 
    4728       TIMING_START (ez_p_hensel_lift);
    4729       gcdfound= Hensel_P (B*lcD, DD, b, lcDD);
    4730       TIMING_END_AND_PRINT (ez_p_hensel_lift, "time for Hensel lift in EZ_P: ");
    4731 
    4732       if (gcdfound == -1) //things became dense
    4733       {
    4734         if (algExtension)
    4735         {
    4736           result= GCD_Fp_extension (F, G, a);
    4737           if (extOfExt)
    4738             result= mapDown (result, primElem, imPrimElem, oldA, dest, source);
    4739           return N (d*result);
    4740         }
    4741         if (CFFactory::gettype() == GaloisFieldDomain)
    4742         {
    4743           result= GCD_GF (F, G);
    4744           if (passToGF)
    4745           {
    4746             CanonicalForm mipo= gf_mipo;
    4747             setCharacteristic (p);
    4748             Variable alpha= rootOf (mipo.mapinto());
    4749             result= GF2FalphaRep (result, alpha);
    4750           }
    4751           if (k > 1)
    4752           {
    4753             result= GFMapDown (result, k);
    4754             setCharacteristic (p, k, gf_name);
    4755           }
    4756           return N (d*result);
    4757         }
    4758         else
    4759         {
    4760           if (p >= cf_getBigPrime(0))
    4761             return N (d*sparseGCDFp (F,G));
    4762           else
    4763             return N (d*GCD_small_p (F,G));
    4764         }
    4765       }
    4766 
    4767       if (gcdfound == 1)
    4768       {
    4769         TIMING_START (termination_test);
    4770         contcand= content (DD[2], Variable (1));
    4771         cand = DD[2] / contcand;
    4772         if (B_is_F)
    4773           gcdfound = fdivides( cand, G ) && cand*(DD[1]/(lcD/contcand)) == F;
    4774         else
    4775           gcdfound = fdivides( cand, F ) && cand*(DD[1]/(lcD/contcand)) == G;
    4776         TIMING_END_AND_PRINT (termination_test,
    4777                               "time for termination test EZ_P: ");
    4778 
    4779         if (passToGF && gcdfound)
    4780         {
    4781           CanonicalForm mipo= gf_mipo;
    4782           setCharacteristic (p);
    4783           Variable alpha= rootOf (mipo.mapinto());
    4784           cand= GF2FalphaRep (cand, alpha);
    4785         }
    4786         if (k > 1 && gcdfound)
    4787         {
    4788           cand= GFMapDown (cand, k);
    4789           setCharacteristic (p, k, gf_name);
    4790         }
    4791         if (extOfExt && gcdfound)
    4792           cand= mapDown (cand, primElem, imPrimElem, oldA, dest, source);
    4793       }
    4794     }
    4795     delta--;
    4796     goodPointCount= 0;
    4797   }
    4798   return N (d*cand);
    4799 }
    4800 
    4801 
    48023858#endif

  • factory/cf_gcd_smallp.h

    rbaa2c3a rf1e3bf5  
    100100}
    101101
    102 /// extended Zassenhaus GCD
    103 CanonicalForm
    104 EZGCD_P (const CanonicalForm& A, const CanonicalForm& B);
    105 
    106102CanonicalForm
    107103randomIrredpoly (int i, const Variable & x) ;

  • factory/fac_util.h

    rbaa2c3a rf1e3bf5  
    4747bool gcd_test_one ( const CanonicalForm & f, const CanonicalForm & g, bool swap, int & d );
    4848
    49 CanonicalForm ezgcd ( const CanonicalForm & f, const CanonicalForm & g );
    50 
    5149void extgcd ( const CanonicalForm & a, const CanonicalForm & b, CanonicalForm & S, CanonicalForm & T, const modpk & pk );
    5250
Note: See TracChangeset for help on using the changeset viewer.