Changeset 018577 in git


Ignore:
Timestamp:
Sep 13, 2010, 5:18:53 PM (13 years ago)
Author:
Martin Lee <martinlee84@…>
Branches:
(u'spielwiese', '0d6b7fcd9813a1ca1ed4220cfa2b104b97a0a003')
Children:
72292fbb6682a392b300d1837489860c3f412907
Parents:
231b0fd0173552d615fba69018a64ff5510211be
Message:
new sparse modular gcd over finite fields, new EZ gcd over finite fields, which is now switched on by default


git-svn-id: file:///usr/local/Singular/svn/trunk@13185 2c84dea3-7e68-4137-9b89-c4e89433aadc
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • Singular/tesths.cc

    r231b0fd r018577  
    134134  On(SW_USE_EZGCD);
    135135  On(SW_USE_CHINREM_GCD);
     136  //On(SW_USE_FF_MOD_GCD);
    136137  //On(SW_USE_fieldGCD);
    137   //On(SW_USE_EZGCD_P);
     138  On(SW_USE_EZGCD_P);
    138139  On(SW_USE_QGCD);
    139140  Off(SW_USE_NTL_SORT); // may be changed by an command line option

  • factory/cf_gcd.cc

    r231b0fd r018577  
    2929bool isPurePoly(const CanonicalForm & );
    3030static CanonicalForm gcd_univar_ntl0( const CanonicalForm &, const CanonicalForm & );
    31 static CanonicalForm gcd_univar_ntlp( const CanonicalForm &, const CanonicalForm & );
    3231#endif
    3332
    … …  
    367366
    368367static CanonicalForm
    369 gcd_poly_p( const CanonicalForm & f, const CanonicalForm & g )
    370 {
    371     CanonicalForm pi, pi1;
    372     CanonicalForm C, Ci, Ci1, Hi, bi, pi2;
    373     bool bpure;
    374     int delta = degree( f ) - degree( g );
    375 
    376     if ( delta >= 0 )
    377     {
    378         pi = f; pi1 = g;
    379     }
    380     else
    381     {
    382         pi = g; pi1 = f; delta = -delta;
    383     }
    384     Ci = content( pi ); Ci1 = content( pi1 );
    385     pi1 = pi1 / Ci1; pi = pi / Ci;
    386     C = gcd( Ci, Ci1 );
    387     if ( !( pi.isUnivariate() && pi1.isUnivariate() ) )
    388     {
    389         //out_cf("F:",f,"\n");
    390         //out_cf("G:",g,"\n");
    391         //out_cf("newGCD:",newGCD(f,g),"\n");
    392         if (isOn(SW_USE_GCD_P) && (getCharacteristic()>0))
    393         {
    394           return newGCD(f,g);
    395         }
    396         if ( gcd_test_one( pi1, pi, true ) )
    397         {
    398           C=abs(C);
    399           //out_cf("GCD:",C,"\n");
    400           return C;
    401         }
    402         bpure = false;
    403     }
    404     else
    405     {
    406         bpure = isPurePoly(pi) && isPurePoly(pi1);
    407 #ifdef HAVE_NTL
    408         if ( isOn(SW_USE_NTL_GCD_P) && bpure && (CFFactory::gettype() != GaloisFieldDomain))
    409             return gcd_univar_ntlp(pi, pi1 ) * C;
    410 #endif
    411     }
    412     Variable v = f.mvar();
    413     Hi = power( LC( pi1, v ), delta );
    414     if ( (delta+1) % 2 )
    415         bi = 1;
    416     else
    417         bi = -1;
    418     while ( degree( pi1, v ) > 0 )
    419     {
    420         pi2 = psr( pi, pi1, v );
    421         pi2 = pi2 / bi;
    422         pi = pi1; pi1 = pi2;
    423         if ( degree( pi1, v ) > 0 )
    424         {
    425             delta = degree( pi, v ) - degree( pi1, v );
    426             if ( (delta+1) % 2 )
    427                 bi = LC( pi, v ) * power( Hi, delta );
    428             else
    429                 bi = -LC( pi, v ) * power( Hi, delta );
    430             Hi = power( LC( pi1, v ), delta ) / power( Hi, delta-1 );
    431         }
    432     }
    433     if ( degree( pi1, v ) == 0 )
    434     {
    435       C=abs(C);
    436       //out_cf("GCD:",C,"\n");
    437       return C;
    438     }
    439     pi /= content( pi );
    440     if ( bpure )
    441         pi /= pi.lc();
    442     C=abs(C*pi);
    443     //out_cf("GCD:",C,"\n");
    444     return C;
    445 }
    446 
    447 static CanonicalForm
    448368gcd_poly_0( const CanonicalForm & f, const CanonicalForm & g )
    449369{
    … …  
    566486    else if (isOn( SW_USE_EZGCD_P ) && (!fc_and_gc_Univariate))
    567487    {
    568       if ( pe == 1 )
    569         fc = fin_ezgcd( fc, gc );
    570       else if ( pe > 0 )// no variable at position 1
    571       {
    572         fc = replacevar( fc, Variable(pe), Variable(1) );
    573         gc = replacevar( gc, Variable(pe), Variable(1) );
    574         fc = replacevar( fin_ezgcd( fc, gc ), Variable(1), Variable(pe) );
    575       }
    576       else
    577       {
    578         pe = -pe;
    579         fc = swapvar( fc, Variable(pe), Variable(1) );
    580         gc = swapvar( gc, Variable(pe), Variable(1) );
    581         fc = swapvar( fin_ezgcd( fc, gc ), Variable(1), Variable(pe) );
    582       }
     488      fc= EZGCD_P (fc, gc);
    583489    }
    584490    else if (isOn(SW_USE_GCD_P))
    … …  
    926832}
    927833
    928 static CanonicalForm
    929 gcd_univar_ntlp( const CanonicalForm & F, const CanonicalForm & G )
    930 {
    931   if (fac_NTL_char!=getCharacteristic())
    932   {
    933     fac_NTL_char=getCharacteristic();
    934     #ifdef NTL_ZZ
    935     ZZ r;
    936     r=getCharacteristic();
    937     ZZ_pContext ccc(r);
    938     #else
    939     zz_pContext ccc(getCharacteristic());
    940     #endif
    941     ccc.restore();
    942     #ifdef NTL_ZZ
    943     ZZ_p::init(r);
    944     #else
    945     zz_p::init(getCharacteristic());
    946     #endif
    947   }
    948   #ifdef NTL_ZZ
    949   ZZ_pX F1=convertFacCF2NTLZZpX(F);
    950   ZZ_pX G1=convertFacCF2NTLZZpX(G);
    951   ZZ_pX R=GCD(F1,G1);
    952   return  convertNTLZZpX2CF(R,F.mvar());
    953   #else
    954   zz_pX F1=convertFacCF2NTLzzpX(F);
    955   zz_pX G1=convertFacCF2NTLzzpX(G);
    956   zz_pX R=GCD(F1,G1);
    957   return  convertNTLzzpX2CF(R,F.mvar());
    958   #endif
    959 }
    960 
    961834#endif
    962835

  • factory/cf_gcd_smallp.cc

    r231b0fd r018577  
    3030#include "ftmpl_functions.h"
    3131#include "cf_random.h"
    32 
    33 // iinline helper functions:
     32#include "cf_algorithm.h"
     33#include "ffreval.h"
     34#include "fac_iterfor.h"
     35#include "cf_binom.h"
     36#include "facHensel.h"
     37#include "facFqFactorizeUtil.h"
     38
     39// inline helper functions:
    3440#include "cf_map_ext.h"
    3541
    … …  
    4652TIMING_DEFINE_PRINT(newton_interpolation);
    4753
     54CanonicalForm
     55gcd_univar_ntlp( const CanonicalForm & F, const CanonicalForm & G )
     56{
     57  if (fac_NTL_char!=getCharacteristic())
     58  {
     59    fac_NTL_char=getCharacteristic();
     60    #ifdef NTL_ZZ
     61    ZZ r;
     62    r=getCharacteristic();
     63    ZZ_pContext ccc(r);
     64    #else
     65    zz_pContext ccc(getCharacteristic());
     66    #endif
     67    ccc.restore();
     68    #ifdef NTL_ZZ
     69    ZZ_p::init(r);
     70    #else
     71    zz_p::init(getCharacteristic());
     72    #endif
     73  }
     74  #ifdef NTL_ZZ
     75  ZZ_pX F1=convertFacCF2NTLZZpX(F);
     76  ZZ_pX G1=convertFacCF2NTLZZpX(G);
     77  ZZ_pX R=GCD(F1,G1);
     78  return  convertNTLZZpX2CF(R,F.mvar());
     79  #else
     80  zz_pX F1=convertFacCF2NTLzzpX(F);
     81  zz_pX G1=convertFacCF2NTLzzpX(G);
     82  zz_pX R=GCD(F1,G1);
     83  return  convertNTLzzpX2CF(R,F.mvar());
     84  #endif
     85}
     86
     87#endif
     88
     89CanonicalForm
     90gcd_poly_p( const CanonicalForm & f, const CanonicalForm & g)
     91{
     92    CanonicalForm pi, pi1;
     93    CanonicalForm C, Ci, Ci1, Hi, bi, pi2;
     94    bool bpure;
     95    int delta = degree( f ) - degree( g );
     96
     97    if ( delta >= 0 )
     98    {
     99        pi = f; pi1 = g;
     100    }
     101    else
     102    {
     103        pi = g; pi1 = f; delta = -delta;
     104    }
     105    Ci = content( pi ); Ci1 = content( pi1 );
     106    pi1 = pi1 / Ci1; pi = pi / Ci;
     107    C = gcd( Ci, Ci1 );
     108    if ( !( pi.isUnivariate() && pi1.isUnivariate() ) )
     109    {
     110        //out_cf("F:",f,"\n");
     111        //out_cf("G:",g,"\n");
     112        //out_cf("newGCD:",newGCD(f,g),"\n");
     113        if ( gcd_test_one( pi1, pi, true ) )
     114        {
     115          C=abs(C);
     116          //out_cf("GCD:",C,"\n");
     117          return C;
     118        }
     119        bpure = false;
     120    }
     121    else
     122    {
     123        bpure = isPurePoly(pi) && isPurePoly(pi1);
     124#ifdef HAVE_NTL
     125        if ( isOn(SW_USE_NTL_GCD_P) && bpure && (CFFactory::gettype() != GaloisFieldDomain))
     126            return gcd_univar_ntlp(pi, pi1 ) * C;
     127#endif
     128    }
     129    Variable v = f.mvar();
     130    Hi = power( LC( pi1, v ), delta );
     131    if ( (delta+1) % 2 )
     132        bi = 1;
     133    else
     134        bi = -1;
     135    while ( degree( pi1, v ) > 0 )
     136    {
     137        pi2 = psr( pi, pi1, v );
     138        pi2 = pi2 / bi;
     139        pi = pi1; pi1 = pi2;
     140        if ( degree( pi1, v ) > 0 )
     141        {
     142            delta = degree( pi, v ) - degree( pi1, v );
     143            if ( (delta+1) % 2 )
     144                bi = LC( pi, v ) * power( Hi, delta );
     145            else
     146                bi = -LC( pi, v ) * power( Hi, delta );
     147            Hi = power( LC( pi1, v ), delta ) / power( Hi, delta-1 );
     148        }
     149    }
     150    if ( degree( pi1, v ) == 0 )
     151    {
     152      C=abs(C);
     153      //out_cf("GCD:",C,"\n");
     154      return C;
     155    }
     156    pi /= content( pi );
     157    if ( bpure )
     158        pi /= pi.lc();
     159    C=abs(C*pi);
     160    //out_cf("GCD:",C,"\n");
     161    return C;
     162}
     163
     164#ifdef HAVE_NTL
    48165/// compressing two polynomials F and G, M is used for compressing,
    49166/// N to reverse the compression
    50167static inline
    51168int myCompress (const CanonicalForm& F, const CanonicalForm& G, CFMap & M,
    52                 CFMap & N, bool& top_level)
     169                CFMap & N, const bool& topLevel)
    53170{
    54171  int n= tmax (F.level(), G.level());
    … …  
    67184  int both_zero= 0;
    68185
    69   if (top_level)
     186  if (topLevel)
    70187  {
    71188    for (int i= 1; i <= n; i++)
    … …  
    198315  delete [] degsg;
    199316
    200   return 1;
     317  if (topLevel)
     318    return both_non_zero;
     319  else
     320    return 1;
    201321}
     322
     323CanonicalForm
     324uniCoeff (const CanonicalForm& F)
     325{
     326  if (F.inCoeffDomain())
     327    return F;
     328  if (F.isUnivariate() && F.level() == 1)
     329    return F;
     330  else if (F.isUnivariate () && F.level() != 1)
     331    return 0;
     332  CanonicalForm result= 0;
     333  if (F.level() == 2)
     334  {
     335    for (CFIterator i= F; i.hasTerms(); i++)
     336      result += i.coeff();
     337  }
     338  else
     339  {
     340    for (CFIterator i= F; i.hasTerms(); i++)
     341      result += uniCoeff (i.coeff());
     342  }
     343  return result;
     344}
     345
     346int
     347substituteCheck (const CanonicalForm& F, const CanonicalForm& G)
     348{
     349  if (F.inBaseDomain() || G.inBaseDomain())
     350    return 0;
     351  Variable x= Variable (1);
     352  CanonicalForm f= swapvar (F, F.mvar(), x); //TODO swapping seems to be pretty expensive
     353  CanonicalForm g= swapvar (G, G.mvar(), x);
     354  int sizef= 0;
     355  int sizeg= 0;
     356  for (CFIterator i= f; i.hasTerms(); i++, sizef++)
     357  {
     358    if (i.exp() == 1)
     359      return 0;
     360  }
     361  for (CFIterator i= g; i.hasTerms(); i++, sizeg++)
     362  {
     363    if (i.exp() == 1)
     364      return 0;
     365  }
     366  int * expf= new int [sizef];
     367  int * expg= new int [sizeg];
     368  int j= 0;
     369  for (CFIterator i= f; i.hasTerms(); i++, j++)
     370  {
     371    expf [j]= i.exp();
     372  }
     373  j= 0;
     374  for (CFIterator i= g; i.hasTerms(); i++, j++)
     375  {
     376    expg [j]= i.exp();
     377  }
     378 
     379  int indf= sizef - 1;
     380  int indg= sizeg - 1;
     381  if (expf[indf] == 0)
     382    indf--;
     383  if (expg[indg] == 0)
     384    indg--;
     385   
     386  if (expg[indg] != expf [indf] || (expg[indg] == 1 && expf[indf] == 1))
     387    return 0;
     388 
     389  // expf[indg] == expf[indf] after here
     390  int result= expg[indg];
     391  for (int i= indf - 1; i >= 0; i--)
     392  {
     393    if (expf [i]%result != 0)
     394      return 0;
     395  }
     396 
     397  for (int i= indg - 1; i >= 0; i--)
     398  {
     399    if (expg [i]%result != 0)
     400      return 0;
     401  }
     402
     403  delete [] expg;
     404  delete [] expf;
     405  return result;
     406}
     407
     408// substiution
     409void
     410subst (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& A,
     411       CanonicalForm& B, const int d
     412      )
     413{
     414  if (d == 1)
     415  {
     416    A= F;
     417    B= G;
     418    return;
     419  }
     420 
     421  CanonicalForm C= 0;
     422  CanonicalForm D= 0; 
     423  Variable x= Variable (1);
     424  CanonicalForm f= swapvar (F, x, F.mvar());
     425  CanonicalForm g= swapvar (G, x, G.mvar());
     426  for (CFIterator i= f; i.hasTerms(); i++)
     427    C += i.coeff()*power (f.mvar(), i.exp()/ d);
     428  for (CFIterator i= g; i.hasTerms(); i++)
     429    D += i.coeff()*power (g.mvar(), i.exp()/ d);
     430  A= swapvar (C, x, F.mvar());
     431  B= swapvar (D, x, G.mvar());
     432}
     433
     434CanonicalForm
     435reverseSubst (const CanonicalForm& F, const int d)
     436{
     437  if (d == 1)
     438    return F;
     439 
     440  Variable x= Variable (1);
     441  CanonicalForm f= swapvar (F, x, F.mvar());
     442  CanonicalForm result= 0;
     443  for (CFIterator i= f; i.hasTerms(); i++)
     444    result += i.coeff()*power (f.mvar(), d*i.exp());
     445  return swapvar (result, x, F.mvar());
     446 
     447}
     448
     449static inline CanonicalForm
     450uni_content (const CanonicalForm & F);
     451
     452CanonicalForm
     453uni_content (const CanonicalForm& F, const Variable& x)
     454{
     455  if (F.inBaseDomain())
     456    return F.genOne();
     457  if (F.level() == x.level() && F.isUnivariate())
     458    return F;
     459  if (F.level() != x.level() && F.isUnivariate())
     460    return F.genOne();
     461   
     462  if (x.level() != 1)
     463  {
     464    CanonicalForm f= swapvar (F, x, Variable (1));
     465    CanonicalForm result= uni_content (f);
     466    return swapvar (result, x, Variable (1));
     467  }
     468  else
     469    return uni_content (F);
     470}
    202471
    203472/// compute the content of F, where F is considered as an element of
    … …  
    212481  if (F.level() != 1 && F.isUnivariate())
    213482    return F.genOne();
    214   if (degree (F,1) == 0) return F.genOne();
     483  if (degree (F, 1) == 0) return F.genOne();
    215484
    216485  int l= F.level();
    … …  
    233502}
    234503
     504CanonicalForm
     505extractContents (const CanonicalForm& F, const CanonicalForm& G,
     506                 CanonicalForm& contentF, CanonicalForm& contentG,
     507                 CanonicalForm& ppF, CanonicalForm& ppG, const int d)
     508{
     509  CanonicalForm uniContentF, uniContentG, gcdcFcG;
     510  contentF= 1;
     511  contentG= 1;
     512  ppF= F;
     513  ppG= G;
     514  CanonicalForm result= 1;
     515  for (int i= 1; i <= d; i++)
     516  {
     517    uniContentF= uni_content (F, Variable (i));
     518    uniContentG= uni_content (G, Variable (i));
     519    gcdcFcG= gcd (uniContentF, uniContentG);
     520    contentF *= uniContentF;
     521    contentG *= uniContentG;
     522    ppF /= uniContentF;
     523    ppG /= uniContentG;
     524    result *= gcdcFcG;
     525  }
     526  return result;
     527}
     528
    235529/// compute the leading coefficient of F, where F is considered as an element
    236530/// of \f$ R[x_{1}][x_{2},\ldots ,x_{n}] \f$, order on Mon (x_{2},\ldots ,x_{n})
    237 /// is dp.
     531/// is Dp.
    238532static inline
    239533CanonicalForm uni_lcoeff (const CanonicalForm& F)
    … …  
    286580  int p= getCharacteristic ();
    287581  int d= degree (mipo);
    288   int bound= ipower (p, d);
     582  double bound= pow ((double) p, d);
    289583  do
    290584  {
    … …  
    326620  ZZ_pX NTLirredpoly;
    327621  //TODO: replace d by max_{i} (deg_x{i}(f))
    328   int i= (int) (log ((double) ipower (d + 1, num_vars))/log ((double) p));
     622  int i= (int) (log ((double) pow ((double) d + 1, (double) num_vars))/log ((double) p));
    329623  int m= degree (getMipo (beta));
    330624  if (i <= 1)
    … …  
    336630
    337631/// GCD of F and G over \f$ F_{p}(\alpha ) \f$ ,
    338 /// l and top_level are only used internally, output is monic
     632/// l and topLevel are only used internally, output is monic
    339633/// based on Alg. 7.2. as described in "Algorithms for
    340634/// Computer Algebra" by Geddes, Czapor, Labahn
    341635CanonicalForm
    342636GCD_Fp_extension (const CanonicalForm& F, const CanonicalForm& G,
    343                   Variable & alpha, CFList& l, bool& top_level)
     637                  Variable & alpha, CFList& l, bool& topLevel)
    344638{
    345639  CanonicalForm A= F;
    … …  
    354648 
    355649  CFMap M,N;
    356   int best_level= myCompress (A, B, M, N, top_level);
     650  int best_level= myCompress (A, B, M, N, topLevel);
    357651
    358652  if (best_level == 0) return B.genOne();
    … …  
    366660  if (A.isUnivariate() && B.isUnivariate())
    367661    return N (gcd(A,B));
    368  
     662
     663  int substitute= substituteCheck (A, B);
     664 
     665  if (substitute > 1)
     666    subst (A, B, A, B, substitute);
     667   
    369668  CanonicalForm cA, cB;    // content of A and B
    370669  CanonicalForm ppA, ppB;    // primitive part of A and B
    371670  CanonicalForm gcdcAcB;
    372 
    373   cA = uni_content (A);
    374   cB = uni_content (B);
    375 
    376   gcdcAcB= gcd (cA, cB);
    377 
    378   ppA= A/cA;
    379   ppB= B/cB;
     671  if (topLevel)
     672  {
     673    if (best_level <= 2)
     674      gcdcAcB= extractContents (A, B, cA, cB, ppA, ppB, best_level);
     675    else
     676      gcdcAcB= extractContents (A, B, cA, cB, ppA, ppB, 2); // TODO needs some more tuning
     677  }
     678  else
     679  {
     680    cA = uni_content (A);
     681    cB = uni_content (B);
     682    gcdcAcB= gcd (cA, cB);
     683    ppA= A/cA;
     684    ppB= B/cB;
     685  }
    380686
    381687  CanonicalForm lcA, lcB;  // leading coefficients of A and B
    … …  
    399705  int d= totaldegree (ppA, Variable(2), Variable (ppA.level()));
    400706
    401   if (d == 0) return N(gcdcAcB);
     707  if (d == 0)
     708  {
     709    if (substitute > 1)
     710      return N (reverseSubst (gcdcAcB, substitute));
     711    else
     712      return N(gcdcAcB);
     713  }
    402714  int d0= totaldegree (ppB, Variable(2), Variable (ppB.level()));
    403715  if (d0 < d)
    404716    d= d0;
    405   if (d == 0) return N(gcdcAcB);
     717  if (d == 0)
     718  {
     719    if (substitute > 1)
     720      return N (reverseSubst (gcdcAcB, substitute));
     721    else
     722      return N(gcdcAcB);
     723  }
    406724
    407725  CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH;
    … …  
    413731  H= 0;
    414732  bool fail= false;
    415   top_level= false;
     733  topLevel= false;
    416734  bool inextension= false;
    417735  Variable V_buf= alpha;
    … …  
    427745      int num_vars= tmin (getNumVars(A), getNumVars(B));
    428746      num_vars--;
    429 
     747     
    430748      choose_extension (d, num_vars, V_buf);
    431749      bool prim_fail= false;
    432750      Variable V_buf2;
    433751      prim_elem= primitiveElement (alpha, V_buf2, prim_fail);
    434 
    435752      ASSERT (!prim_fail, "failure in integer factorizer");
    436753      if (prim_fail)
    … …  
    438755      else
    439756        im_prim_elem= mapPrimElem (prim_elem, alpha, V_buf);
    440 
     757     
    441758      DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (alpha));
    442759      DEBOUTLN (cerr, "getMipo (V_buf2)= " << getMipo (V_buf2));
    … …  
    461778      G_random_element=
    462779      GCD_Fp_extension (ppA (random_element, x), ppB (random_element, x), V_buf,
    463                         list, top_level);
     780                        list, topLevel);
    464781      TIMING_END_AND_PRINT (gcd_recursion,
    465782                            "time for recursive call: ");
    … …  
    472789      G_random_element=
    473790      GCD_Fp_extension (ppA(random_element, x), ppB(random_element, x), V_buf,
    474                         list, top_level);
     791                        list, topLevel);
    475792      TIMING_END_AND_PRINT (gcd_recursion,
    476793                            "time for recursive call: ");
    … …  
    480797    d0= totaldegree (G_random_element, Variable(2),
    481798                     Variable (G_random_element.level()));
    482     if (d0 == 0) return N(gcdcAcB);
     799    if (d0 == 0)
     800    {
     801      if (substitute > 1)
     802        return N (reverseSubst (gcdcAcB, substitute));
     803      else
     804        return N(gcdcAcB);
     805    }
    483806    if (d0 >  d)
    484807    {
    … …  
    521844        ppH /= Lc(ppH);
    522845        DEBOUTLN (cerr, "ppH after mapDown= " << ppH);
    523         if (fdivides (ppH, A) && fdivides (ppH, B))         
     846        if (fdivides (ppH, A) && fdivides (ppH, B))
     847        {
     848          if (substitute > 1)
     849          {
     850            ppH= reverseSubst (ppH, substitute);
     851            gcdcAcB= reverseSubst (gcdcAcB, substitute);
     852          }
    524853          return N(gcdcAcB*ppH);
    525       }
    526       else if (fdivides (ppH, A) && fdivides (ppH, B))
    527         return N(gcdcAcB*(ppH/Lc(ppH)));
     854        }
     855      }
     856      else if (fdivides (ppH, A) && fdivides (ppH, B))
     857      {
     858        if (substitute > 1)
     859        {
     860          ppH= reverseSubst (ppH, substitute);
     861          gcdcAcB= reverseSubst (gcdcAcB, substitute);
     862        }
     863        return N(gcdcAcB*ppH);
     864      }
    528865    }
    529866
    … …  
    549886  int p= getCharacteristic();
    550887  int d= getGFDegree();
    551   int bound= ipower (p, d);
     888  double bound= pow ((double) p, (double) d);
    552889  do
    553890  {
    … …  
    580917/// the size of the base field is bounded by 2^16, output is monic
    581918CanonicalForm GCD_GF (const CanonicalForm& F, const CanonicalForm& G,
    582         CFList& l, bool& top_level)
     919        CFList& l, bool& topLevel)
    583920{
    584921  CanonicalForm A= F;
    … …  
    593930 
    594931  CFMap M,N;
    595   int best_level= myCompress (A, B, M, N, top_level);
     932  int best_level= myCompress (A, B, M, N, topLevel);
    596933
    597934  if (best_level == 0) return B.genOne();
    … …  
    652989  H= 0;
    653990  bool fail= false;
    654   top_level= false;
     991  topLevel= false;
    655992  bool inextension= false;
    656993  int p;
    … …  
    6671004      num_vars--;
    6681005      p= getCharacteristic();
    669       expon= (int) floor ((log ((double) ipower (d + 1, num_vars))/log ((double) p)));
     1006      expon= (int) floor ((log ((double) pow ((double) d + 1, num_vars))/log ((double) p)));
    6701007      if (expon < 2)
    6711008        expon= 2;
    6721009      kk= getGFDegree();
    673       if (ipower (p, kk*expon) < (1 << 16))
     1010      if (pow ( (double) p, kk*expon) < (1 << 16))
    6741011        setCharacteristic (p, kk*(int)expon, 'b');
    6751012      else
    … …  
    6941031      TIMING_START (gcd_recursion);
    6951032      G_random_element= GCD_GF (ppA(random_element, x), ppB(random_element, x),
    696                                 list, top_level);
     1033                                list, topLevel);
    6971034      TIMING_END_AND_PRINT (gcd_recursion,
    6981035                            "time for recursive call: ");
    … …  
    7041041      TIMING_START (gcd_recursion);
    7051042      G_random_element= GCD_GF (ppA(random_element, x), ppB(random_element, x),
    706                                 list, top_level);
     1043                                list, topLevel);
    7071044      TIMING_END_AND_PRINT (gcd_recursion,
    7081045                            "time for recursive call: ");
    … …  
    8291166
    8301167CanonicalForm GCD_small_p (const CanonicalForm& F, const CanonicalForm&  G,
    831                            bool& top_level, CFList& l)
     1168                           bool& topLevel, CFList& l)
    8321169{
    8331170  CanonicalForm A= F;
    … …  
    8421179
    8431180  CFMap M,N;
    844   int best_level= myCompress (A, B, M, N, top_level);
     1181  int best_level= myCompress (A, B, M, N, topLevel);
    8451182
    8461183  if (best_level == 0) return B.genOne();
    … …  
    8481185  A= M(A);
    8491186  B= M(B);
    850 
    851   Variable x= Variable (1);
    8521187
    8531188  //univariate case
    8541189  if (A.isUnivariate() && B.isUnivariate())
    8551190    return N (gcd (A, B));
     1191   
     1192  int substitute= substituteCheck (A, B);
     1193 
     1194  if (substitute > 1)
     1195    subst (A, B, A, B, substitute);
     1196   
     1197  Variable x= Variable (1);
    8561198
    8571199  CanonicalForm cA, cB;    // content of A and B
    8581200  CanonicalForm ppA, ppB;    // primitive part of A and B
    8591201  CanonicalForm gcdcAcB;
    860   cA = uni_content (A);
    861   cB = uni_content (B);
    862   gcdcAcB= gcd (cA, cB);
    863   ppA= A/cA;
    864   ppB= B/cB;
     1202  if (topLevel)
     1203  {
     1204    if (best_level <= 2)
     1205      gcdcAcB= extractContents (A, B, cA, cB, ppA, ppB, best_level);
     1206    else
     1207      gcdcAcB= extractContents (A, B, cA, cB, ppA, ppB, 2);
     1208  }
     1209  else
     1210  {
     1211    cA = uni_content (A);
     1212    cB = uni_content (B);
     1213    gcdcAcB= gcd (cA, cB);
     1214    ppA= A/cA;
     1215    ppB= B/cB;
     1216  }
    8651217
    8661218  CanonicalForm lcA, lcB;  // leading coefficients of A and B
    … …  
    8791231      return G/Lc(G);
    8801232  }
    881 
     1233 
    8821234  gcdlcAlcB= gcd (lcA, lcB);
    8831235 
    … …  
    8851237  int d0;
    8861238
    887   if (d == 0) return N(gcdcAcB);
     1239  if (d == 0)
     1240  {
     1241    if (substitute > 1)
     1242      return N (reverseSubst (gcdcAcB, substitute));
     1243    else
     1244      return N(gcdcAcB);
     1245  }
    8881246  d0= totaldegree (ppB, Variable (2), Variable (ppB.level()));
    8891247
    … …  
    8911249    d= d0;
    8921250
    893   if (d == 0) return N(gcdcAcB);
     1251  if (d == 0)
     1252  {
     1253    if (substitute > 1)
     1254      return N (reverseSubst (gcdcAcB, substitute));
     1255    else
     1256      return N(gcdcAcB);
     1257  }
    8941258
    8951259  CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH;
    … …  
    8991263  G_m= 0;
    9001264  Variable alpha, V_buf;
    901   int p, k;
    9021265  bool fail= false;
    9031266  bool inextension= false;
    9041267  bool inextensionextension= false;
    905   top_level= false;
     1268  topLevel= false;
    9061269  CFList source, dest;
    9071270  do
    … …  
    9171280      TIMING_START (gcd_recursion);
    9181281      G_random_element=
    919       GCD_small_p (ppA (random_element,x), ppB (random_element,x), top_level,
     1282      GCD_small_p (ppA (random_element,x), ppB (random_element,x), topLevel,
    9201283      list);
    9211284      TIMING_END_AND_PRINT (gcd_recursion,
    … …  
    9291292      G_random_element=
    9301293      GCD_Fp_extension (ppA (random_element, x), ppB (random_element, x), alpha,
    931                         list, top_level);
     1294                        list, topLevel);
    9321295      TIMING_END_AND_PRINT (gcd_recursion,
    9331296                            "time for recursive call: ");
    … …  
    9421305      int deg= 2;
    9431306      do {
     1307        l= CFList();
    9441308        mipo= randomIrredpoly (deg, x);
    9451309        alpha= rootOf (mipo);
    … …  
    9481312        random_element= randomElement (m, alpha, l, fail);
    9491313        deg++;
    950       } while (fail); 
     1314      } while (fail);
    9511315      list= CFList();
    9521316      TIMING_START (gcd_recursion);
    9531317      G_random_element=
    9541318      GCD_Fp_extension (ppA (random_element, x), ppB (random_element, x), alpha,
    955                         list, top_level);
     1319                        list, topLevel);
    9561320      TIMING_END_AND_PRINT (gcd_recursion,
    9571321                            "time for recursive call: ");
    … …  
    9701334      CanonicalForm prim_elem, im_prim_elem;
    9711335      prim_elem= primitiveElement (alpha, V_buf2, prim_fail);
    972        
    9731336      ASSERT (!prim_fail, "failure in integer factorizer");
    9741337      if (prim_fail)
    … …  
    9991362      G_random_element=
    10001363      GCD_Fp_extension (ppA (random_element, x), ppB (random_element, x), V_buf,
    1001                         list, top_level);
     1364                        list, topLevel);
    10021365      TIMING_END_AND_PRINT (gcd_recursion,
    10031366                            "time for recursive call: ");
    … …  
    10111374    if (d0 == 0)
    10121375    {
    1013       return N(gcdcAcB);
     1376      if (substitute > 1)
     1377        return N (reverseSubst (gcdcAcB, substitute));
     1378      else
     1379        return N(gcdcAcB);
    10141380    }
    10151381    if (d0 >  d)
    … …  
    10471413      ppH /= Lc (ppH);
    10481414      DEBOUTLN (cerr, "ppH= " << ppH);
    1049       if (fdivides (ppH, A) && fdivides (ppH, B))
     1415      if (fdivides (ppH, A) && fdivides (ppH, B))
     1416      {
     1417        if (substitute > 1)
     1418        {
     1419          ppH= reverseSubst (ppH, substitute);
     1420          gcdcAcB= reverseSubst (gcdcAcB, substitute);
     1421        }
    10501422        return N(gcdcAcB*ppH);
     1423      }
    10511424    }
    10521425
    … …  
    10591432}
    10601433
     1434CFArray
     1435solveVandermonde (const CFArray& M, const CFArray& A)
     1436{
     1437  int r= M.size();
     1438  ASSERT (A.size() == r, "vector does not have right size");
     1439
     1440  if (r == 1)
     1441  {
     1442    CFArray result= CFArray (1);
     1443    result [0]= A [0] / M [0];
     1444    return result;
     1445  }
     1446  // check solvability
     1447  bool notDistinct= false;
     1448  for (int i= 0; i < r - 1; i++)
     1449  {
     1450    for (int j= i + 1; j < r; j++)
     1451    {
     1452      if (M [i] == M [j])
     1453      {
     1454        notDistinct= true;
     1455        break;
     1456      }
     1457    }
     1458  }
     1459  if (notDistinct)
     1460    return CFArray();
     1461
     1462  CanonicalForm master= 1;
     1463  Variable x= Variable (1);
     1464  for (int i= 0; i < r; i++)
     1465    master *= x - M [i];
     1466  CFList Pj;
     1467  CanonicalForm tmp;
     1468  for (int i= 0; i < r; i++)
     1469  {
     1470    tmp= master/(x - M [i]);
     1471    tmp /= tmp (M [i], 1);
     1472    Pj.append (tmp);
     1473  } 
     1474  CFArray result= CFArray (r);
     1475 
     1476  CFListIterator j= Pj;
     1477  for (int i= 1; i <= r; i++, j++)
     1478  {
     1479    tmp= 0;
     1480    for (int l= 0; l < A.size(); l++)
     1481      tmp += A[l]*j.getItem()[l];
     1482    result[i - 1]= tmp;
     1483  }
     1484  return result;
     1485}
     1486
     1487/// solve M*x=L, where M is a Vandermonde matrix
     1488CFArray
     1489solveVandermonde (const CFMatrix& M, const CFArray& A)
     1490{
     1491  int r= M.rows();
     1492  int c= M.columns();
     1493  ASSERT (c == r, "number of columns and rows not equal");
     1494  ASSERT (A.size() == r, "vector does not have right size");
     1495
     1496  if (c == 1)
     1497  {
     1498    CFArray result= CFArray (1);
     1499    result [0]= A [0] / M (1,1);
     1500    return result;
     1501  }
     1502  // check solvability
     1503  CFMatrix secondRow= M (2, 2, 1, c);
     1504  bool notDistinct= false;
     1505  for (int i= 1; i < r; i++)
     1506  {
     1507    for (int j= i + 1; j <= r; j++)
     1508    {
     1509      if (secondRow (1, i) == secondRow (1, j))
     1510      {
     1511        notDistinct= true;
     1512        break;
     1513      }
     1514    }
     1515  }
     1516  if (notDistinct)
     1517    return CFArray();
     1518
     1519  CanonicalForm master= 1;
     1520  Variable x= Variable (1);
     1521  for (int i= 1; i <= r; i++)
     1522    master *= x - secondRow (1, i);
     1523  CFList Pj;
     1524  CanonicalForm tmp;
     1525  for (int i= 1; i <= r; i++)
     1526  {
     1527    tmp= master/(x - secondRow (1, i));
     1528    tmp /= tmp (secondRow (1, i), 1);
     1529    Pj.append (tmp);
     1530  } 
     1531  CFArray result= CFArray (r);
     1532 
     1533  CFListIterator j= Pj;
     1534  for (int i= 1; i <= r; i++, j++)
     1535  {
     1536    tmp= 0;
     1537    for (int l= 0; l < A.size(); l++)
     1538      tmp += A[l]*j.getItem()[l];
     1539    result[i - 1]= tmp;
     1540  }
     1541  return result;
     1542}
     1543
     1544CFArray
     1545solveGeneralVandermonde (const CFArray& M, const CFArray& A)
     1546{
     1547  int r= M.size();
     1548  ASSERT (c == r, "number of columns and rows not equal");
     1549  ASSERT (A.size() == r, "vector does not have right size");
     1550  if (r == 1)
     1551  {
     1552    CFArray result= CFArray (1);
     1553    result [0]= A[0] / M [0];
     1554    return result;
     1555  }
     1556  // check solvability
     1557  bool notDistinct= false;
     1558  for (int i= 0; i < r - 1; i++)
     1559  {
     1560    for (int j= i + 1; j < r; j++)
     1561    {
     1562      if (M [i] == M [j])
     1563      {
     1564        notDistinct= true;
     1565        break;
     1566      }
     1567    }
     1568  }
     1569  if (notDistinct)
     1570    return CFArray();
     1571
     1572  CanonicalForm master= 1;
     1573  Variable x= Variable (1);
     1574  for (int i= 0; i < r; i++)
     1575    master *= x - M [i];
     1576  master *= x;
     1577  CFList Pj;
     1578  CanonicalForm tmp;
     1579  for (int i= 0; i < r; i++)
     1580  {
     1581    tmp= master/(x - M [i]);
     1582    tmp /= tmp (M [i], 1);
     1583    Pj.append (tmp);
     1584  } 
     1585 
     1586  CFArray result= CFArray (r);
     1587
     1588  CFListIterator j= Pj;
     1589  for (int i= 1; i <= r; i++, j++)
     1590  {
     1591    tmp= 0;
     1592
     1593    for (int l= 1; l <= A.size(); l++)
     1594      tmp += A[l - 1]*j.getItem()[l];
     1595    result[i - 1]= tmp;
     1596  }
     1597  return result;
     1598}
     1599
     1600CFArray
     1601solveGeneralVandermonde (const CFMatrix& M, const CFArray& A)
     1602{
     1603  int r= M.rows();
     1604  int c= M.columns();
     1605  ASSERT (c == r, "number of columns and rows not equal");
     1606  ASSERT (A.size() == r, "vector does not have right size");
     1607  if (c == 1)
     1608  {
     1609    CFArray result= CFArray (1);
     1610    result [0]= A[0] / M (1,1);
     1611    return result;
     1612  }
     1613  // check solvability
     1614  CFMatrix firstRow= M (1, 1, 1, c);
     1615  bool notDistinct= false;
     1616  for (int i= 1; i < r; i++)
     1617  {
     1618    for (int j= i + 1; j <= r; j++)
     1619    {
     1620      if (firstRow (1, i) == firstRow (1, j))
     1621      {
     1622        notDistinct= true;
     1623        break;
     1624      }
     1625    }
     1626  }
     1627  if (notDistinct)
     1628    return CFArray();
     1629
     1630  CanonicalForm master= 1;
     1631  Variable x= Variable (1);
     1632  for (int i= 1; i <= r; i++)
     1633    master *= x - firstRow (1, i);
     1634  master *= x;
     1635  CFList Pj;
     1636  CanonicalForm tmp;
     1637  for (int i= 1; i <= r; i++)
     1638  {
     1639    tmp= master/(x - firstRow (1, i));
     1640    tmp /= tmp (firstRow (1, i), 1);
     1641    Pj.append (tmp);
     1642  } 
     1643 
     1644  CFArray result= CFArray (r);
     1645  CFListIterator j= Pj;
     1646  for (int i= 1; i <= r; i++, j++)
     1647  {
     1648    tmp= 0;
     1649    for (int l= 1; l <= A.size(); l++)
     1650      tmp += A[l - 1]*j.getItem()[l];
     1651    result[i - 1]= tmp;
     1652  }
     1653  return result;
     1654}
     1655
     1656/// M in row echolon form, rk rank of M
     1657CFArray
     1658readOffSolution (const CFMatrix& M, const long rk)
     1659{
     1660  CFArray result= CFArray (rk);
     1661  CanonicalForm tmp1, tmp2, tmp3;
     1662  for (int i= rk; i >= 1; i--)
     1663  {
     1664    tmp3= 0;
     1665    tmp1= M (i, M.columns());
     1666    for (int j= M.columns() - 1; j >= 1; j--)
     1667    {
     1668      tmp2= M (i, j);
     1669      if (j == i)
     1670        break;
     1671      else
     1672        tmp3 += tmp2*result[j - 1];
     1673    }
     1674    result[i - 1]= (tmp1 - tmp3)/tmp2;
     1675  }
     1676  return result;
     1677}
     1678
     1679CFArray
     1680readOffSolution (const CFMatrix& M, const CFArray& L, const CFArray& partialSol)
     1681{
     1682  CFArray result= CFArray (M.rows());
     1683  CanonicalForm tmp1, tmp2, tmp3;
     1684  int k;
     1685  for (int i= M.rows(); i >= 1; i--)
     1686  {
     1687    tmp3= 0;
     1688    tmp1= L[i - 1];
     1689    k= 0;
     1690    for (int j= M.columns(); j >= 1; j--, k++)
     1691    {
     1692      tmp2= M (i, j);
     1693      if (j == i)
     1694        break;
     1695      else
     1696      {
     1697        if (k > partialSol.size() - 1)
     1698          tmp3 += tmp2*result[j - 1];
     1699        else
     1700          tmp3 += tmp2*partialSol[partialSol.size() - k - 1];
     1701      }
     1702    }
     1703    result [i - 1]= (tmp1 - tmp3)/tmp2;
     1704  }
     1705  return result;
     1706}
     1707
     1708long
     1709gaussianElimFp (CFMatrix& M, CFArray& L)
     1710{
     1711  ASSERT (L.size() <= M.rows(), "dimension exceeded");
     1712  CFMatrix *N;
     1713  N= new CFMatrix (M.rows(), M.columns() + 1);
     1714
     1715  for (int i= 1; i <= M.rows(); i++)
     1716    for (int j= 1; j <= M.columns(); j++)
     1717      (*N) (i, j)= M (i, j);
     1718 
     1719  int j= 1;
     1720  for (int i= 0; i < L.size(); i++, j++)
     1721    (*N) (j, M.columns() + 1)= L[i];
     1722  int p= getCharacteristic ();
     1723  zz_p::init (p);
     1724  mat_zz_p *NTLN= convertFacCFMatrix2NTLmat_zz_p(*N);
     1725  long rk= gauss (*NTLN);
     1726 
     1727  N= convertNTLmat_zz_p2FacCFMatrix (*NTLN);
     1728 
     1729  L= CFArray (M.rows());
     1730  for (int i= 0; i < M.rows(); i++)
     1731    L[i]= (*N) (i + 1, M.columns() + 1);
     1732  M= (*N) (1, M.rows(), 1, M.columns());
     1733  return rk;
     1734}
     1735
     1736long
     1737gaussianElimFq (CFMatrix& M, CFArray& L, const Variable& alpha)
     1738{
     1739  ASSERT (L.size() <= M.rows(), "dimension exceeded");
     1740  CFMatrix *N;
     1741  N= new CFMatrix (M.rows(), M.columns() + 1);
     1742
     1743  for (int i= 1; i <= M.rows(); i++)
     1744    for (int j= 1; j <= M.columns(); j++)
     1745      (*N) (i, j)= M (i, j);
     1746   
     1747  int j= 1;
     1748  for (int i= 0; i < L.size(); i++, j++)
     1749    (*N) (j, M.columns() + 1)= L[i];
     1750  int p= getCharacteristic ();
     1751  zz_p::init (p);
     1752  zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
     1753  zz_pE::init (NTLMipo);
     1754  mat_zz_pE *NTLN= convertFacCFMatrix2NTLmat_zz_pE(*N);
     1755  long rk= gauss (*NTLN);
     1756
     1757  N= convertNTLmat_zz_pE2FacCFMatrix (*NTLN, alpha);
     1758
     1759  M= (*N) (1, M.rows(), 1, M.columns());
     1760  L= CFArray (M.rows());
     1761  for (int i= 0; i < M.rows(); i++)
     1762    L[i]= (*N) (i + 1, M.columns() + 1);
     1763  return rk;
     1764}
     1765
     1766CFArray
     1767solveSystemFp (const CFMatrix& M, const CFArray& L)
     1768{
     1769  ASSERT (L.size() <= M.rows(), "dimension exceeded");
     1770  CFMatrix *N;
     1771  N= new CFMatrix (M.rows(), M.columns() + 1);
     1772
     1773  for (int i= 1; i <= M.rows(); i++)
     1774    for (int j= 1; j <= M.columns(); j++)
     1775      (*N) (i, j)= M (i, j);
     1776
     1777  int j= 1;
     1778  for (int i= 0; i < L.size(); i++, j++)
     1779    (*N) (j, M.columns() + 1)= L[i];
     1780  int p= getCharacteristic ();
     1781  zz_p::init (p);
     1782  mat_zz_p *NTLN= convertFacCFMatrix2NTLmat_zz_p(*N);
     1783  long rk= gauss (*NTLN);   
     1784  if (rk != M.columns()) 
     1785    return CFArray();
     1786   
     1787  N= convertNTLmat_zz_p2FacCFMatrix (*NTLN);
     1788
     1789  CFArray A= readOffSolution (*N, rk);
     1790
     1791  return A;
     1792}
     1793
     1794CFArray
     1795solveSystemFq (const CFMatrix& M, const CFArray& L, const Variable& alpha)
     1796{
     1797  ASSERT (L.size() <= M.rows(), "dimension exceeded");
     1798  CFMatrix *N;
     1799  N= new CFMatrix (M.rows(), M.columns() + 1);
     1800
     1801  for (int i= 1; i <= M.rows(); i++)
     1802    for (int j= 1; j <= M.columns(); j++)
     1803      (*N) (i, j)= M (i, j);
     1804  int j= 1;
     1805  for (int i= 0; i < L.size(); i++, j++)
     1806    (*N) (j, M.columns() + 1)= L[i];
     1807  int p= getCharacteristic ();
     1808  zz_p::init (p);
     1809  zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
     1810  zz_pE::init (NTLMipo);
     1811  mat_zz_pE *NTLN= convertFacCFMatrix2NTLmat_zz_pE(*N);
     1812  long rk= gauss (*NTLN);
     1813  if (rk != M.columns()) 
     1814    return CFArray();
     1815 
     1816  N= convertNTLmat_zz_pE2FacCFMatrix (*NTLN, alpha);
     1817 
     1818  CFArray A= readOffSolution (*N, rk);
     1819 
     1820  return A;
     1821}
     1822
     1823CFArray
     1824getMonoms (const CanonicalForm& F)
     1825{
     1826  if (F.inCoeffDomain())
     1827  {
     1828    CFArray result= CFArray (1);
     1829    result [0]= 1;
     1830    return result;
     1831  }
     1832  if (F.isUnivariate())
     1833  {
     1834    CFArray result= CFArray (size(F));
     1835    int j= 0;
     1836    for (CFIterator i= F; i.hasTerms(); i++, j++)
     1837      result[j]= power (F.mvar(), i.exp());
     1838    return result;
     1839  }
     1840  int numMon= size (F);
     1841  CFArray result= CFArray (numMon);
     1842  int j= 0;
     1843  CFArray recResult;
     1844  Variable x= F.mvar();
     1845  CanonicalForm powX;
     1846  for (CFIterator i= F; i.hasTerms(); i++)
     1847  {
     1848    powX= power (x, i.exp());
     1849    recResult= getMonoms (i.coeff());
     1850    for (int k= 0; k < recResult.size(); k++)
     1851      result[j+k]= powX*recResult[k];
     1852    j += recResult.size();
     1853  }
     1854  return result;
     1855}
     1856
     1857CFArray
     1858evaluateMonom (const CanonicalForm& F, const CFList& evalPoints)
     1859{
     1860  if (F.inCoeffDomain())
     1861  {
     1862    CFArray result= CFArray (1);
     1863    result [0]= F;
     1864    return result;
     1865  }
     1866  if (F.isUnivariate())
     1867  {
     1868    ASSERT (evalPoints.length() == 1,
     1869            "expected an eval point with only one component");
     1870    CFArray result= CFArray (size(F));
     1871    int j= 0;
     1872    CanonicalForm evalPoint= evalPoints.getLast();
     1873    for (CFIterator i= F; i.hasTerms(); i++, j++)
     1874      result[j]= power (evalPoint, i.exp());
     1875    return result;
     1876  }
     1877  int numMon= size (F);
     1878  CFArray result= CFArray (numMon);
     1879  int j= 0;
     1880  CanonicalForm evalPoint= evalPoints.getLast();
     1881  CFList buf= evalPoints;
     1882  buf.removeLast();
     1883  CFArray recResult;
     1884  CanonicalForm powEvalPoint;
     1885  for (CFIterator i= F; i.hasTerms(); i++)
     1886  {
     1887    powEvalPoint= power (evalPoint, i.exp());
     1888    recResult= evaluateMonom (i.coeff(), buf);
     1889    for (int k= 0; k < recResult.size(); k++)
     1890      result[j+k]= powEvalPoint*recResult[k];
     1891    j += recResult.size();
     1892  }
     1893  return result;
     1894}
     1895
     1896CFArray
     1897evaluate (const CFArray& A, const CFList& evalPoints)
     1898{
     1899  CFArray result= A.size();
     1900  CanonicalForm tmp;
     1901  int k;
     1902  for (int i= 0; i < A.size(); i++)
     1903  {
     1904    tmp= A[i];
     1905    k= 1;
     1906    for (CFListIterator j= evalPoints; j.hasItem(); j++, k++)
     1907      tmp= tmp (j.getItem(), k);
     1908    result[i]= tmp;
     1909  }
     1910  return result;
     1911}
     1912
     1913CFList
     1914evaluationPoints (const CanonicalForm& F, const CanonicalForm& G,
     1915                  CanonicalForm& Feval, CanonicalForm& Geval,
     1916                  const CanonicalForm& LCF, const bool& GF,
     1917                  const Variable& alpha, bool& fail, CFList& list
     1918                 )
     1919{
     1920  int k= tmax (F.level(), G.level()) - 1;
     1921  Variable x= Variable (1);
     1922  CFList result;
     1923  FFRandom genFF;
     1924  GFRandom genGF;
     1925  int p= getCharacteristic ();
     1926  double bound;
     1927  if (alpha != Variable (1))
     1928  {
     1929    bound= pow ((double) p, (double) degree (getMipo(alpha)));
     1930    bound= pow ((double) bound, (double) k);
     1931  }
     1932  else if (GF)
     1933  {
     1934    bound= pow ((double) p, (double) getGFDegree());
     1935    bound= pow ((double) bound, (double) k);
     1936  }
     1937  else
     1938    bound= pow ((double) p, (double) k);
     1939
     1940  CanonicalForm random;
     1941  int j;
     1942  bool zeroOneOccured= false;
     1943  bool allEqual= false;
     1944  CanonicalForm buf;
     1945  do
     1946  {
     1947    random= 0;
     1948    // possible overflow if list.length() does not fit into a int
     1949    if (list.length() >= bound)
     1950    {
     1951      fail= true;
     1952      break;
     1953    }
     1954    for (int i= 0; i < k; i++)
     1955    {
     1956      if (GF)
     1957      {
     1958        result.append (genGF.generate());
     1959        random += result.getLast()*power (x, i);
     1960      }
     1961      else if (alpha != Variable(1))
     1962      {
     1963        AlgExtRandomF genAlgExt (alpha);
     1964        result.append (genAlgExt.generate());
     1965        random += result.getLast()*power (x, i);
     1966      }
     1967      else
     1968      {
     1969        result.append (genFF.generate());
     1970        random += result.getLast()*power (x, i);
     1971      }
     1972      if (result.getLast().isOne() || result.getLast().isZero())
     1973        zeroOneOccured= true;
     1974    }
     1975    if (find (list, random))
     1976    {
     1977      zeroOneOccured= false;
     1978      allEqual= false;
     1979      result= CFList();
     1980      continue;
     1981    }
     1982    if (zeroOneOccured)
     1983    {
     1984      list.append (random);
     1985      zeroOneOccured= false;
     1986      allEqual= false;
     1987      result= CFList();
     1988      continue;
     1989    }
     1990    // no zero at this point
     1991    if (k > 1)
     1992    {
     1993      allEqual= true;
     1994      CFIterator iter= random;
     1995      buf= iter.coeff();
     1996      iter++;
     1997      for (; iter.hasTerms(); iter++)
     1998        if (buf != iter.coeff())
     1999          allEqual= false;
     2000    }
     2001    if (allEqual)
     2002    {
     2003      list.append (random);
     2004      allEqual= false;
     2005      zeroOneOccured= false;
     2006      result= CFList();
     2007      continue;
     2008    }
     2009   
     2010    Feval= F;
     2011    Geval= G;
     2012    CanonicalForm LCeval= LCF;
     2013    j= 1;
     2014    for (CFListIterator i= result; i.hasItem(); i++, j++)
     2015    {
     2016      Feval= Feval (i.getItem(), j);
     2017      Geval= Geval (i.getItem(), j);
     2018      LCeval= LCeval (i.getItem(), j);
     2019    }
     2020
     2021    if (LCeval.isZero())
     2022    {
     2023      if (!find (list, random))
     2024        list.append (random);
     2025      zeroOneOccured= false;
     2026      allEqual= false;
     2027      result= CFList();
     2028      continue;
     2029    }
     2030
     2031    if (list.length() >= bound)
     2032    {
     2033      fail= true;
     2034      break;
     2035    }
     2036  } while (find (list, random));
     2037
     2038  return result;
     2039}
     2040
     2041/*CFList
     2042evalPointsSmallChar (const CanonicalForm& F, const CanonicalForm& G,
     2043                     CanonicalForm& Feval, CanonicalForm& Geval,
     2044                     const CanonicalForm& LCF, const bool& GF,
     2045                     const Variable& alpha, bool& fail, CFList& list
     2046                    )
     2047{
     2048  int k= tmax (F.level(), G.level()) - 1;
     2049  Variable x= Variable (1);
     2050  CFList result;
     2051  FFRandom genFF;
     2052  GFRandom genGF;
     2053  int p= getCharacteristic ();
     2054  double bound;
     2055  if (alpha != Variable (1))
     2056  {
     2057    bound= pow ((double) p, (double) degree (getMipo(alpha)));
     2058    bound= pow ((double) bound, (double) k);
     2059  }
     2060  else if (GF)
     2061  {
     2062    bound= pow ((double) p, (double) getGFDegree());
     2063    bound= pow ((double) bound, (double) k);
     2064  }
     2065  else
     2066    bound= pow ((double) p, (double) k);
     2067
     2068  CanonicalForm random;
     2069  int j;
     2070  do
     2071  {
     2072    random= 0;
     2073    // possible overflow if list.length() does not fit into a int
     2074    if (list.length() >= bound)
     2075    {
     2076      fail= true;
     2077      break;
     2078    }
     2079    for (int i= 0; i < k; i++)
     2080    {
     2081      if (GF)
     2082      {
     2083        result.append (genGF.generate());
     2084        random += result.getLast()*power (x, i);
     2085      }
     2086      else if (alpha != Variable(1))
     2087      {
     2088        AlgExtRandomF genAlgExt (alpha);
     2089        result.append (genAlgExt.generate());
     2090        random += result.getLast()*power (x, i);
     2091      }
     2092      else
     2093      {
     2094        result.append (genFF.generate());
     2095        random += result.getLast()*power (x, i);
     2096      }
     2097    }
     2098    if (find (list, random))
     2099    {
     2100      result= CFList();
     2101      continue;
     2102    }
     2103   
     2104    Feval= F;
     2105    Geval= G;
     2106    CanonicalForm LCeval= LCF;
     2107    j= 1;
     2108    for (CFListIterator i= result; i.hasItem(); i++, j++)
     2109    {
     2110      Feval= Feval (i.getItem(), j);
     2111      Geval= Geval (i.getItem(), j);
     2112      LCeval= LCeval (i.getItem(), j);
     2113    }
     2114
     2115    if (LCeval.isZero())
     2116    {
     2117      if (!find (list, random))
     2118        list.append (random);
     2119      result= CFList();
     2120      continue;
     2121    }
     2122
     2123    if (list.length() >= bound)
     2124    {
     2125      fail= true;
     2126      break;
     2127    }
     2128  } while (find (list, random));
     2129
     2130  return result;
     2131}*/
     2132
     2133
     2134/// multiply two lists componentwise
     2135void mult (CFList& L1, const CFList& L2)
     2136{
     2137  ASSERT (L1.length() == L2.length(), "lists of the same size expected");
     2138 
     2139  CFListIterator j= L2;
     2140  for (CFListIterator i= L1; i.hasItem(); i++, j++)
     2141    i.getItem() *= j.getItem();
     2142}
     2143
     2144void eval (const CanonicalForm& A, const CanonicalForm& B, CanonicalForm& Aeval,
     2145           CanonicalForm& Beval, const CFList& L)
     2146{
     2147  Aeval= A;
     2148  Beval= B;
     2149  int j= 1;
     2150  for (CFListIterator i= L; i.hasItem(); i++, j++)
     2151  {
     2152    Aeval= Aeval (i.getItem(), j);
     2153    Beval= Beval (i.getItem(), j);
     2154  }
     2155}
     2156
     2157CanonicalForm
     2158monicSparseInterpol (const CanonicalForm& F, const CanonicalForm& G,
     2159                     const CanonicalForm& skeleton, const Variable& alpha,
     2160                     bool& fail, CFArray*& coeffMonoms, CFArray& Monoms
     2161                    )
     2162{
     2163  CanonicalForm A= F;
     2164  CanonicalForm B= G;
     2165  if (F.isZero() && degree(G) > 0) return G/Lc(G);
     2166  else if (G.isZero() && degree (F) > 0) return F/Lc(F);
     2167  else if (F.isZero() && G.isZero()) return F.genOne();
     2168  if (F.inBaseDomain() || G.inBaseDomain()) return F.genOne();
     2169  if (F.isUnivariate() && fdivides(F, G)) return F/Lc(F);
     2170  if (G.isUnivariate() && fdivides(G, F)) return G/Lc(G);
     2171  if (F == G) return F/Lc(F);
     2172 
     2173  CFMap M,N;
     2174  int best_level= myCompress (A, B, M, N, false);
     2175 
     2176  if (best_level == 0)
     2177    return B.genOne();
     2178
     2179  A= M(A);
     2180  B= M(B);
     2181
     2182  Variable x= Variable (1);
     2183  ASSERT (degree (A, y) == 0, "expected degree (F, 1) == 0");
     2184  ASSERT (degree (B, y) == 0, "expected degree (G, 1) == 0");
     2185 
     2186  //univariate case
     2187  if (A.isUnivariate() && B.isUnivariate())
     2188    return N (gcd (A, B));
     2189   
     2190  CanonicalForm skel= M(skeleton);
     2191  CanonicalForm cA, cB;    // content of A and B
     2192  CanonicalForm ppA, ppB;    // primitive part of A and B
     2193  CanonicalForm gcdcAcB;
     2194  cA = uni_content (A);
     2195  cB = uni_content (B);
     2196  gcdcAcB= gcd (cA, cB);
     2197  ppA= A/cA;
     2198  ppB= B/cB;
     2199
     2200  CanonicalForm lcA, lcB;  // leading coefficients of A and B
     2201  CanonicalForm gcdlcAlcB;
     2202  lcA= uni_lcoeff (ppA);
     2203  lcB= uni_lcoeff (ppB);
     2204
     2205  if (fdivides (lcA, lcB))
     2206  {
     2207    if (fdivides (A, B))
     2208      return F/Lc(F);
     2209  }   
     2210  if (fdivides (lcB, lcA))
     2211  {
     2212    if (fdivides (B, A))
     2213      return G/Lc(G);
     2214  }
     2215
     2216  gcdlcAlcB= gcd (lcA, lcB);
     2217  int skelSize= size (skel, skel.mvar());
     2218
     2219  int j= 0;
     2220  int biggestSize= 0;
     2221 
     2222  for (CFIterator i= skel; i.hasTerms(); i++, j++)
     2223    biggestSize= tmax (biggestSize, size (i.coeff()));
     2224
     2225  CanonicalForm g, Aeval, Beval;
     2226
     2227  CFList evalPoints;
     2228  bool evalFail= false;
     2229  CFList list;
     2230  bool GF= false;
     2231  CanonicalForm LCA= LC (A);
     2232  CanonicalForm tmp;
     2233  CFArray gcds= CFArray (biggestSize);
     2234  CFList * pEvalPoints= new CFList [biggestSize];
     2235  Variable V_buf= alpha;
     2236  CFList source, dest;
     2237  CanonicalForm prim_elem, im_prim_elem;
     2238  for (int i= 0; i < biggestSize; i++)
     2239  {
     2240    if (i == 0)
     2241      evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf, evalFail,
     2242                                    list);
     2243    else
     2244    {
     2245      mult (evalPoints, pEvalPoints [0]);
     2246      eval (A, B, Aeval, Beval, evalPoints);
     2247    }
     2248   
     2249    if (evalFail)
     2250    {
     2251      if (V_buf != Variable (1))
     2252      {
     2253        do
     2254        {
     2255          int num_vars= tmin (getNumVars(A), getNumVars(B));
     2256          int d= totaldegree (A, Variable(2), Variable (A.level()));
     2257          d= tmin (d, totaldegree (B, Variable(2), Variable (B.level())));
     2258          Variable V_buf2= V_buf;
     2259          choose_extension (d, num_vars, V_buf2);
     2260          source= CFList();
     2261          dest= CFList();
     2262       
     2263          bool prim_fail= false;
     2264          Variable V_buf3;
     2265          prim_elem= primitiveElement (V_buf, V_buf3, prim_fail);
     2266       
     2267          ASSERT (!prim_fail, "failure in integer factorizer");
     2268          if (prim_fail)
     2269            ; //ERROR
     2270          else
     2271            im_prim_elem= mapPrimElem (prim_elem, V_buf, V_buf2);
     2272
     2273          DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf));
     2274          DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf2));
     2275
     2276          for (CFListIterator j= list; j.hasItem(); j++)
     2277            j.getItem()= mapUp (j.getItem(), V_buf, V_buf2, prim_elem,
     2278                                im_prim_elem, source, dest);
     2279          for (int k= 0; k < i; k++)
     2280          {
     2281            for (CFListIterator j= pEvalPoints[k]; j.hasItem(); j++)
     2282              j.getItem()= mapUp (j.getItem(), V_buf, V_buf2, prim_elem,
     2283                                  im_prim_elem, source, dest);
     2284            gcds[k]= mapUp (gcds[k], V_buf, V_buf2, prim_elem, im_prim_elem,
     2285                            source, dest); 
     2286          }                   
     2287                                   
     2288          if (alpha != Variable (1))
     2289          {                   
     2290            A= mapUp (A, V_buf, V_buf2, prim_elem, im_prim_elem, source,dest);
     2291            B= mapUp (B, V_buf, V_buf2, prim_elem, im_prim_elem, source,dest);
     2292          }                   
     2293          evalFail= false;
     2294          evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf,
     2295                                        evalFail, list);
     2296        } while (evalFail);
     2297      }
     2298      else
     2299      {
     2300        CanonicalForm mipo;
     2301        int deg= 2;
     2302        do {
     2303          mipo= randomIrredpoly (deg, x);
     2304          V_buf= rootOf (mipo);
     2305          evalFail= false;
     2306          evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf,
     2307                                        evalFail, list);
     2308          deg++;
     2309        } while (evalFail);
     2310      }
     2311    }
     2312   
     2313    g= gcd (Aeval, Beval);
     2314    g /= Lc (g);
     2315
     2316    if (degree (g) != degree (skel) || (skelSize != size (g)))
     2317    {
     2318      delete[] pEvalPoints;
     2319      fail= true;
     2320      return 0;
     2321    }
     2322    CFIterator l= skel;
     2323    for (CFIterator k= g; k.hasTerms(); k++, l++)
     2324    {
     2325      if (k.exp() != l.exp())
     2326      {
     2327        delete[] pEvalPoints;
     2328        fail= true;
     2329        return 0;
     2330      }
     2331    }
     2332    pEvalPoints[i]= evalPoints;
     2333    gcds[i]= g;
     2334
     2335    tmp= 0;
     2336    int j= 0;
     2337    for (CFListIterator k= evalPoints; k.hasItem(); k++, j++)
     2338      tmp += k.getItem()*power (x, j);
     2339    list.append (tmp);
     2340  }
     2341
     2342  if (Monoms.size() == 0)
     2343    Monoms= getMonoms (skel);
     2344  if (coeffMonoms == NULL)
     2345    coeffMonoms= new CFArray [skelSize];
     2346  j= 0;
     2347  for (CFIterator i= skel; i.hasTerms(); i++, j++)
     2348    coeffMonoms[j]= getMonoms (i.coeff());
     2349
     2350  CFArray* pL= new CFArray [skelSize];
     2351  CFArray* pM= new CFArray [skelSize];
     2352  for (int i= 0; i < biggestSize; i++)
     2353  {
     2354    CFIterator l= gcds [i];
     2355    evalPoints= pEvalPoints [i];
     2356    for (int k= 0; k < skelSize; k++, l++)
     2357    {
     2358      if (i == 0)
     2359        pL[k]= CFArray (biggestSize);
     2360      pL[k] [i]= l.coeff();
     2361
     2362      if (i == 0)
     2363        pM[k]= evaluate (coeffMonoms [k], evalPoints);
     2364    }
     2365  }
     2366 
     2367  CFArray solution;
     2368  CanonicalForm result= 0;
     2369  int ind= 0;
     2370  CFArray bufArray;
     2371  CFMatrix Mat;
     2372  for (int k= 0; k < skelSize; k++)
     2373  {
     2374    if (biggestSize != coeffMonoms[k].size())
     2375    {
     2376      bufArray= CFArray (coeffMonoms[k].size());
     2377      for (int i= 0; i < coeffMonoms[k].size(); i++)
     2378        bufArray [i]= pL[k] [i];
     2379      solution= solveGeneralVandermonde (pM [k], bufArray);
     2380    }
     2381    else
     2382      solution= solveGeneralVandermonde (pM [k], pL[k]);
     2383
     2384    if (solution.size() == 0)
     2385    {
     2386      delete[] pEvalPoints;
     2387      delete[] pM;
     2388      delete[] pL;
     2389      delete[] coeffMonoms;
     2390      fail= true;
     2391      return 0;
     2392    }
     2393    for (int l= 0; l < solution.size(); l++)
     2394      result += solution[l]*Monoms [ind + l];
     2395    ind += solution.size();
     2396  }
     2397
     2398  delete[] pEvalPoints;
     2399  delete[] pM;
     2400  delete[] pL;
     2401 
     2402  if (alpha != Variable (1) && V_buf != alpha)
     2403  {
     2404    CFList u, v;
     2405    result= mapDown (result, prim_elem, im_prim_elem, alpha, u, v);
     2406  }
     2407 
     2408  result= N(result);
     2409  if (fdivides (result, F) && fdivides (result, G))
     2410    return result;
     2411  else
     2412  {
     2413    delete[] coeffMonoms;
     2414    fail= true;
     2415    return 0;
     2416  }
     2417}
     2418
     2419CanonicalForm
     2420nonMonicSparseInterpol (const CanonicalForm& F, const CanonicalForm& G,
     2421                        const CanonicalForm& skeleton, const Variable& alpha,
     2422                        bool& fail, CFArray*& coeffMonoms, CFArray& Monoms
     2423                       )
     2424{
     2425  CanonicalForm A= F;
     2426  CanonicalForm B= G;
     2427  if (F.isZero() && degree(G) > 0) return G/Lc(G);
     2428  else if (G.isZero() && degree (F) > 0) return F/Lc(F);
     2429  else if (F.isZero() && G.isZero()) return F.genOne();
     2430  if (F.inBaseDomain() || G.inBaseDomain()) return F.genOne();
     2431  if (F.isUnivariate() && fdivides(F, G)) return F/Lc(F);
     2432  if (G.isUnivariate() && fdivides(G, F)) return G/Lc(G);
     2433  if (F == G) return F/Lc(F);
     2434 
     2435  CFMap M,N;
     2436  int best_level= myCompress (A, B, M, N, false);
     2437 
     2438  if (best_level == 0)
     2439    return B.genOne();
     2440
     2441  A= M(A);
     2442  B= M(B);
     2443 
     2444  Variable x= Variable (1);
     2445  ASSERT (degree (A, y) == 0, "expected degree (F, 1) == 0");
     2446  ASSERT (degree (B, y) == 0, "expected degree (G, 1) == 0");
     2447 
     2448  //univariate case
     2449  if (A.isUnivariate() && B.isUnivariate())
     2450    return N (gcd (A, B));
     2451
     2452  CanonicalForm skel= M(skeleton);
     2453
     2454  CanonicalForm cA, cB;    // content of A and B
     2455  CanonicalForm ppA, ppB;    // primitive part of A and B
     2456  CanonicalForm gcdcAcB;
     2457  cA = uni_content (A);
     2458  cB = uni_content (B);
     2459  gcdcAcB= gcd (cA, cB);
     2460  ppA= A/cA;
     2461  ppB= B/cB;
     2462
     2463  CanonicalForm lcA, lcB;  // leading coefficients of A and B
     2464  CanonicalForm gcdlcAlcB;
     2465  lcA= uni_lcoeff (ppA);
     2466  lcB= uni_lcoeff (ppB);
     2467
     2468  if (fdivides (lcA, lcB))
     2469  {
     2470    if (fdivides (A, B))
     2471      return F/Lc(F);
     2472  }   
     2473  if (fdivides (lcB, lcA))
     2474  {
     2475    if (fdivides (B, A))
     2476      return G/Lc(G);
     2477  }
     2478
     2479  gcdlcAlcB= gcd (lcA, lcB);
     2480  int skelSize= size (skel, skel.mvar());
     2481
     2482  int j= 0;
     2483  int biggestSize= 0;
     2484  int bufSize;
     2485  int numberUni= 0;
     2486  for (CFIterator i= skel; i.hasTerms(); i++, j++)
     2487  {
     2488    bufSize= size (i.coeff());
     2489    biggestSize= tmax (biggestSize, bufSize);
     2490    numberUni += bufSize;
     2491  }
     2492  numberUni--;
     2493  numberUni= (int) ceil ( (double) numberUni / (skelSize - 1));
     2494  biggestSize= tmax (biggestSize , numberUni);
     2495
     2496  numberUni= biggestSize + size (LC(skel)) - 1;
     2497  int biggestSize2= tmax (numberUni, biggestSize);
     2498
     2499  CanonicalForm g, Aeval, Beval;
     2500
     2501  CFList evalPoints;
     2502  CFArray coeffEval;
     2503  bool evalFail= false;
     2504  CFList list;
     2505  bool GF= false;
     2506  CanonicalForm LCA= LC (A);
     2507  CanonicalForm tmp;
     2508  CFArray gcds= CFArray (biggestSize);
     2509  CFList * pEvalPoints= new CFList [biggestSize];
     2510  Variable V_buf= alpha;
     2511  CFList source, dest;
     2512  CanonicalForm prim_elem, im_prim_elem;
     2513  for (int i= 0; i < biggestSize; i++)
     2514  {
     2515    if (i == 0)
     2516    {
     2517      if (getCharacteristic() > 3)
     2518        evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf,
     2519                                    evalFail, list);
     2520      else
     2521        evalFail= true;
     2522     
     2523      if (evalFail)
     2524      {
     2525        if (V_buf != Variable (1))
     2526        {
     2527          do
     2528          {
     2529            int num_vars= tmin (getNumVars(A), getNumVars(B));
     2530            int d= totaldegree (A, Variable(2), Variable (A.level()));
     2531            d= tmin (d, totaldegree (B, Variable(2), Variable (B.level())));
     2532            Variable V_buf2= V_buf;
     2533            choose_extension (d, num_vars, V_buf2);
     2534            source= CFList();
     2535            dest= CFList();
     2536         
     2537            bool prim_fail= false;
     2538            Variable V_buf3;
     2539            prim_elem= primitiveElement (V_buf, V_buf3, prim_fail);
     2540       
     2541            ASSERT (!prim_fail, "failure in integer factorizer");
     2542            if (prim_fail)
     2543              ; //ERROR
     2544            else
     2545              im_prim_elem= mapPrimElem (prim_elem, V_buf, V_buf2);
     2546
     2547            DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf));
     2548            DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf2));
     2549
     2550            for (CFListIterator i= list; i.hasItem(); i++)
     2551              i.getItem()= mapUp (i.getItem(), V_buf, V_buf2, prim_elem,
     2552                                im_prim_elem, source, dest);
     2553            evalFail= false;
     2554            evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf,
     2555                                          evalFail, list);
     2556          } while (evalFail);
     2557        }
     2558        else
     2559        {
     2560          CanonicalForm mipo;
     2561          int deg= 2;
     2562          do {
     2563            mipo= randomIrredpoly (deg, x);
     2564            V_buf= rootOf (mipo);
     2565            evalFail= false;
     2566            evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf,
     2567                                          evalFail, list);
     2568            deg++;
     2569          } while (evalFail);
     2570        }
     2571      }
     2572    }
     2573    else
     2574    {
     2575      mult (evalPoints, pEvalPoints[0]);
     2576      eval (A, B, Aeval, Beval, evalPoints);
     2577    }
     2578
     2579    g= gcd (Aeval, Beval);
     2580    g /= Lc (g);
     2581
     2582    if (degree (g) != degree (skel) || (skelSize != size (g)))
     2583    {
     2584      delete[] pEvalPoints;
     2585      fail= true;
     2586      return 0;
     2587    }
     2588    CFIterator l= skel;
     2589    for (CFIterator k= g; k.hasTerms(); k++, l++)
     2590    {
     2591      if (k.exp() != l.exp())
     2592      {
     2593        delete[] pEvalPoints;
     2594        fail= true;
     2595        return 0;
     2596      }
     2597    }
     2598    pEvalPoints[i]= evalPoints;
     2599    gcds[i]= g;
     2600
     2601    tmp= 0;
     2602    int j= 0;
     2603    for (CFListIterator k= evalPoints; k.hasItem(); k++, j++)
     2604      tmp += k.getItem()*power (x, j);
     2605    list.append (tmp);
     2606  }
     2607
     2608  if (Monoms.size() == 0)
     2609    Monoms= getMonoms (skel);
     2610
     2611  if (coeffMonoms == NULL)
     2612    coeffMonoms= new CFArray [skelSize];
     2613
     2614  j= 0;
     2615  for (CFIterator i= skel; i.hasTerms(); i++, j++)
     2616    coeffMonoms[j]= getMonoms (i.coeff());
     2617   
     2618  int minimalColumnsIndex;
     2619  if (skelSize > 1)
     2620    minimalColumnsIndex= 1;
     2621  else
     2622    minimalColumnsIndex= 0;
     2623  int minimalColumns;
     2624
     2625  CFArray* pM= new CFArray [skelSize];
     2626  CFMatrix Mat;
     2627  // find the Matrix with minimal number of columns
     2628  for (int i= 0; i < skelSize; i++)
     2629  {
     2630    pM[i]= CFArray (coeffMonoms[i].size());
     2631    if (i == 1)
     2632      minimalColumns= coeffMonoms[i].size();
     2633    if (i > 1)
     2634    {
     2635      if (minimalColumns > coeffMonoms[i].size())
     2636      {
     2637        minimalColumns= coeffMonoms[i].size();
     2638        minimalColumnsIndex= i;
     2639      }
     2640    } 
     2641  }
     2642  CFMatrix* pMat= new CFMatrix [2];
     2643  pMat[0]= CFMatrix (biggestSize, coeffMonoms[0].size());
     2644  pMat[1]= CFMatrix (biggestSize, coeffMonoms[minimalColumnsIndex].size());
     2645  CFArray* pL= new CFArray [skelSize];
     2646  for (int i= 0; i < biggestSize; i++)
     2647  {
     2648    CFIterator l= gcds [i];
     2649    evalPoints= pEvalPoints [i];
     2650    for (int k= 0; k < skelSize; k++, l++)
     2651    {
     2652      if (i == 0)
     2653        pL[k]= CFArray (biggestSize);
     2654      pL[k] [i]= l.coeff();
     2655
     2656      if (i == 0 && k != 0 && k != minimalColumnsIndex)
     2657        pM[k]= evaluate (coeffMonoms[k], evalPoints);
     2658      else if (i == 0 && (k == 0 || k == minimalColumnsIndex))
     2659        coeffEval= evaluate (coeffMonoms[k], evalPoints);
     2660      if (i > 0 && (k == 0 || k == minimalColumnsIndex))
     2661        coeffEval= evaluate (coeffMonoms[k], evalPoints);
     2662     
     2663      if (k == 0)
     2664      {
     2665        if (pMat[k].rows() >= i + 1)
     2666        {
     2667          for (int ind= 1; ind < coeffEval.size() + 1; ind++)
     2668            pMat[k] (i + 1, ind)= coeffEval[ind - 1];                 
     2669        }
     2670      }
     2671      if (k == minimalColumnsIndex)
     2672      {
     2673        if (pMat[1].rows() >= i + 1)
     2674        {
     2675          for (int ind= 1; ind < coeffEval.size() + 1; ind++)
     2676            pMat[1] (i + 1, ind)= coeffEval[ind - 1];                 
     2677        }
     2678      }
     2679    }
     2680  }
     2681 
     2682  CFArray solution;
     2683  CanonicalForm result= 0;
     2684  int ind= 1;
     2685  int matRows, matColumns;
     2686  matRows= pMat[1].rows();
     2687  matColumns= pMat[0].columns() - 1;
     2688  matColumns += pMat[1].columns();
     2689
     2690  Mat= CFMatrix (matRows, matColumns);
     2691  for (int i= 1; i <= matRows; i++)
     2692    for (int j= 1; j <= pMat[1].columns(); j++)
     2693      Mat (i, j)= pMat[1] (i, j);
     2694
     2695  for (int j= pMat[1].columns() + 1; j <= matColumns;
     2696       j++, ind++)
     2697  {
     2698    for (int i= 1; i <= matRows; i++)
     2699      Mat (i, j)= (-pMat [0] (i, ind + 1))*pL[minimalColumnsIndex] [i - 1];
     2700  }
     2701
     2702  CFArray firstColumn= CFArray (pMat[0].rows());
     2703  for (int i= 0; i < pMat[0].rows(); i++)
     2704    firstColumn [i]= pMat[0] (i + 1, 1);
     2705
     2706  CFArray bufArray= pL[minimalColumnsIndex];
     2707 
     2708  for (int i= 0; i < biggestSize; i++)
     2709    pL[minimalColumnsIndex] [i] *= firstColumn[i];
     2710
     2711  CFMatrix bufMat= pMat[1];
     2712  pMat[1]= Mat;
     2713
     2714  if (V_buf != x)
     2715    solution= solveSystemFq (pMat[1],
     2716                             pL[minimalColumnsIndex], V_buf);
     2717  else
     2718    solution= solveSystemFp (pMat[1],
     2719                             pL[minimalColumnsIndex]);
     2720                         
     2721  if (solution.size() == 0) //Javadi's method failed, try deKleine, Monagan, Wittkopf instead
     2722  {
     2723    CFMatrix bufMat0= pMat[0];
     2724    delete [] pMat;
     2725    pMat= new CFMatrix [skelSize];
     2726    pL[minimalColumnsIndex]= bufArray;
     2727    CFList* bufpEvalPoints;
     2728    CFArray bufGcds;
     2729    if (biggestSize != biggestSize2)
     2730    {
     2731      bufpEvalPoints= pEvalPoints;
     2732      pEvalPoints= new CFList [biggestSize2];
     2733      bufGcds= gcds;
     2734      gcds= CFArray (biggestSize2);
     2735      for (int i= 0; i < biggestSize; i++)
     2736      {
     2737        pEvalPoints[i]= bufpEvalPoints [i];
     2738        gcds[i]= bufGcds[i];
     2739      } 
     2740      for (int i= 0; i < biggestSize2 - biggestSize; i++)
     2741      {
     2742        mult (evalPoints, pEvalPoints[0]);
     2743        eval (A, B, Aeval, Beval, evalPoints);
     2744        g= gcd (Aeval, Beval);
     2745        g /= Lc (g);
     2746        if (degree (g) != degree (skel) || (skelSize != size (g)))
     2747        {
     2748          delete[] pEvalPoints;
     2749          delete[] pMat;
     2750          delete[] pL;
     2751          delete[] coeffMonoms;
     2752          delete[] pM;
     2753          fail= true;
     2754          return 0;
     2755        }
     2756        CFIterator l= skel;
     2757        for (CFIterator k= g; k.hasTerms(); k++, l++)
     2758        {
     2759          if (k.exp() != l.exp())
     2760          {
     2761            delete[] pEvalPoints;
     2762            delete[] pMat;
     2763            delete[] pL;
     2764            delete[] coeffMonoms;
     2765            delete[] pM;
     2766            fail= true;
     2767            return 0;
     2768          }
     2769        }
     2770        pEvalPoints[i + biggestSize]= evalPoints;
     2771        gcds[i + biggestSize]= g;
     2772      }
     2773    }
     2774      for (int i= 0; i < biggestSize; i++)
     2775      {
     2776        CFIterator l= gcds [i];
     2777        evalPoints= pEvalPoints [i];
     2778        for (int k= 1; k < skelSize; k++, l++)
     2779        {
     2780          if (i == 0)
     2781            pMat[k]= CFMatrix (biggestSize2, coeffMonoms[k].size() + biggestSize2 - 1);
     2782          if (k == minimalColumnsIndex)
     2783            continue;
     2784          coeffEval= evaluate (coeffMonoms[k], evalPoints);
     2785          if (pMat[k].rows() >= i + 1)
     2786          {
     2787            for (int ind= 1; ind < coeffEval.size() + 1; ind++)
     2788              pMat[k] (i + 1, ind)= coeffEval[ind - 1];
     2789          }
     2790        }
     2791      }
     2792      Mat= bufMat0;
     2793      pMat[0]= CFMatrix (biggestSize2, coeffMonoms[0].size() + biggestSize2 - 1);
     2794      for (int j= 1; j <= Mat.rows(); j++)
     2795        for (int k= 1; k <= Mat.columns(); k++)
     2796          pMat [0] (j,k)= Mat (j,k);
     2797      Mat= bufMat;
     2798      for (int j= 1; j <= Mat.rows(); j++)
     2799        for (int k= 1; k <= Mat.columns(); k++)
     2800          pMat [minimalColumnsIndex] (j,k)= Mat (j,k);
     2801      // write old matrix entries into new matrices
     2802      for (int i= 0; i < skelSize; i++)
     2803      {
     2804        bufArray= pL[i];
     2805        pL[i]= CFArray (biggestSize2);
     2806        for (int j= 0; j < bufArray.size(); j++)
     2807          pL[i] [j]= bufArray [j];
     2808      }
     2809      //write old vector entries into new and add new entries to old matrices
     2810      for (int i= 0; i < biggestSize2 - biggestSize; i++)
     2811      {
     2812        CFIterator l= gcds [i + biggestSize];
     2813        evalPoints= pEvalPoints [i + biggestSize];
     2814        for (int k= 0; k < skelSize; k++, l++)
     2815        {
     2816          pL[k] [i + biggestSize]= l.coeff();
     2817          coeffEval= evaluate (coeffMonoms[k], evalPoints);
     2818          if (pMat[k].rows() >= i + biggestSize + 1)
     2819          {
     2820            for (int ind= 1; ind < coeffEval.size() + 1; ind++)
     2821              pMat[k] (i + biggestSize + 1, ind)= coeffEval[ind - 1];
     2822          }
     2823        }
     2824      }
     2825      // begin new
     2826      for (int i= 0; i < skelSize; i++)
     2827      {
     2828        if (pL[i].size() > 1)
     2829        {
     2830          for (int j= 2; j <= pMat[i].rows(); j++)
     2831            pMat[i] (j, coeffMonoms[i].size() + j - 1)=
     2832               -pL[i] [j - 1];
     2833        }
     2834      }
     2835   
     2836   
     2837    long rk;
     2838    matColumns= biggestSize2 - 1;
     2839    matRows= 0;
     2840    for (int i= 0; i < skelSize; i++)
     2841    {
     2842      if (V_buf == x)
     2843        rk= gaussianElimFp (pMat[i], pL[i]);
     2844      else
     2845        rk= gaussianElimFq (pMat[i], pL[i], V_buf);
     2846
     2847      if (pMat[i] (coeffMonoms[i].size(), coeffMonoms[i].size()) == 0)
     2848      {
     2849        delete[] pEvalPoints;
     2850        delete[] pMat;
     2851        delete[] pL;
     2852        delete[] coeffMonoms;
     2853        delete[] pM;
     2854        fail= true;
     2855        return 0;
     2856      }
     2857      matRows += pMat[i].rows() - coeffMonoms[i].size();
     2858    }
     2859
     2860    CFMatrix bufMat;
     2861    Mat= CFMatrix (matRows, matColumns);
     2862    ind= 0;
     2863    bufArray= CFArray (matRows);
     2864    CFArray bufArray2;
     2865    for (int i= 0; i < skelSize; i++)
     2866    {
     2867      bufMat= pMat[i] (coeffMonoms[i].size() + 1, pMat[i].rows(),
     2868                       coeffMonoms[i].size() + 1, pMat[i].columns());
     2869
     2870      for (int j= 1; j <= bufMat.rows(); j++)
     2871        for (int k= 1; k <= bufMat.columns(); k++)
     2872          Mat (j + ind, k)= bufMat(j, k);
     2873      bufArray2= coeffMonoms[i].size();
     2874      for (int j= 1; j <= pMat[i].rows(); j++)
     2875      {
     2876        if (j > coeffMonoms[i].size())
     2877          bufArray [j-coeffMonoms[i].size() + ind - 1]= pL[i] [j - 1];
     2878        else
     2879          bufArray2 [j - 1]= pL[i] [j - 1];
     2880      }
     2881      pL[i]= bufArray2;
     2882      ind += bufMat.rows();     
     2883      pMat[i]= pMat[i] (1, coeffMonoms[i].size(), 1, pMat[i].columns());
     2884    }
     2885   
     2886    if (V_buf != x)
     2887      solution= solveSystemFq (Mat, bufArray, V_buf);
     2888    else
     2889      solution= solveSystemFp (Mat, bufArray);
     2890   
     2891    if (solution.size() == 0)
     2892    {
     2893      delete[] pEvalPoints;
     2894      delete[] pMat;
     2895      delete[] pL;
     2896      delete[] coeffMonoms;
     2897      delete[] pM;
     2898      fail= true;
     2899      return 0;
     2900    }
     2901   
     2902    ind= 0;
     2903    result= 0;
     2904    CFArray bufSolution;
     2905    for (int i= 0; i < skelSize; i++)
     2906    {
     2907      bufSolution= readOffSolution (pMat[i], pL[i], solution);
     2908      for (int i= 0; i < bufSolution.size(); i++, ind++)
     2909        result += Monoms [ind]*bufSolution[i];
     2910    }
     2911    if (alpha != Variable (1) && V_buf != alpha)
     2912    {
     2913      CFList u, v;
     2914      result= mapDown (result, prim_elem, im_prim_elem, alpha, u, v);
     2915    }
     2916    result= N(result);
     2917    if (fdivides (result, F) && fdivides (result, G))
     2918      return result;
     2919    else
     2920    {
     2921      fail= true;
     2922      return 0;
     2923    }
     2924  } // end of deKleine, Monagan & Wittkopf 
     2925
     2926  result += Monoms[0];
     2927  int ind2= 0, ind3= 2;
     2928  ind= 0;
     2929  for (int l= 0; l < minimalColumnsIndex; l++)
     2930    ind += coeffMonoms[l].size();
     2931  for (int l= solution.size() - pMat[0].columns() + 1; l < solution.size();
     2932       l++, ind2++, ind3++)
     2933  {
     2934    result += solution[l]*Monoms [1 + ind2];
     2935    for (int i= 0; i < pMat[0].rows(); i++)
     2936      firstColumn[i] += solution[l]*pMat[0] (i + 1, ind3);
     2937  }
     2938  for (int l= 0; l < solution.size() + 1 - pMat[0].columns(); l++)
     2939    result += solution[l]*Monoms [ind + l];
     2940  ind= coeffMonoms[0].size();
     2941  for (int k= 1; k < skelSize; k++)
     2942  {
     2943    if (k == minimalColumnsIndex)
     2944    {
     2945      ind += coeffMonoms[k].size();
     2946      continue;
     2947    }
     2948    if (k != minimalColumnsIndex)
     2949    {
     2950      for (int i= 0; i < biggestSize; i++)
     2951        pL[k] [i] *= firstColumn [i];
     2952    }
     2953
     2954    if (biggestSize != coeffMonoms[k].size() && k != minimalColumnsIndex)
     2955    {
     2956      bufArray= CFArray (coeffMonoms[k].size());
     2957      for (int i= 0; i < bufArray.size(); i++)
     2958        bufArray [i]= pL[k] [i];
     2959      solution= solveGeneralVandermonde (pM[k], bufArray);
     2960    }
     2961    else
     2962      solution= solveGeneralVandermonde (pM[k], pL[k]);
     2963
     2964    if (solution.size() == 0)
     2965    {
     2966      delete[] pEvalPoints;
     2967      delete[] pMat;
     2968      delete[] pL;
     2969      delete[] coeffMonoms;
     2970      delete[] pM;
     2971      fail= true;
     2972      return 0;
     2973    }
     2974    if (k != minimalColumnsIndex)
     2975    {
     2976      for (int l= 0; l < solution.size(); l++)
     2977        result += solution[l]*Monoms [ind + l];
     2978      ind += solution.size();
     2979    }
     2980  }
     2981
     2982  delete[] pEvalPoints;
     2983  delete[] pMat;
     2984  delete[] pL;
     2985  delete[] pM;
     2986 
     2987  if (alpha != Variable (1) && V_buf != alpha)
     2988  {
     2989    CFList u, v;
     2990    result= mapDown (result, prim_elem, im_prim_elem, alpha, u, v);
     2991  }
     2992  result= N(result);
     2993
     2994  if (fdivides (result, F) && fdivides (result, G))
     2995    return result;
     2996  else
     2997  {
     2998    delete[] coeffMonoms;
     2999    fail= true;
     3000    return 0;
     3001  }
     3002}
     3003
     3004CanonicalForm sparseGCDFq (const CanonicalForm& F, const CanonicalForm& G,
     3005                           const Variable & alpha, CFList& l, bool& topLevel)
     3006{
     3007  CanonicalForm A= F;
     3008  CanonicalForm B= G;
     3009  if (F.isZero() && degree(G) > 0) return G/Lc(G);
     3010  else if (G.isZero() && degree (F) > 0) return F/Lc(F);
     3011  else if (F.isZero() && G.isZero()) return F.genOne();
     3012  if (F.inBaseDomain() || G.inBaseDomain()) return F.genOne();
     3013  if (F.isUnivariate() && fdivides(F, G)) return F/Lc(F);
     3014  if (G.isUnivariate() && fdivides(G, F)) return G/Lc(G);
     3015  if (F == G) return F/Lc(F);
     3016
     3017  CFMap M,N;
     3018  int best_level= myCompress (A, B, M, N, topLevel);
     3019
     3020  if (best_level == 0) return B.genOne();
     3021
     3022  A= M(A);
     3023  B= M(B);
     3024
     3025  Variable x= Variable (1);
     3026
     3027  //univariate case
     3028  if (A.isUnivariate() && B.isUnivariate())
     3029    return N (gcd (A, B));
     3030
     3031  int substitute= substituteCheck (A, B);
     3032 
     3033  if (substitute > 1)
     3034    subst (A, B, A, B, substitute);
     3035   
     3036  CanonicalForm cA, cB;    // content of A and B
     3037  CanonicalForm ppA, ppB;    // primitive part of A and B
     3038  CanonicalForm gcdcAcB;
     3039  if (topLevel)
     3040  {
     3041    if (best_level <= 2)
     3042      gcdcAcB= extractContents (A, B, cA, cB, ppA, ppB, best_level);
     3043    else
     3044      gcdcAcB= extractContents (A, B, cA, cB, ppA, ppB, 2);
     3045  }
     3046  else
     3047  {
     3048    cA = uni_content (A);
     3049    cB = uni_content (B);
     3050    gcdcAcB= gcd (cA, cB);
     3051    ppA= A/cA;
     3052    ppB= B/cB;
     3053  }
     3054
     3055  CanonicalForm lcA, lcB;  // leading coefficients of A and B
     3056  CanonicalForm gcdlcAlcB;
     3057  lcA= uni_lcoeff (ppA);
     3058  lcB= uni_lcoeff (ppB);
     3059
     3060  if (fdivides (lcA, lcB))
     3061  {
     3062    if (fdivides (A, B))
     3063      return F/Lc(F);
     3064  }   
     3065  if (fdivides (lcB, lcA))
     3066  {
     3067    if (fdivides (B, A))
     3068      return G/Lc(G);
     3069  }
     3070
     3071  gcdlcAlcB= gcd (lcA, lcB);
     3072 
     3073  int d= totaldegree (ppA, Variable (2), Variable (ppA.level()));
     3074  int d0;
     3075
     3076  if (d == 0)
     3077  {
     3078    if (substitute > 1)
     3079      return N(reverseSubst (gcdcAcB, substitute));
     3080    else
     3081      return N(gcdcAcB);
     3082  }
     3083  d0= totaldegree (ppB, Variable (2), Variable (ppB.level()));
     3084
     3085  if (d0 < d)
     3086    d= d0;
     3087
     3088  if (d == 0)
     3089  {
     3090    if (substitute > 1)
     3091      return N(reverseSubst (gcdcAcB, substitute));
     3092    else
     3093      return N(gcdcAcB);
     3094  }
     3095
     3096  CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH, skeleton;
     3097  CanonicalForm newtonPoly= 1;
     3098  m= gcdlcAlcB;
     3099  G_m= 0;
     3100  H= 0;
     3101  bool fail= false;
     3102  topLevel= false;
     3103  bool inextension= false;
     3104  Variable V_buf= alpha;
     3105  CanonicalForm prim_elem, im_prim_elem;
     3106  CFList source, dest;
     3107  do // first do
     3108  {
     3109    random_element= randomElement (m, V_buf, l, fail);
     3110    if (random_element == 0 && !fail)
     3111    {
     3112      if (!find (l, random_element))
     3113        l.append (random_element);
     3114      continue;
     3115    }
     3116    if (fail)
     3117    {
     3118      source= CFList();
     3119      dest= CFList();
     3120      int num_vars= tmin (getNumVars(A), getNumVars(B));
     3121      num_vars--;
     3122
     3123      choose_extension (d, num_vars, V_buf);
     3124      bool prim_fail= false;
     3125      Variable V_buf2;
     3126      prim_elem= primitiveElement (alpha, V_buf2, prim_fail);
     3127
     3128      ASSERT (!prim_fail, "failure in integer factorizer");
     3129      if (prim_fail)
     3130        ; //ERROR
     3131      else
     3132        im_prim_elem= mapPrimElem (prim_elem, alpha, V_buf);
     3133
     3134      DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (alpha));
     3135      DEBOUTLN (cerr, "getMipo (V_buf2)= " << getMipo (V_buf2));
     3136      inextension= true;
     3137      for (CFListIterator i= l; i.hasItem(); i++)
     3138        i.getItem()= mapUp (i.getItem(), alpha, V_buf, prim_elem,
     3139                             im_prim_elem, source, dest);
     3140      m= mapUp (m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
     3141      G_m= mapUp (G_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
     3142      newtonPoly= mapUp (newtonPoly, alpha, V_buf, prim_elem, im_prim_elem,
     3143                          source, dest);
     3144      ppA= mapUp (ppA, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
     3145      ppB= mapUp (ppB, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
     3146      gcdlcAlcB= mapUp (gcdlcAlcB, alpha, V_buf, prim_elem, im_prim_elem,
     3147                         source, dest);
     3148
     3149      fail= false;
     3150      random_element= randomElement (m, V_buf, l, fail );
     3151      DEBOUTLN (cerr, "fail= " << fail);
     3152      CFList list;
     3153      TIMING_START (gcd_recursion);
     3154      G_random_element=
     3155      sparseGCDFq (ppA (random_element, x), ppB (random_element, x), V_buf,
     3156                        list, topLevel);
     3157      TIMING_END_AND_PRINT (gcd_recursion,
     3158                            "time for recursive call: ");
     3159      DEBOUTLN (cerr, "G_random_element= " << G_random_element);
     3160    }
     3161    else
     3162    {
     3163      CFList list;
     3164      TIMING_START (gcd_recursion);
     3165      G_random_element=
     3166      sparseGCDFq (ppA(random_element, x), ppB(random_element, x), V_buf,
     3167                        list, topLevel);
     3168      TIMING_END_AND_PRINT (gcd_recursion,
     3169                            "time for recursive call: ");
     3170      DEBOUTLN (cerr, "G_random_element= " << G_random_element);
     3171    }
     3172
     3173    d0= totaldegree (G_random_element, Variable(2),
     3174                     Variable (G_random_element.level()));
     3175    if (d0 == 0) 
     3176    {
     3177      if (substitute > 1)
     3178        return N(reverseSubst (gcdcAcB, substitute));
     3179      else
     3180        return N(gcdcAcB);
     3181    }
     3182    if (d0 >  d)
     3183    {
     3184      if (!find (l, random_element))
     3185        l.append (random_element);
     3186      continue;
     3187    }
     3188
     3189    G_random_element=
     3190    (gcdlcAlcB(random_element, x)/uni_lcoeff (G_random_element))
     3191    * G_random_element;
     3192   
     3193    skeleton= G_random_element;
     3194    d0= totaldegree (G_random_element, Variable(2),
     3195                     Variable(G_random_element.level()));
     3196    if (d0 <  d)
     3197    {
     3198      m= gcdlcAlcB;
     3199      newtonPoly= 1;
     3200      G_m= 0;
     3201      d= d0;
     3202    }
     3203
     3204    H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
     3205    if (uni_lcoeff (H) == gcdlcAlcB)
     3206    {
     3207      cH= uni_content (H);
     3208      ppH= H/cH;
     3209      if (inextension)
     3210      {
     3211        CFList u, v;
     3212        //maybe it's better to test if ppH is an element of F(\alpha) before
     3213        //mapping down
     3214        DEBOUTLN (cerr, "ppH before mapDown= " << ppH);
     3215        ppH= mapDown (ppH, prim_elem, im_prim_elem, alpha, u, v);
     3216        ppH /= Lc(ppH);
     3217        DEBOUTLN (cerr, "ppH after mapDown= " << ppH);
     3218        if (fdivides (ppH, A) && fdivides (ppH, B))
     3219        {
     3220          if (substitute > 1)
     3221          {
     3222            ppH= reverseSubst (ppH, substitute);
     3223            gcdcAcB= reverseSubst (gcdcAcB, substitute);
     3224          }
     3225          return N(gcdcAcB*ppH);
     3226        }
     3227      }
     3228      else if (fdivides (ppH, A) && fdivides (ppH, B))
     3229      {
     3230        if (substitute > 1)
     3231        {
     3232          ppH= reverseSubst (ppH, substitute);
     3233          gcdcAcB= reverseSubst (gcdcAcB, substitute);
     3234        }
     3235        return N(gcdcAcB*ppH);
     3236      }
     3237    }
     3238    G_m= H;
     3239    newtonPoly= newtonPoly*(x - random_element);
     3240    m= m*(x - random_element);
     3241    if (!find (l, random_element))
     3242      l.append (random_element);
     3243     
     3244    if (getCharacteristic () > 3 && size (skeleton) < 100) // TODO consider a ratio of characteristic and degree?
     3245    {
     3246    CFArray Monoms;
     3247    CFArray *coeffMonoms= NULL;
     3248    do //second do
     3249    {
     3250      random_element= randomElement (m, V_buf, l, fail);
     3251      if (random_element == 0 && !fail)
     3252      {
     3253        if (!find (l, random_element))
     3254          l.append (random_element);
     3255        continue;
     3256      }
     3257      if (fail)
     3258      {
     3259        source= CFList();
     3260        dest= CFList();
     3261        int num_vars= tmin (getNumVars(A), getNumVars(B));
     3262        num_vars--;
     3263
     3264        choose_extension (d, num_vars, V_buf);
     3265        bool prim_fail= false;
     3266        Variable V_buf2;
     3267        prim_elem= primitiveElement (alpha, V_buf2, prim_fail);
     3268
     3269        ASSERT (!prim_fail, "failure in integer factorizer");
     3270        if (prim_fail)
     3271          ; //ERROR
     3272        else
     3273          im_prim_elem= mapPrimElem (prim_elem, alpha, V_buf);
     3274
     3275        DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (alpha));
     3276        DEBOUTLN (cerr, "getMipo (V_buf2)= " << getMipo (V_buf2));
     3277        inextension= true;
     3278        for (CFListIterator i= l; i.hasItem(); i++)
     3279          i.getItem()= mapUp (i.getItem(), alpha, V_buf, prim_elem,
     3280                              im_prim_elem, source, dest);
     3281        m= mapUp (m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
     3282        G_m= mapUp (G_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
     3283        newtonPoly= mapUp (newtonPoly, alpha, V_buf, prim_elem, im_prim_elem,
     3284                            source, dest);
     3285        ppA= mapUp (ppA, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
     3286        ppB= mapUp (ppB, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
     3287
     3288        gcdlcAlcB= mapUp (gcdlcAlcB, alpha, V_buf, prim_elem, im_prim_elem,
     3289                          source, dest);
     3290
     3291        fail= false;
     3292        random_element= randomElement (m, V_buf, l, fail);
     3293        DEBOUTLN (cerr, "fail= " << fail);
     3294        CFList list;
     3295        TIMING_START (gcd_recursion);
     3296
     3297        //sparseInterpolation
     3298        bool sparseFail= false;
     3299        if (LC (skeleton).inCoeffDomain())
     3300          G_random_element=
     3301          monicSparseInterpol (ppA (random_element, x), ppB (random_element, x),
     3302                               skeleton,V_buf, sparseFail, coeffMonoms, Monoms);
     3303        else
     3304          G_random_element=
     3305          nonMonicSparseInterpol (ppA(random_element,x),ppB (random_element,x),
     3306                                  skeleton, V_buf, sparseFail, coeffMonoms,
     3307                                  Monoms);
     3308        if (sparseFail)
     3309          break;
     3310       
     3311        TIMING_END_AND_PRINT (gcd_recursion,
     3312                              "time for recursive call: ");
     3313        DEBOUTLN (cerr, "G_random_element= " << G_random_element);
     3314      }
     3315      else
     3316      {
     3317        CFList list;
     3318        TIMING_START (gcd_recursion);
     3319        bool sparseFail= false;
     3320        if (LC (skeleton).inCoeffDomain())
     3321          G_random_element=
     3322          monicSparseInterpol (ppA (random_element, x), ppB (random_element, x),
     3323                               skeleton,V_buf, sparseFail, coeffMonoms, Monoms);
     3324        else
     3325          G_random_element=
     3326          nonMonicSparseInterpol (ppA(random_element,x), ppB(random_element,x),
     3327                                  skeleton, V_buf, sparseFail, coeffMonoms,
     3328                                  Monoms);
     3329        if (sparseFail)
     3330          break;
     3331
     3332        TIMING_END_AND_PRINT (gcd_recursion,
     3333                              "time for recursive call: ");
     3334        DEBOUTLN (cerr, "G_random_element= " << G_random_element);
     3335      }
     3336
     3337      d0= totaldegree (G_random_element, Variable(2),
     3338                      Variable (G_random_element.level()));
     3339      if (d0 == 0) 
     3340      {
     3341        if (substitute > 1)
     3342          return N(reverseSubst (gcdcAcB, substitute));
     3343        else
     3344          return N(gcdcAcB);
     3345      }
     3346      if (d0 >  d)
     3347      {
     3348        if (!find (l, random_element))
     3349          l.append (random_element);
     3350        continue;
     3351      }
     3352
     3353      G_random_element=
     3354      (gcdlcAlcB(random_element, x)/uni_lcoeff (G_random_element))
     3355      * G_random_element;
     3356
     3357      d0= totaldegree (G_random_element, Variable(2),
     3358                      Variable(G_random_element.level()));
     3359      if (d0 <  d)
     3360      {
     3361        m= gcdlcAlcB;
     3362        newtonPoly= 1;
     3363        G_m= 0;
     3364        d= d0;
     3365      }
     3366
     3367      TIMING_START (newton_interpolation);
     3368      H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
     3369      TIMING_END_AND_PRINT (newton_interpolation,
     3370                            "time for newton interpolation: ");
     3371
     3372      //termination test
     3373      if (uni_lcoeff (H) == gcdlcAlcB)
     3374      {
     3375        cH= uni_content (H);
     3376        ppH= H/cH;
     3377        if (inextension)
     3378        {
     3379          CFList u, v;
     3380          //maybe it's better to test if ppH is an element of F(\alpha) before
     3381          //mapping down
     3382          DEBOUTLN (cerr, "ppH before mapDown= " << ppH);
     3383          ppH= mapDown (ppH, prim_elem, im_prim_elem, alpha, u, v);
     3384          ppH /= Lc(ppH);
     3385          DEBOUTLN (cerr, "ppH after mapDown= " << ppH);
     3386          if (fdivides (ppH, A) && fdivides (ppH, B))
     3387          {
     3388            if (substitute > 1)
     3389            {
     3390              ppH= reverseSubst (ppH, substitute);
     3391              gcdcAcB= reverseSubst (gcdcAcB, substitute);
     3392            }
     3393            return N(gcdcAcB*ppH);
     3394          }
     3395        }
     3396        else if (fdivides (ppH, A) && fdivides (ppH, B))
     3397        {
     3398          if (substitute > 1)
     3399          {
     3400            ppH= reverseSubst (ppH, substitute);
     3401            gcdcAcB= reverseSubst (gcdcAcB, substitute);
     3402          }
     3403          return N(gcdcAcB*ppH);
     3404        }
     3405      }
     3406
     3407      G_m= H;
     3408      newtonPoly= newtonPoly*(x - random_element);
     3409      m= m*(x - random_element);
     3410      if (!find (l, random_element))
     3411        l.append (random_element);
     3412
     3413    } while (1);
     3414  }
     3415  } while (1);
     3416}
     3417
     3418CanonicalForm sparseGCDFp (const CanonicalForm& F, const CanonicalForm& G,
     3419                           bool& topLevel, CFList& l)
     3420{
     3421  CanonicalForm A= F;
     3422  CanonicalForm B= G;
     3423  if (F.isZero() && degree(G) > 0) return G/Lc(G);
     3424  else if (G.isZero() && degree (F) > 0) return F/Lc(F);
     3425  else if (F.isZero() && G.isZero()) return F.genOne();
     3426  if (F.inBaseDomain() || G.inBaseDomain()) return F.genOne();
     3427  if (F.isUnivariate() && fdivides(F, G)) return F/Lc(F);
     3428  if (G.isUnivariate() && fdivides(G, F)) return G/Lc(G);
     3429  if (F == G) return F/Lc(F);
     3430
     3431  CFMap M,N;
     3432  int best_level= myCompress (A, B, M, N, topLevel);
     3433
     3434  if (best_level == 0) return B.genOne();
     3435
     3436  A= M(A);
     3437  B= M(B);
     3438
     3439  Variable x= Variable (1);
     3440
     3441  //univariate case
     3442  if (A.isUnivariate() && B.isUnivariate())
     3443    return N (gcd (A, B));
     3444
     3445  int substitute= substituteCheck (A, B);
     3446 
     3447  if (substitute > 1)
     3448    subst (A, B, A, B, substitute);
     3449
     3450  CanonicalForm cA, cB;    // content of A and B
     3451  CanonicalForm ppA, ppB;    // primitive part of A and B
     3452  CanonicalForm gcdcAcB;
     3453  if (topLevel)
     3454  {
     3455    if (best_level <= 2)
     3456      gcdcAcB= extractContents (A, B, cA, cB, ppA, ppB, best_level);
     3457    else
     3458      gcdcAcB= extractContents (A, B, cA, cB, ppA, ppB, 2);
     3459  }
     3460  else
     3461  {
     3462    cA = uni_content (A);
     3463    cB = uni_content (B);
     3464    gcdcAcB= gcd (cA, cB);
     3465    ppA= A/cA;
     3466    ppB= B/cB;
     3467  }
     3468   
     3469  CanonicalForm lcA, lcB;  // leading coefficients of A and B
     3470  CanonicalForm gcdlcAlcB;
     3471  lcA= uni_lcoeff (ppA);
     3472  lcB= uni_lcoeff (ppB);
     3473
     3474  if (fdivides (lcA, lcB))
     3475  {
     3476    if (fdivides (A, B))
     3477      return F/Lc(F);
     3478  }   
     3479  if (fdivides (lcB, lcA))
     3480  {
     3481    if (fdivides (B, A))
     3482      return G/Lc(G);
     3483  }
     3484
     3485  gcdlcAlcB= gcd (lcA, lcB);
     3486 
     3487  int d= totaldegree (ppA, Variable (2), Variable (ppA.level()));
     3488  int d0;
     3489
     3490  if (d == 0)
     3491  {
     3492    if (substitute > 1)
     3493      return N(reverseSubst (gcdcAcB, substitute));
     3494    else
     3495      return N(gcdcAcB);
     3496  }
     3497  d0= totaldegree (ppB, Variable (2), Variable (ppB.level()));
     3498
     3499  if (d0 < d)
     3500    d= d0;
     3501
     3502  if (d == 0)
     3503  {
     3504    if (substitute > 1)
     3505      return N(reverseSubst (gcdcAcB, substitute));
     3506    else
     3507      return N(gcdcAcB);
     3508  }
     3509   
     3510  CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH, skeleton;
     3511  CanonicalForm newtonPoly= 1;
     3512  m= gcdlcAlcB;
     3513  G_m= 0;
     3514  H= 0;
     3515  bool fail= false;
     3516  bool topLevel2= topLevel;
     3517  int loops= 0;
     3518  topLevel= false;
     3519  bool inextension= false;
     3520  bool inextensionextension= false;
     3521  Variable V_buf, alpha;
     3522  CanonicalForm prim_elem, im_prim_elem;
     3523  CFList source, dest;
     3524  do //first do
     3525  {
     3526    if (inextension)
     3527      random_element= randomElement (m, alpha, l, fail);
     3528    else
     3529      random_element= FpRandomElement (m, l, fail);
     3530    if (random_element == 0 && !fail)
     3531    {
     3532      if (!find (l, random_element))
     3533        l.append (random_element);
     3534      continue;
     3535    }
     3536
     3537    if (!fail && !inextension)
     3538    {
     3539      CFList list;
     3540      TIMING_START (gcd_recursion);
     3541      G_random_element=
     3542      sparseGCDFp (ppA (random_element,x), ppB (random_element,x), topLevel,
     3543                   list);
     3544      TIMING_END_AND_PRINT (gcd_recursion,
     3545                            "time for recursive call: ");
     3546      DEBOUTLN (cerr, "G_random_element= " << G_random_element);
     3547    }
     3548    else if (!fail && inextension)
     3549    {
     3550      CFList list;
     3551      TIMING_START (gcd_recursion);
     3552      G_random_element=
     3553      sparseGCDFq (ppA (random_element, x), ppB (random_element, x), alpha,
     3554                        list, topLevel);
     3555      TIMING_END_AND_PRINT (gcd_recursion,
     3556                            "time for recursive call: ");
     3557      DEBOUTLN (cerr, "G_random_element= " << G_random_element);   
     3558    }
     3559    else if (fail && !inextension)
     3560    {
     3561      source= CFList();
     3562      dest= CFList();
     3563      CFList list;
     3564      CanonicalForm mipo;
     3565      int deg= 3;
     3566      do
     3567      {
     3568        mipo= randomIrredpoly (deg, x);
     3569        alpha= rootOf (mipo);
     3570        inextension= true;
     3571        fail= false;
     3572        random_element= randomElement (m, alpha, l, fail);
     3573        deg++;
     3574      } while (fail);
     3575      list= CFList();
     3576      TIMING_START (gcd_recursion);
     3577      G_random_element=
     3578      sparseGCDFq (ppA (random_element, x), ppB (random_element, x), alpha,
     3579                        list, topLevel);
     3580      TIMING_END_AND_PRINT (gcd_recursion,
     3581                            "time for recursive call: ");
     3582      DEBOUTLN (cerr, "G_random_element= " << G_random_element);
     3583    }
     3584    else if (fail && inextension)
     3585    {
     3586      source= CFList();
     3587      dest= CFList();
     3588      int num_vars= tmin (getNumVars(A), getNumVars(B));
     3589      num_vars--;
     3590      V_buf= alpha;
     3591      choose_extension (d, num_vars, V_buf);
     3592      bool prim_fail= false;
     3593      Variable V_buf2;
     3594      CanonicalForm prim_elem, im_prim_elem;
     3595      prim_elem= primitiveElement (alpha, V_buf2, prim_fail);
     3596       
     3597      ASSERT (!prim_fail, "failure in integer factorizer");
     3598      if (prim_fail)
     3599        ; //ERROR
     3600      else
     3601        im_prim_elem= mapPrimElem (prim_elem, alpha, V_buf);
     3602
     3603      DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (alpha));
     3604      DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf2));
     3605
     3606      inextensionextension= true;
     3607      for (CFListIterator i= l; i.hasItem(); i++)
     3608        i.getItem()= mapUp (i.getItem(), alpha, V_buf, prim_elem,
     3609                             im_prim_elem, source, dest);
     3610      m= mapUp (m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
     3611      G_m= mapUp (G_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
     3612      newtonPoly= mapUp (newtonPoly, alpha, V_buf, prim_elem, im_prim_elem,
     3613                          source, dest);
     3614      ppA= mapUp (ppA, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
     3615      ppB= mapUp (ppB, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
     3616      gcdlcAlcB= mapUp (gcdlcAlcB, alpha, V_buf, prim_elem, im_prim_elem,
     3617                         source, dest);
     3618      fail= false;
     3619      random_element= randomElement (m, V_buf, l, fail );
     3620      DEBOUTLN (cerr, "fail= " << fail);
     3621      CFList list;
     3622      TIMING_START (gcd_recursion);
     3623      G_random_element=
     3624      sparseGCDFq (ppA (random_element, x), ppB (random_element, x), V_buf,
     3625                        list, topLevel);
     3626      TIMING_END_AND_PRINT (gcd_recursion,
     3627                            "time for recursive call: ");
     3628      DEBOUTLN (cerr, "G_random_element= " << G_random_element);
     3629      alpha= V_buf;
     3630    }
     3631   
     3632    d0= totaldegree (G_random_element, Variable(2),
     3633                     Variable (G_random_element.level()));
     3634    if (d0 == 0) 
     3635    {
     3636      if (substitute > 1)
     3637        return N(reverseSubst (gcdcAcB, substitute));
     3638      else
     3639        return N(gcdcAcB);
     3640    }
     3641    if (d0 >  d)
     3642    {
     3643      if (!find (l, random_element))
     3644        l.append (random_element);
     3645      continue;
     3646    }
     3647
     3648    G_random_element=
     3649    (gcdlcAlcB(random_element, x)/uni_lcoeff (G_random_element))
     3650    * G_random_element;
     3651   
     3652    skeleton= G_random_element;
     3653   
     3654    if (size (skeleton) < 3 && getCharacteristic() <= 3)
     3655    {
     3656      return N (gcdcAcB*gcd_poly_p (ppA, ppB));
     3657    }
     3658   
     3659    d0= totaldegree (G_random_element, Variable(2),
     3660                     Variable(G_random_element.level()));
     3661    if (d0 <  d)
     3662    {
     3663      m= gcdlcAlcB;
     3664      newtonPoly= 1;
     3665      G_m= 0;
     3666      d= d0;
     3667    }
     3668
     3669    H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
     3670
     3671    if (uni_lcoeff (H) == gcdlcAlcB)
     3672    {
     3673      cH= uni_content (H);
     3674      ppH= H/cH;
     3675      ppH /= Lc (ppH);
     3676      DEBOUTLN (cerr, "ppH= " << ppH);
     3677     
     3678      if (fdivides (ppH, A) && fdivides (ppH, B))
     3679      {
     3680        if (substitute > 1)
     3681        {
     3682          ppH= reverseSubst (ppH, substitute);
     3683          gcdcAcB= reverseSubst (gcdcAcB, substitute);
     3684        }
     3685        return N(gcdcAcB*ppH);
     3686      }
     3687    }
     3688    G_m= H;
     3689    newtonPoly= newtonPoly*(x - random_element);
     3690    m= m*(x - random_element);
     3691    if (!find (l, random_element))
     3692      l.append (random_element);
     3693     
     3694    if (getCharacteristic () > 3 && size (skeleton) < 100) // TODO consider a ratio of characteristic and degree?
     3695    {
     3696    CFArray Monoms;
     3697    CFArray* coeffMonoms= NULL;
     3698
     3699    do //second do
     3700    {
     3701      if (inextension)
     3702        random_element= randomElement (m, alpha, l, fail);
     3703      else
     3704        random_element= FpRandomElement (m, l, fail);
     3705      if (random_element == 0 && !fail)
     3706      {
     3707        if (!find (l, random_element))
     3708          l.append (random_element);
     3709        continue;
     3710      }
     3711
     3712      bool sparseFail= false;
     3713      if (!fail && !inextension)
     3714      {
     3715        CFList list;
     3716        TIMING_START (gcd_recursion);
     3717        if (LC (skeleton).inCoeffDomain())
     3718          G_random_element=
     3719          monicSparseInterpol (ppA (random_element, x), ppB (random_element, x),
     3720                               skeleton, Variable(1), sparseFail, coeffMonoms,
     3721                               Monoms);
     3722        else
     3723          G_random_element=
     3724          nonMonicSparseInterpol (ppA(random_element,x), ppB(random_element,x),
     3725                                  skeleton, Variable (1), sparseFail,
     3726                                  coeffMonoms, Monoms);
     3727        TIMING_END_AND_PRINT (gcd_recursion,
     3728                              "time for recursive call: ");
     3729        DEBOUTLN (cerr, "G_random_element= " << G_random_element);
     3730      }
     3731      else if (!fail && inextension)
     3732      {
     3733        CFList list;
     3734        TIMING_START (gcd_recursion);
     3735        if (LC (skeleton).inCoeffDomain())
     3736          G_random_element=
     3737          monicSparseInterpol (ppA (random_element, x), ppB (random_element, x),
     3738                               skeleton, alpha, sparseFail, coeffMonoms,
     3739                               Monoms);
     3740        else
     3741          G_random_element=
     3742          nonMonicSparseInterpol (ppA(random_element,x), ppB(random_element,x),
     3743                                  skeleton, alpha, sparseFail, coeffMonoms,
     3744                                  Monoms);
     3745        TIMING_END_AND_PRINT (gcd_recursion,
     3746                              "time for recursive call: ");
     3747        DEBOUTLN (cerr, "G_random_element= " << G_random_element);   
     3748      }
     3749      else if (fail && !inextension)
     3750      {
     3751        source= CFList();
     3752        dest= CFList();
     3753        CFList list;
     3754        CanonicalForm mipo;
     3755        int deg= 3;
     3756        do
     3757        {
     3758          mipo= randomIrredpoly (deg, x);
     3759          alpha= rootOf (mipo);
     3760          inextension= true;
     3761          fail= false;
     3762          random_element= randomElement (m, alpha, l, fail);
     3763          deg++;
     3764        } while (fail);
     3765        list= CFList();
     3766        TIMING_START (gcd_recursion);
     3767        if (LC (skeleton).inCoeffDomain())
     3768          G_random_element=
     3769          monicSparseInterpol (ppA (random_element, x), ppB (random_element, x),
     3770                               skeleton, alpha, sparseFail, coeffMonoms,
     3771                               Monoms);
     3772        else
     3773          G_random_element=
     3774          nonMonicSparseInterpol (ppA(random_element,x), ppB(random_element,x),
     3775                                  skeleton, alpha, sparseFail, coeffMonoms,
     3776                                  Monoms);
     3777        TIMING_END_AND_PRINT (gcd_recursion,
     3778                              "time for recursive call: ");
     3779        DEBOUTLN (cerr, "G_random_element= " << G_random_element);
     3780      }
     3781      else if (fail && inextension)
     3782      {
     3783        source= CFList();
     3784        dest= CFList();
     3785        int num_vars= tmin (getNumVars(A), getNumVars(B));
     3786        num_vars--;
     3787        V_buf= alpha;
     3788        choose_extension (d, num_vars, V_buf);
     3789        bool prim_fail= false;
     3790        Variable V_buf2;
     3791        CanonicalForm prim_elem, im_prim_elem;
     3792        prim_elem= primitiveElement (alpha, V_buf2, prim_fail);
     3793       
     3794        ASSERT (!prim_fail, "failure in integer factorizer");
     3795        if (prim_fail)
     3796          ; //ERROR
     3797        else
     3798          im_prim_elem= mapPrimElem (prim_elem, alpha, V_buf);
     3799
     3800        DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (alpha));
     3801        DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf2));
     3802
     3803        inextensionextension= true;
     3804        for (CFListIterator i= l; i.hasItem(); i++)
     3805          i.getItem()= mapUp (i.getItem(), alpha, V_buf, prim_elem,
     3806                              im_prim_elem, source, dest);
     3807        m= mapUp (m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
     3808        G_m= mapUp (G_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
     3809        newtonPoly= mapUp (newtonPoly, alpha, V_buf, prim_elem, im_prim_elem,
     3810                            source, dest);
     3811        ppA= mapUp (ppA, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
     3812        ppB= mapUp (ppB, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
     3813        gcdlcAlcB= mapUp (gcdlcAlcB, alpha, V_buf, prim_elem, im_prim_elem,
     3814                          source, dest);
     3815        fail= false;
     3816        random_element= randomElement (m, V_buf, l, fail );
     3817        DEBOUTLN (cerr, "fail= " << fail);
     3818        CFList list;
     3819        TIMING_START (gcd_recursion);
     3820        if (LC (skeleton).inCoeffDomain())
     3821          G_random_element=
     3822          monicSparseInterpol (ppA (random_element, x), ppB (random_element, x),
     3823                               skeleton, V_buf, sparseFail, coeffMonoms,
     3824                               Monoms);
     3825        else
     3826          G_random_element=
     3827          nonMonicSparseInterpol (ppA(random_element,x), ppB(random_element,x),
     3828                                  skeleton, V_buf, sparseFail,
     3829                                  coeffMonoms, Monoms);
     3830        TIMING_END_AND_PRINT (gcd_recursion,
     3831                              "time for recursive call: ");
     3832        DEBOUTLN (cerr, "G_random_element= " << G_random_element);
     3833        alpha= V_buf;
     3834      }
     3835
     3836      if (sparseFail)
     3837        break;
     3838
     3839      d0= totaldegree (G_random_element, Variable(2),
     3840                      Variable (G_random_element.level()));
     3841      if (d0 == 0) 
     3842      {
     3843        if (substitute > 1)
     3844          return N(reverseSubst (gcdcAcB, substitute));
     3845        else
     3846          return N(gcdcAcB);
     3847      } 
     3848      if (d0 >  d)
     3849      {
     3850        if (!find (l, random_element))
     3851          l.append (random_element);
     3852        continue;
     3853      }
     3854
     3855      G_random_element=
     3856      (gcdlcAlcB(random_element, x)/uni_lcoeff (G_random_element))
     3857      * G_random_element;
     3858
     3859      d0= totaldegree (G_random_element, Variable(2),
     3860                      Variable(G_random_element.level()));
     3861      if (d0 <  d)
     3862      {
     3863        m= gcdlcAlcB;
     3864        newtonPoly= 1;
     3865        G_m= 0;
     3866        d= d0;
     3867      }
     3868
     3869      TIMING_START (newton_interpolation);
     3870      H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
     3871      TIMING_END_AND_PRINT (newton_interpolation,
     3872                            "time for newton interpolation: ");
     3873
     3874      //termination test
     3875      if (uni_lcoeff (H) == gcdlcAlcB)
     3876      {
     3877        cH= uni_content (H);
     3878        ppH= H/cH;
     3879        ppH /= Lc (ppH);
     3880        DEBOUTLN (cerr, "ppH= " << ppH);
     3881        if (fdivides (ppH, A) && fdivides (ppH, B))
     3882        {
     3883          if (substitute > 1)
     3884          {
     3885            ppH= reverseSubst (ppH, substitute);
     3886            gcdcAcB= reverseSubst (gcdcAcB, substitute);
     3887          }
     3888          return N(gcdcAcB*ppH);
     3889        }
     3890      }
     3891
     3892      G_m= H;
     3893      newtonPoly= newtonPoly*(x - random_element);
     3894      m= m*(x - random_element);
     3895      if (!find (l, random_element))
     3896        l.append (random_element);
     3897
     3898    } while (1); //end of second do
     3899  }
     3900  } while (1); //end of first do
     3901}
     3902
     3903static inline
     3904int compress4EZGCD (const CanonicalForm& F, const CanonicalForm& G, CFMap & M,
     3905                    CFMap & N)
     3906{
     3907  int n= tmax (F.level(), G.level());
     3908  int * degsf= new int [n + 1];
     3909  int * degsg= new int [n + 1];
     3910
     3911  for (int i = 0; i <= n; i++)
     3912    degsf[i]= degsg[i]= 0;
     3913   
     3914  degsf= degrees (F, degsf);
     3915  degsg= degrees (G, degsg);
     3916 
     3917  int both_non_zero= 0;
     3918  int f_zero= 0;
     3919  int g_zero= 0;
     3920  int both_zero= 0;
     3921
     3922  for (int i= 1; i <= n; i++)
     3923  {
     3924    if (degsf[i] != 0 && degsg[i] != 0)
     3925    {
     3926      both_non_zero++;
     3927      continue;
     3928    }
     3929    if (degsf[i] == 0 && degsg[i] != 0 && i <= G.level())
     3930    {
     3931      f_zero++;
     3932      continue;
     3933    }
     3934    if (degsg[i] == 0 && degsf[i] && i <= F.level())
     3935    {
     3936      g_zero++;
     3937      continue;
     3938    }
     3939  }
     3940
     3941  if (both_non_zero == 0) return 0;
     3942
     3943  // map Variables which do not occur in both polynomials to higher levels
     3944  int k= 1;
     3945  int l= 1;
     3946  for (int i= 1; i <= n; i++)
     3947  {
     3948    if (degsf[i] != 0 && degsg[i] == 0 && i <= F.level())
     3949    {
     3950      if (k + both_non_zero != i)
     3951      {
     3952        M.newpair (Variable (i), Variable (k + both_non_zero));
     3953        N.newpair (Variable (k + both_non_zero), Variable (i));
     3954      }
     3955      k++;
     3956    }
     3957    if (degsf[i] == 0 && degsg[i] != 0 && i <= G.level())
     3958    {
     3959      if (l + g_zero + both_non_zero != i)
     3960      {
     3961        M.newpair (Variable (i), Variable (l + g_zero + both_non_zero));
     3962        N.newpair (Variable (l + g_zero + both_non_zero), Variable (i));
     3963      }
     3964      l++;
     3965    }
     3966  }
     3967
     3968  // sort Variables x_{i} in decreasing order of
     3969  // min(deg_{x_{i}}(f),deg_{x_{i}}(g))
     3970  int m= tmin (F.level(), G.level());
     3971  int max_min_deg;
     3972  k= both_non_zero;
     3973  l= 0;
     3974  int i= 1;
     3975  while (k > 0)
     3976  {
     3977    max_min_deg= tmin (degsf[i], degsg[i]);
     3978    while (max_min_deg == 0)
     3979    {
     3980      i++;
     3981      max_min_deg= tmin (degsf[i], degsg[i]);
     3982    }
     3983    for (int j= i + 1; j <=  m; j++)
     3984    {
     3985      if (tmin (degsf[j],degsg[j]) < max_min_deg)
     3986      {
     3987        max_min_deg= tmin (degsf[j], degsg[j]);
     3988        l= j;
     3989      }
     3990    }
     3991    if (l != 0)
     3992    {
     3993      if (l != k)
     3994      {
     3995        M.newpair (Variable (l), Variable(k));
     3996        N.newpair (Variable (k), Variable(l));
     3997        degsf[l]= 0;
     3998        degsg[l]= 0;
     3999        l= 0;
     4000      }
     4001      else
     4002      {
     4003        degsf[l]= 0;
     4004        degsg[l]= 0;
     4005        l= 0;
     4006      }
     4007    } 
     4008    else if (l == 0)
     4009    {
     4010      if (i != k)
     4011      {
     4012        M.newpair (Variable (i), Variable (k));
     4013        N.newpair (Variable (k), Variable (i));
     4014        degsf[i]= 0;
     4015        degsg[i]= 0;
     4016      }
     4017      else
     4018      {
     4019        degsf[i]= 0;
     4020        degsg[i]= 0;
     4021      }
     4022      i++;
     4023    }
     4024    k--;
     4025  }
     4026
     4027  delete [] degsf;
     4028  delete [] degsg;
     4029
     4030  return both_non_zero;
     4031}
     4032
     4033static
     4034CanonicalForm myShift2Zero (const CanonicalForm& F, CFList& Feval,
     4035                            const CFList& evaluation)
     4036{
     4037  CanonicalForm A= F;
     4038  int k= 2;
     4039  for (CFListIterator i= evaluation; i.hasItem(); i++, k++)
     4040    A= A (Variable (k) + i.getItem(), k);
     4041
     4042  CanonicalForm buf= A;
     4043  Feval= CFList();
     4044  Feval.append (buf);
     4045  for (k= evaluation.length() + 1; k > 2; k--)
     4046  {
     4047    buf= mod (buf, Variable (k));
     4048    Feval.insert (buf);
     4049  }
     4050  return A;
     4051}
     4052
     4053static CanonicalForm myReverseShift (const CanonicalForm& F,
     4054                                     const CFList& evaluation)
     4055{
     4056  int l= evaluation.length() + 1;
     4057  CanonicalForm result= F;
     4058  CFListIterator j= evaluation;
     4059  for (int i= 2; i < l + 1; i++, j++)
     4060  {
     4061    if (F.level() < i)
     4062      continue;
     4063    result= result (Variable (i) - j.getItem(), i);
     4064  }
     4065  return result;
     4066}
     4067
     4068static
     4069bool Hensel_P (const CanonicalForm & U, CFArray & G, const Evaluation & A,
     4070                const Variable & x, const CFArray& LCs )
     4071{
     4072  CFList factors;
     4073  factors.append (G[1]);
     4074  factors.append (G[2]);
     4075  CFList evaluation;
     4076  for (int i= A.min(); i <= A.max(); i++)
     4077    evaluation.append (A [i]);
     4078  CFList UEval;
     4079  CanonicalForm shiftedU= myShift2Zero (U, UEval, evaluation);
     4080  CFArray shiftedLCs= CFArray (2);
     4081  CFList shiftedLCsEval1, shiftedLCsEval2;
     4082  shiftedLCs[0]= myShift2Zero (LCs[1], shiftedLCsEval1, evaluation);
     4083  shiftedLCs[1]= myShift2Zero (LCs[2], shiftedLCsEval2, evaluation);
     4084  CanonicalForm LcUEval= LC (UEval.getFirst(), x);
     4085  factors.insert (1);
     4086  int liftBound= degree (UEval.getFirst(), 2) + 1;
     4087  CFArray Pi;
     4088  CFMatrix M= CFMatrix (liftBound, factors.length() - 1);
     4089  CFList diophant;
     4090  CFArray lcs= CFArray (2);
     4091  lcs [0]= shiftedLCsEval1.getFirst();
     4092  lcs [1]= shiftedLCsEval2.getFirst();
     4093  henselLift122 (UEval.getFirst(), factors, liftBound, Pi, diophant, M,
     4094                 lcs, false);
     4095
     4096  bool noChange= true;
     4097  for (CFListIterator i= factors; i.hasItem(); i++)
     4098  {
     4099    if (degree (i.getItem(), 2) != 0)
     4100      noChange= false;
     4101  }
     4102  if (noChange)
     4103    return !noChange;
     4104  int * liftBounds;
     4105  if (U.level() > 2)
     4106  {
     4107    liftBounds= new int [U.level() - 1];
     4108    liftBounds[0]= liftBound;
     4109    for (int i= 1; i < U.level() - 1; i++)
     4110      liftBounds[i]= degree (shiftedU, Variable (i + 2)) + 1;     
     4111    factors= henselLift2 (UEval, factors, liftBounds, U.level() - 1, false,
     4112                          shiftedLCsEval1, shiftedLCsEval2, Pi, diophant);
     4113  }
     4114  G[1]= factors.getFirst();
     4115  G[2]= factors.getLast();
     4116  G[1]= myReverseShift (G[1], evaluation);
     4117  G[2]= myReverseShift (G[2], evaluation);
     4118  return true;
     4119}
     4120
     4121static
     4122bool findeval_P (const CanonicalForm & F, const CanonicalForm & G,
     4123                 CanonicalForm & Fb, CanonicalForm & Gb, CanonicalForm & Db,
     4124                 FFREvaluation & b, int delta, int degF, int degG, int maxeval,
     4125                 int & count )
     4126{
     4127  if( count == 0 && delta != 0)
     4128  {
     4129    if( count++ > maxeval )
     4130      return false;
     4131  }
     4132  if (count > 0)
     4133  {
     4134    b.nextpoint();
     4135    if (count++ > maxeval)
     4136      return false;
     4137  }
     4138  while( true )
     4139  {
     4140    Fb = b( F );
     4141    if( degree( Fb, 1 ) == degF )
     4142    {
     4143      Gb = b( G );
     4144      if( degree( Gb, 1 ) == degG )
     4145      {
     4146        Db = gcd( Fb, Gb );
     4147        if( delta > 0 )
     4148        {
     4149          if( degree( Db, 1 ) <= delta )
     4150            return true;
     4151        }
     4152        else
     4153          return true;
     4154      }
     4155    }
     4156    b.nextpoint();
     4157    if( count++ > maxeval )
     4158      return false;
     4159  }
     4160}
     4161
     4162/// Extended Zassenhaus GCD for finite fields
     4163CanonicalForm EZGCD_P( const CanonicalForm & FF, const CanonicalForm & GG )
     4164{
     4165  if (FF.isZero() && degree(GG) > 0) return GG/Lc(GG);
     4166  else if (GG.isZero() && degree (FF) > 0) return FF/Lc(FF);
     4167  else if (FF.isZero() && GG.isZero()) return FF.genOne();
     4168  if (FF.inBaseDomain() || GG.inBaseDomain()) return FF.genOne();
     4169  if (FF.isUnivariate() && fdivides(FF, GG)) return FF/Lc(FF);
     4170  if (GG.isUnivariate() && fdivides(GG, FF)) return GG/Lc(GG);
     4171  if (FF == GG) return FF/Lc(FF);
     4172 
     4173  CanonicalForm F, G, f, g, d, Fb, Gb, Db, Fbt, Gbt, Dbt, B0, B, D0, lcF, lcG,
     4174                lcD;
     4175  CFArray DD( 1, 2 ), lcDD( 1, 2 );
     4176  int degF, degG, delta, count;
     4177  int maxeval;
     4178  // bound on the number of eval. to use
     4179  maxeval = (int) ceil (getCharacteristic()/log (getCharacteristic()))*2;
     4180  count = 0; // number of eval. used
     4181  FFREvaluation b, bt;
     4182  bool gcdfound = false;
     4183  Variable x = Variable(1);
     4184
     4185  F= FF;
     4186  G= GG;
     4187 
     4188  CFMap M,N;
     4189  int best_level= compress4EZGCD (F, G, M, N);
     4190 
     4191  if (best_level == 0) return G.genOne();
     4192 
     4193  F= M (F);
     4194  G= M (G);
     4195 
     4196  f = content( F, x ); g = content( G, x ); d = gcd( f, g );
     4197  F /= f; G /= g;
     4198 
     4199  if( F.isUnivariate() && G.isUnivariate() )
     4200  {
     4201    if( F.mvar() == G.mvar() )
     4202      d *= gcd( F, G );
     4203    return N (d);
     4204  }
     4205
     4206  int maxNumVars= tmax (getNumVars (F), getNumVars (G));
     4207  Variable a,bv;
     4208  if (size (F)/maxNumVars > 500 && size (G)/maxNumVars > 500) //TODO find a
     4209                                                              //better bound?
     4210  {
     4211    if (hasFirstAlgVar (F, a) || hasFirstAlgVar (G, a))
     4212      return N (d*GCD_Fp_extension (F, G, a));
     4213    else if (CFFactory::gettype() == GaloisFieldDomain)
     4214      return N (d*GCD_GF (F, G));
     4215    else
     4216      return N (d*GCD_small_p (F, G));
     4217  }
     4218 
     4219  if( gcd_test_one( F, G, false ) )
     4220  {
     4221    return N (d);
     4222  }
     4223 
     4224  lcF = LC( F, x ); lcG = LC( G, x );
     4225  lcD = gcd( lcF, lcG );
     4226 
     4227  if (fdivides (lcF, lcG)) {
     4228    if (fdivides (F, G))
     4229      return FF/Lc(FF);   
     4230  }
     4231  if (fdivides (lcG, lcF))
     4232  {
     4233    if (fdivides (G, F))
     4234      return GG/Lc(GG);
     4235  }
     4236 
     4237  delta = 0;
     4238  degF = degree( F, x ); degG = degree( G, x );
     4239  if(hasFirstAlgVar(F,a))
     4240  {
     4241    if(hasFirstAlgVar(G,bv))
     4242    {
     4243      if(bv.level() > a.level())
     4244        a = bv;
     4245    }
     4246    b = FFREvaluation( 2, tmax(F.level(), G.level()), AlgExtRandomF( a ) );
     4247  }
     4248  else // F not in extension
     4249  {
     4250    if(hasFirstAlgVar(G,a))
     4251      b = FFREvaluation( 2, tmax(F.level(), G.level()), AlgExtRandomF( a ) );
     4252    else
     4253    { // both not in extension given by algebraic variable
     4254      if (CFFactory::gettype() != GaloisFieldDomain)
     4255        b = FFREvaluation( 2, tmax(F.level(), G.level()), FFRandom() );
     4256      else
     4257        b = FFREvaluation( 2, tmax(F.level(), G.level()), GFRandom() );
     4258    }
     4259  }
     4260  CanonicalForm cand;
     4261  while( !gcdfound )
     4262  {
     4263    if( !findeval_P( F, G, Fb, Gb, Db, b, delta, degF, degG, maxeval, count ))
     4264    { // too many eval. used --> try another method
     4265      if (hasFirstAlgVar (F, a) || hasFirstAlgVar (G, a))
     4266        return N (d*sparseGCDFq (F, G, a));
     4267      if (CFFactory::gettype() == GaloisFieldDomain)
     4268        return N (d*GCD_GF (F, G));
     4269      else
     4270        return N (d*sparseGCDFp (F,G));
     4271    }
     4272    delta = degree( Db );
     4273    if( delta == 0 )
     4274      return N (d);
     4275    while( true )
     4276    {
     4277      bt = b;
     4278      if( !findeval_P( F, G, Fbt, Gbt, Dbt, bt, delta, degF, degG, maxeval, count))
     4279      { // too many eval. used --> try another method
     4280        if (hasFirstAlgVar (F, a) || hasFirstAlgVar (G, a))
     4281          return N (d*sparseGCDFq (F, G, a));
     4282        if (CFFactory::gettype() == GaloisFieldDomain)
     4283          return N (d*GCD_GF (F, G));
     4284        else
     4285          return N (d*sparseGCDFp (F,G));
     4286      }
     4287      int dd = degree( Dbt );
     4288      if( dd == 0 )
     4289        return N (d);
     4290      if( dd == delta )
     4291        break;
     4292      if( dd < delta )
     4293      {
     4294        delta = dd;
     4295        b = bt;
     4296        Db = Dbt; Fb = Fbt; Gb = Gbt;
     4297      }
     4298    }
     4299    if( degF <= degG && delta == degF && fdivides( F, G ) )
     4300      return N (d*F);
     4301    if( degG < degF && delta == degG && fdivides( G, F ) )
     4302      return N (d*G);
     4303    if( delta != degF && delta != degG )
     4304    {
     4305      bool B_is_F;
     4306      CanonicalForm xxx1, xxx2;
     4307      CanonicalForm buf;
     4308      DD[1] = Fb / Db;
     4309      buf= Gb/Db;
     4310      xxx1 = gcd( DD[1], Db );
     4311      xxx2 = gcd( buf, Db );
     4312      if (((xxx1.inCoeffDomain() && xxx2.inCoeffDomain()) &&
     4313          (size (F) <= size (G)))
     4314          || (xxx1.inCoeffDomain() && !xxx2.inCoeffDomain()))
     4315      {
     4316        B = F;
     4317        DD[2] = Db;
     4318        lcDD[1] = lcF;
     4319        lcDD[2] = lcD;
     4320        B_is_F = true;
     4321      }
     4322      else if (((xxx1.inCoeffDomain() && xxx2.inCoeffDomain()) &&
     4323               (size (G) < size (F)))
     4324               || (!xxx1.inCoeffDomain() && xxx2.inCoeffDomain()))
     4325      {
     4326        DD[1] = buf;
     4327        B = G;
     4328        DD[2] = Db;
     4329        lcDD[1] = lcG;
     4330        lcDD[2] = lcD;
     4331        B_is_F = false;
     4332      }
     4333      else // special case handling
     4334      {
     4335        if (hasFirstAlgVar (F, a) || hasFirstAlgVar (G, a))
     4336          return N (d*sparseGCDFq (F, G, a));
     4337        if (CFFactory::gettype() == GaloisFieldDomain)
     4338          return N (d*GCD_GF (F, G));
     4339        else
     4340          return N (d*sparseGCDFp (F,G));
     4341      }
     4342      DD[2] = DD[2] * ( b( lcDD[2] ) / lc( DD[2] ) );
     4343      DD[1] = DD[1] * ( b( lcDD[1] ) / lc( DD[1] ) );
     4344
     4345      gcdfound= Hensel_P (B*lcD, DD, b, x, lcDD);
     4346     
     4347      if( gcdfound )
     4348      {
     4349        cand = DD[2] / content( DD[2], Variable(1) );
     4350        if( B_is_F )
     4351          gcdfound = fdivides( cand, G );
     4352        else
     4353          gcdfound = fdivides( cand, F );
     4354      }
     4355    }
     4356    delta= 0;
     4357  }
     4358  return N (d*cand);
     4359}
     4360
    10614361#endif

  • factory/cf_gcd_smallp.h

    r231b0fd r018577  
    6464
    6565CanonicalForm
    66 randomIrredpoly (int i, const Variable & x) ;
     66randomIrredpoly (int i, const Variable & x);
     67
     68CanonicalForm
     69sparseGCDFp (const CanonicalForm& F, const CanonicalForm& G,
     70             bool& topLevel, CFList& l);
     71             
     72static inline CanonicalForm sparseGCDFp (const CanonicalForm& A, const CanonicalForm& B)
     73{
     74  CFList list;
     75  bool top_level= true;
     76  return sparseGCDFp (A, B, top_level, list);
     77}
     78
     79CanonicalForm
     80sparseGCDFq (const CanonicalForm& F, const CanonicalForm& G, const Variable& alpha, CFList& l,
     81             bool& topLevel);
     82             
     83static inline CanonicalForm sparseGCDFq (const CanonicalForm& A, const CanonicalForm& B, const Variable& alpha)
     84{
     85  CFList list;
     86  bool top_level= true;
     87  return sparseGCDFq (A, B, alpha, list, top_level);
     88}
     89
     90/// Extended Zassenhaus GCD over finite fields described in "Algorithms for
     91/// Computer Algebra" by Geddes, Czapor & Labahn
     92CanonicalForm EZGCD_P (const CanonicalForm& F, const CanonicalForm& G);
     93
     94/// pseudo remainder GCD over finite fields
     95CanonicalForm
     96gcd_poly_p( const CanonicalForm & f, const CanonicalForm & g);
     97
     98CanonicalForm
     99gcd_univar_ntlp( const CanonicalForm & F, const CanonicalForm & G );
     100
     101
    67102#endif

  • factory/cf_switches.cc

    r231b0fd r018577  
    3131  On(SW_USE_NTL);
    3232  On(SW_USE_CHINREM_GCD);
     33  On(SW_USE_EZGCD_P);
    3334  //Off(SW_USE_NTL_GCD_0);
    3435  //Off(SW_USE_NTL_GCD_P);
    … …  
    3637#endif
    3738  On(SW_USE_EZGCD);
    38   //On(SW_USE_EZGCD_P); // still testing
    3939  On(SW_USE_QGCD);
    4040  //On(SW_USE_fieldGCD);

  • factory/facHensel.cc

    r231b0fd r018577  
    924924void
    925925henselLift12 (const CanonicalForm& F, CFList& factors, int l, CFArray& Pi,
    926               CFList& diophant, CFMatrix& M, bool sort)
    927 {
    928   if (sort)
    929     sortList (factors, Variable (1));
     926              CFList& diophant, CFMatrix& M)
     927{
     928  sortList (factors, Variable (1));
    930929  Pi= CFArray (factors.length() - 1);
    931930  CFListIterator j= factors;
Note: See TracChangeset for help on using the changeset viewer.