Changeset c4f4fd in git


Ignore:
Timestamp:
Sep 16, 2010, 5:34:13 PM (13 years ago)
Author:
Hans Schoenemann <hannes@…>
Branches:
(u'spielwiese', '8d54773d6c9e2f1d2593a28bc68b7eeab54ed529')
Children:
8e7db4cc8c23b1a45ca71f3217a02819e04d9c84
Parents:
7b8a4166957e36a8fb67484b68a9449846e1995d
Message:
going back to r13182 (facstd.tst fails)

git-svn-id: file:///usr/local/Singular/svn/trunk@13206 2c84dea3-7e68-4137-9b89-c4e89433aadc
Location:
factory
Files:
8 edited

Legend:

Unmodified
Added
Removed

  • factory/cf_gcd.cc

    r7b8a416 rc4f4fd  
    2929bool isPurePoly(const CanonicalForm & );
    3030static CanonicalForm gcd_univar_ntl0( const CanonicalForm &, const CanonicalForm & );
     31static CanonicalForm gcd_univar_ntlp( const CanonicalForm &, const CanonicalForm & );
    3132#endif
    3233
    … …  
    366367
    367368static CanonicalForm
     369gcd_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
     447static CanonicalForm
    368448gcd_poly_0( const CanonicalForm & f, const CanonicalForm & g )
    369449{
    … …  
    486566    else if (isOn( SW_USE_EZGCD_P ) && (!fc_and_gc_Univariate))
    487567    {
    488       fc= EZGCD_P (fc, gc);
     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      }
    489583    }
    490584    else if (isOn(SW_USE_GCD_P))
    … …  
    832926}
    833927
     928static CanonicalForm
     929gcd_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
    834961#endif
    835962

  • factory/cf_gcd_smallp.cc

    r7b8a416 rc4f4fd  
    3030#include "ftmpl_functions.h"
    3131#include "cf_random.h"
    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:
     32
     33// iinline helper functions:
    4034#include "cf_map_ext.h"
    4135
    … …  
    5246TIMING_DEFINE_PRINT(newton_interpolation);
    5347
    54 CanonicalForm
    55 gcd_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 
    89 CanonicalForm
    90 gcd_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
    16548/// compressing two polynomials F and G, M is used for compressing,
    16649/// N to reverse the compression
    16750static inline
    16851int myCompress (const CanonicalForm& F, const CanonicalForm& G, CFMap & M,
    169                 CFMap & N, const bool& topLevel)
     52                CFMap & N, bool& top_level)
    17053{
    17154  int n= tmax (F.level(), G.level());
    … …  
    18467  int both_zero= 0;
    18568
    186   if (topLevel)
     69  if (top_level)
    18770  {
    18871    for (int i= 1; i <= n; i++)
    … …  
    315198  delete [] degsg;
    316199
    317   if (topLevel)
    318     return both_non_zero;
    319   else
    320     return 1;
     200  return 1;
    321201}
    322 
    323 CanonicalForm
    324 uniCoeff (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 
    346 int
    347 substituteCheck (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
    409 void
    410 subst (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 
    434 CanonicalForm
    435 reverseSubst (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 
    449 static inline CanonicalForm
    450 uni_content (const CanonicalForm & F);
    451 
    452 CanonicalForm
    453 uni_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 }
    471202
    472203/// compute the content of F, where F is considered as an element of
    … …  
    481212  if (F.level() != 1 && F.isUnivariate())
    482213    return F.genOne();
    483   if (degree (F, 1) == 0) return F.genOne();
     214  if (degree (F,1) == 0) return F.genOne();
    484215
    485216  int l= F.level();
    … …  
    502233}
    503234
    504 CanonicalForm
    505 extractContents (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 
    529235/// compute the leading coefficient of F, where F is considered as an element
    530236/// of \f$ R[x_{1}][x_{2},\ldots ,x_{n}] \f$, order on Mon (x_{2},\ldots ,x_{n})
    531 /// is Dp.
     237/// is dp.
    532238static inline
    533239CanonicalForm uni_lcoeff (const CanonicalForm& F)
    … …  
    580286  int p= getCharacteristic ();
    581287  int d= degree (mipo);
    582   double bound= pow ((double) p, d);
     288  int bound= ipower (p, d);
    583289  do
    584290  {
    … …  
    620326  ZZ_pX NTLirredpoly;
    621327  //TODO: replace d by max_{i} (deg_x{i}(f))
    622   int i= (int) (log ((double) pow ((double) d + 1, (double) num_vars))/log ((double) p));
     328  int i= (int) (log ((double) ipower (d + 1, num_vars))/log ((double) p));
    623329  int m= degree (getMipo (beta));
    624330  if (i <= 1)
    … …  
    630336
    631337/// GCD of F and G over \f$ F_{p}(\alpha ) \f$ ,
    632 /// l and topLevel are only used internally, output is monic
     338/// l and top_level are only used internally, output is monic
    633339/// based on Alg. 7.2. as described in "Algorithms for
    634340/// Computer Algebra" by Geddes, Czapor, Labahn
    635341CanonicalForm
    636342GCD_Fp_extension (const CanonicalForm& F, const CanonicalForm& G,
    637                   Variable & alpha, CFList& l, bool& topLevel)
     343                  Variable & alpha, CFList& l, bool& top_level)
    638344{
    639345  CanonicalForm A= F;
    … …  
    648354 
    649355  CFMap M,N;
    650   int best_level= myCompress (A, B, M, N, topLevel);
     356  int best_level= myCompress (A, B, M, N, top_level);
    651357
    652358  if (best_level == 0) return B.genOne();
    … …  
    660366  if (A.isUnivariate() && B.isUnivariate())
    661367    return N (gcd(A,B));
    662 
    663   int substitute= substituteCheck (A, B);
    664368 
    665   if (substitute > 1)
    666     subst (A, B, A, B, substitute);
    667    
    668369  CanonicalForm cA, cB;    // content of A and B
    669370  CanonicalForm ppA, ppB;    // primitive part of A and B
    670371  CanonicalForm gcdcAcB;
    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   }
     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;
    686380
    687381  CanonicalForm lcA, lcB;  // leading coefficients of A and B
    … …  
    705399  int d= totaldegree (ppA, Variable(2), Variable (ppA.level()));
    706400
    707   if (d == 0)
    708   {
    709     if (substitute > 1)
    710       return N (reverseSubst (gcdcAcB, substitute));
    711     else
    712       return N(gcdcAcB);
    713   }
     401  if (d == 0) return N(gcdcAcB);
    714402  int d0= totaldegree (ppB, Variable(2), Variable (ppB.level()));
    715403  if (d0 < d)
    716404    d= d0;
    717   if (d == 0)
    718   {
    719     if (substitute > 1)
    720       return N (reverseSubst (gcdcAcB, substitute));
    721     else
    722       return N(gcdcAcB);
    723   }
     405  if (d == 0) return N(gcdcAcB);
    724406
    725407  CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH;
    … …  
    731413  H= 0;
    732414  bool fail= false;
    733   topLevel= false;
     415  top_level= false;
    734416  bool inextension= false;
    735417  Variable V_buf= alpha;
    … …  
    745427      int num_vars= tmin (getNumVars(A), getNumVars(B));
    746428      num_vars--;
    747      
     429
    748430      choose_extension (d, num_vars, V_buf);
    749431      bool prim_fail= false;
    750432      Variable V_buf2;
    751433      prim_elem= primitiveElement (alpha, V_buf2, prim_fail);
     434
    752435      ASSERT (!prim_fail, "failure in integer factorizer");
    753436      if (prim_fail)
    … …  
    755438      else
    756439        im_prim_elem= mapPrimElem (prim_elem, alpha, V_buf);
    757      
     440
    758441      DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (alpha));
    759442      DEBOUTLN (cerr, "getMipo (V_buf2)= " << getMipo (V_buf2));
    … …  
    778461      G_random_element=
    779462      GCD_Fp_extension (ppA (random_element, x), ppB (random_element, x), V_buf,
    780                         list, topLevel);
     463                        list, top_level);
    781464      TIMING_END_AND_PRINT (gcd_recursion,
    782465                            "time for recursive call: ");
    … …  
    789472      G_random_element=
    790473      GCD_Fp_extension (ppA(random_element, x), ppB(random_element, x), V_buf,
    791                         list, topLevel);
     474                        list, top_level);
    792475      TIMING_END_AND_PRINT (gcd_recursion,
    793476                            "time for recursive call: ");
    … …  
    797480    d0= totaldegree (G_random_element, Variable(2),
    798481                     Variable (G_random_element.level()));
    799     if (d0 == 0)
    800     {
    801       if (substitute > 1)
    802         return N (reverseSubst (gcdcAcB, substitute));
    803       else
    804         return N(gcdcAcB);
    805     }
     482    if (d0 == 0) return N(gcdcAcB);
    806483    if (d0 >  d)
    807484    {
    … …  
    844521        ppH /= Lc(ppH);
    845522        DEBOUTLN (cerr, "ppH after mapDown= " << ppH);
    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           }
     523        if (fdivides (ppH, A) && fdivides (ppH, B))         
    853524          return N(gcdcAcB*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       }
     525      }
     526      else if (fdivides (ppH, A) && fdivides (ppH, B))
     527        return N(gcdcAcB*(ppH/Lc(ppH)));
    865528    }
    866529
    … …  
    886549  int p= getCharacteristic();
    887550  int d= getGFDegree();
    888   double bound= pow ((double) p, (double) d);
     551  int bound= ipower (p, d);
    889552  do
    890553  {
    … …  
    917580/// the size of the base field is bounded by 2^16, output is monic
    918581CanonicalForm GCD_GF (const CanonicalForm& F, const CanonicalForm& G,
    919         CFList& l, bool& topLevel)
     582        CFList& l, bool& top_level)
    920583{
    921584  CanonicalForm A= F;
    … …  
    930593 
    931594  CFMap M,N;
    932   int best_level= myCompress (A, B, M, N, topLevel);
     595  int best_level= myCompress (A, B, M, N, top_level);
    933596
    934597  if (best_level == 0) return B.genOne();
    … …  
    989652  H= 0;
    990653  bool fail= false;
    991   topLevel= false;
     654  top_level= false;
    992655  bool inextension= false;
    993656  int p;
    … …  
    1004667      num_vars--;
    1005668      p= getCharacteristic();
    1006       expon= (int) floor ((log ((double) pow ((double) d + 1, num_vars))/log ((double) p)));
     669      expon= (int) floor ((log ((double) ipower (d + 1, num_vars))/log ((double) p)));
    1007670      if (expon < 2)
    1008671        expon= 2;
    1009672      kk= getGFDegree();
    1010       if (pow ( (double) p, kk*expon) < (1 << 16))
     673      if (ipower (p, kk*expon) < (1 << 16))
    1011674        setCharacteristic (p, kk*(int)expon, 'b');
    1012675      else
    … …  
    1031694      TIMING_START (gcd_recursion);
    1032695      G_random_element= GCD_GF (ppA(random_element, x), ppB(random_element, x),
    1033                                 list, topLevel);
     696                                list, top_level);
    1034697      TIMING_END_AND_PRINT (gcd_recursion,
    1035698                            "time for recursive call: ");
    … …  
    1041704      TIMING_START (gcd_recursion);
    1042705      G_random_element= GCD_GF (ppA(random_element, x), ppB(random_element, x),
    1043                                 list, topLevel);
     706                                list, top_level);
    1044707      TIMING_END_AND_PRINT (gcd_recursion,
    1045708                            "time for recursive call: ");
    … …  
    1166829
    1167830CanonicalForm GCD_small_p (const CanonicalForm& F, const CanonicalForm&  G,
    1168                            bool& topLevel, CFList& l)
     831                           bool& top_level, CFList& l)
    1169832{
    1170833  CanonicalForm A= F;
    … …  
    1179842
    1180843  CFMap M,N;
    1181   int best_level= myCompress (A, B, M, N, topLevel);
     844  int best_level= myCompress (A, B, M, N, top_level);
    1182845
    1183846  if (best_level == 0) return B.genOne();
    … …  
    1185848  A= M(A);
    1186849  B= M(B);
     850
     851  Variable x= Variable (1);
    1187852
    1188853  //univariate case
    1189854  if (A.isUnivariate() && B.isUnivariate())
    1190855    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);
    1198856
    1199857  CanonicalForm cA, cB;    // content of A and B
    1200858  CanonicalForm ppA, ppB;    // primitive part of A and B
    1201859  CanonicalForm gcdcAcB;
    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   }
     860  cA = uni_content (A);
     861  cB = uni_content (B);
     862  gcdcAcB= gcd (cA, cB);
     863  ppA= A/cA;
     864  ppB= B/cB;
    1217865
    1218866  CanonicalForm lcA, lcB;  // leading coefficients of A and B
    … …  
    1231879      return G/Lc(G);
    1232880  }
    1233  
     881
    1234882  gcdlcAlcB= gcd (lcA, lcB);
    1235883 
    … …  
    1237885  int d0;
    1238886
    1239   if (d == 0)
    1240   {
    1241     if (substitute > 1)
    1242       return N (reverseSubst (gcdcAcB, substitute));
    1243     else
    1244       return N(gcdcAcB);
    1245   }
     887  if (d == 0) return N(gcdcAcB);
    1246888  d0= totaldegree (ppB, Variable (2), Variable (ppB.level()));
    1247889
    … …  
    1249891    d= d0;
    1250892
    1251   if (d == 0)
    1252   {
    1253     if (substitute > 1)
    1254       return N (reverseSubst (gcdcAcB, substitute));
    1255     else
    1256       return N(gcdcAcB);
    1257   }
     893  if (d == 0) return N(gcdcAcB);
    1258894
    1259895  CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH;
    … …  
    1263899  G_m= 0;
    1264900  Variable alpha, V_buf;
     901  int p, k;
    1265902  bool fail= false;
    1266903  bool inextension= false;
    1267904  bool inextensionextension= false;
    1268   topLevel= false;
     905  top_level= false;
     906  CanonicalForm CF_buf;
    1269907  CFList source, dest;
     908  CanonicalForm gcdcheck;
    1270909  do
    1271910  {
    … …  
    1280919      TIMING_START (gcd_recursion);
    1281920      G_random_element=
    1282       GCD_small_p (ppA (random_element,x), ppB (random_element,x), topLevel,
     921      GCD_small_p (ppA (random_element,x), ppB (random_element,x), top_level,
    1283922      list);
    1284923      TIMING_END_AND_PRINT (gcd_recursion,
    … …  
    1292931      G_random_element=
    1293932      GCD_Fp_extension (ppA (random_element, x), ppB (random_element, x), alpha,
    1294                         list, topLevel);
     933                        list, top_level);
    1295934      TIMING_END_AND_PRINT (gcd_recursion,
    1296935                            "time for recursive call: ");
    … …  
    1303942      CFList list;
    1304943      CanonicalForm mipo;
    1305       int deg= 2;
     944      int deg= 3;
    1306945      do {
    1307         l= CFList();
    1308946        mipo= randomIrredpoly (deg, x);
    1309947        alpha= rootOf (mipo);
    … …  
    1312950        random_element= randomElement (m, alpha, l, fail);
    1313951        deg++;
    1314       } while (fail);
     952      } while (fail); 
    1315953      list= CFList();
    1316954      TIMING_START (gcd_recursion);
    1317955      G_random_element=
    1318956      GCD_Fp_extension (ppA (random_element, x), ppB (random_element, x), alpha,
    1319                         list, topLevel);
     957                        list, top_level);
    1320958      TIMING_END_AND_PRINT (gcd_recursion,
    1321959                            "time for recursive call: ");
    … …  
    1334972      CanonicalForm prim_elem, im_prim_elem;
    1335973      prim_elem= primitiveElement (alpha, V_buf2, prim_fail);
     974       
    1336975      ASSERT (!prim_fail, "failure in integer factorizer");
    1337976      if (prim_fail)
    … …  
    13621001      G_random_element=
    13631002      GCD_Fp_extension (ppA (random_element, x), ppB (random_element, x), V_buf,
    1364                         list, topLevel);
     1003                        list, top_level);
    13651004      TIMING_END_AND_PRINT (gcd_recursion,
    13661005                            "time for recursive call: ");
    … …  
    13741013    if (d0 == 0)
    13751014    {
    1376       if (substitute > 1)
    1377         return N (reverseSubst (gcdcAcB, substitute));
    1378       else
    1379         return N(gcdcAcB);
     1015      return N(gcdcAcB);
    13801016    }
    13811017    if (d0 >  d)
    … …  
    14131049      ppH /= Lc (ppH);
    14141050      DEBOUTLN (cerr, "ppH= " << ppH);
    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         }
     1051      if (fdivides (ppH, A) && fdivides (ppH, B))
    14221052        return N(gcdcAcB*ppH);
    1423       }
    14241053    }
    14251054
    … …  
    14321061}
    14331062
    1434 CFArray
    1435 solveVandermonde (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
    1488 CFArray
    1489 solveVandermonde (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 
    1544 CFArray
    1545 solveGeneralVandermonde (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 
    1600 CFArray
    1601 solveGeneralVandermonde (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
    1657 CFArray
    1658 readOffSolution (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 
    1679 CFArray
    1680 readOffSolution (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 
    1708 long
    1709 gaussianElimFp (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 
    1736 long
    1737 gaussianElimFq (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 
    1766 CFArray
    1767 solveSystemFp (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 
    1794 CFArray
    1795 solveSystemFq (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 
    1823 CFArray
    1824 getMonoms (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 
    1857 CFArray
    1858 evaluateMonom (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 
    1896 CFArray
    1897 evaluate (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 
    1913 CFList
    1914 evaluationPoints (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
    2042 evalPointsSmallChar (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
    2135 void 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 
    2144 void 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 
    2157 CanonicalForm
    2158 monicSparseInterpol (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 
    2419 CanonicalForm
    2420 nonMonicSparseInterpol (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 
    3004 CanonicalForm 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 
    3418 CanonicalForm 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 
    3903 static inline
    3904 int 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 
    4033 static
    4034 CanonicalForm 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 
    4053 static 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 
    4068 static
    4069 bool 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 
    4121 static
    4122 bool 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
    4163 CanonicalForm 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 = (getCharacteristic()/ilog2 (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 
    43611063#endif

  • factory/cf_gcd_smallp.h

    r7b8a416 rc4f4fd  
    6464
    6565CanonicalForm
    66 randomIrredpoly (int i, const Variable & x);
    67 
    68 CanonicalForm
    69 sparseGCDFp (const CanonicalForm& F, const CanonicalForm& G,
    70              bool& topLevel, CFList& l);
    71              
    72 static 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 
    79 CanonicalForm
    80 sparseGCDFq (const CanonicalForm& F, const CanonicalForm& G, const Variable& alpha, CFList& l,
    81              bool& topLevel);
    82              
    83 static 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
    92 CanonicalForm EZGCD_P (const CanonicalForm& F, const CanonicalForm& G);
    93 
    94 /// pseudo remainder GCD over finite fields
    95 CanonicalForm
    96 gcd_poly_p( const CanonicalForm & f, const CanonicalForm & g);
    97 
    98 CanonicalForm
    99 gcd_univar_ntlp( const CanonicalForm & F, const CanonicalForm & G );
    100 
    101 
     66randomIrredpoly (int i, const Variable & x) ;
    10267#endif

  • factory/cf_map_ext.cc

    r7b8a416 rc4f4fd  
    340340      return 0;
    341341  } while (1);
    342  
    343   zz_pX alpha_mipo= convertFacCF2NTLzzpX (mipo);
    344   zz_pE::init (alpha_mipo);
    345   zz_pEX NTL_beta_mipo= to_zz_pEX (NTL_mipo);
    346   zz_pE root= FindRoot (NTL_beta_mipo);
     342  zz_pE::init (NTL_mipo);
     343  zz_pEX NTL_alpha_mipo= convertFacCF2NTLzz_pEX (mipo, NTL_mipo);
     344  zz_pE root= FindRoot (NTL_alpha_mipo);
    347345  return convertNTLzzpE2CF (root, alpha);
    348346}

  • factory/cf_switches.cc

    r7b8a416 rc4f4fd  
    3131  On(SW_USE_NTL);
    3232  On(SW_USE_CHINREM_GCD);
    33   On(SW_USE_EZGCD_P);
    3433  //Off(SW_USE_NTL_GCD_0);
    3534  //Off(SW_USE_NTL_GCD_P);
    … …  
    3736#endif
    3837  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

    r7b8a416 rc4f4fd  
    2626#include "facHensel.h"
    2727#include "cf_util.h"
    28 #include "fac_util.h"
    29 #include "cf_algorithm.h"
    3028
    3129#ifdef HAVE_NTL
    3230#include <NTL/lzz_pEX.h>
    3331#include "NTLconvert.h"
    34 
    35 CanonicalForm replaceLC (const CanonicalForm& F, const CanonicalForm& c);
    3632
    3733static inline
    … …  
    937933  if (j.hasItem())
    938934    j++;
    939   for (; j.hasItem(); j++, i++)
     935  for (j; j.hasItem(); j++, i++)
    940936  {
    941937    Pi [i]= mulNTL (Pi [i - 1], j.getItem());
    … …  
    958954    k.getItem()= bufFactors[i];
    959955  factors.removeFirst();
    960   return;
    961 }
    962 
    963 void
    964 henselStep122 (const CanonicalForm& F, const CFList& factors,
    965               CFArray& bufFactors, const CFList& diophant, CFMatrix& M,
    966               CFArray& Pi, int j, const CFArray& LCs)
    967 {
    968   Variable x= F.mvar();
    969   CanonicalForm xToJ= power (x, j);
    970  
    971   CanonicalForm E;
    972   // compute the error
    973   if (degree (Pi [factors.length() - 2], x) > 0)
    974     E= F[j] - Pi [factors.length() - 2] [j];
    975   else
    976     E= F[j];
    977 
    978   CFArray buf= CFArray (diophant.length());
    979 
    980   int k= 0;
    981   CanonicalForm remainder;
    982   // actual lifting
    983   for (CFListIterator i= diophant; i.hasItem(); i++, k++)
    984   {
    985     if (degree (bufFactors[k], x) > 0)
    986       remainder= modNTL (E, bufFactors[k] [0]);
    987     else
    988       remainder= modNTL (E, bufFactors[k]);
    989     buf[k]= mulNTL (i.getItem(), remainder);
    990     if (degree (bufFactors[k], x) > 0)
    991       buf[k]= modNTL (buf[k], bufFactors[k] [0]);
    992     else
    993       buf[k]= modNTL (buf[k], bufFactors[k]);
    994   }
    995  
    996   for (k= 0; k < factors.length(); k++)
    997     bufFactors[k] += xToJ*buf[k];
    998    
    999   // update Pi [0]
    1000   int degBuf0= degree (bufFactors[0], x);
    1001   int degBuf1= degree (bufFactors[1], x);
    1002   if (degBuf0 > 0 && degBuf1 > 0)
    1003   {
    1004     M (j + 1, 1)= mulNTL (bufFactors[0] [j], bufFactors[1] [j]);
    1005     if (j + 2 <= M.rows())
    1006       M (j + 2, 1)= mulNTL (bufFactors[0] [j + 1], bufFactors[1] [j + 1]);
    1007   }
    1008   CanonicalForm uIZeroJ;
    1009   if (degBuf0 > 0 && degBuf1 > 0)
    1010     uIZeroJ= mulNTL (bufFactors[0] [0], buf [1]) + mulNTL (bufFactors[1] [0], buf[0]);
    1011   else if (degBuf0 > 0)
    1012     uIZeroJ= mulNTL (buf[0], bufFactors[1]);
    1013   else if (degBuf1 > 0)
    1014     uIZeroJ= mulNTL (bufFactors[0], buf [1]);
    1015   else
    1016     uIZeroJ= 0;
    1017   Pi [0] += xToJ*uIZeroJ;
    1018 
    1019   CFArray tmp= CFArray (factors.length() - 1);
    1020   for (k= 0; k < factors.length() - 1; k++)
    1021     tmp[k]= 0;
    1022   CFIterator one, two;
    1023   one= bufFactors [0];
    1024   two= bufFactors [1];
    1025   if (degBuf0 > 0 && degBuf1 > 0)
    1026   {
    1027     while (one.exp() > j) one++;
    1028     while (two.exp() > j) two++;
    1029     for (k= 1; k <= (int) ceil (j/2.0); k++)
    1030     { 
    1031       if (k != j - k + 1)
    1032       {
    1033         if (one.exp() == j - k + 1 && two.exp() == j - k + 1)
    1034         {
    1035           tmp[0] += mulNTL ((bufFactors[0] [k] + one.coeff()), (bufFactors[1] [k] +
    1036                      two.coeff())) - M (k + 1, 1) - M (j - k + 2, 1);
    1037           one++;
    1038           two++;
    1039         }
    1040         else if (one.exp() == j - k + 1)
    1041         {
    1042           tmp[0] += mulNTL ((bufFactors[0] [k] + one.coeff()), bufFactors[1] [k]) -
    1043                      M (k + 1, 1);
    1044           one++;
    1045         }
    1046         else if (two.exp() == j - k + 1)
    1047         {
    1048           tmp[0] += mulNTL (bufFactors[0] [k], (bufFactors[1] [k] + two.coeff())) -
    1049                     M (k + 1, 1);
    1050           two++;
    1051         }
    1052       }
    1053       else
    1054         tmp[0] += M (k + 1, 1);
    1055     }
    1056 
    1057     if (j + 2 <= M.rows())
    1058     {
    1059       if (degBuf0 >= j + 1 && degBuf1 >= j + 1)
    1060         tmp [0] += mulNTL ((bufFactors [0] [j + 1]+ bufFactors [0] [0]), (bufFactors [1] [j + 1] + bufFactors [1] [0]))
    1061                - M(1,1) - M (j + 2,1);
    1062       else if (degBuf0 >= j + 1)
    1063         tmp[0] += mulNTL (bufFactors [0] [j+1], bufFactors [1] [0]);
    1064       else if (degBuf1 >= j + 1)
    1065         tmp[0] += mulNTL (bufFactors [0] [0], bufFactors [1] [j + 1]);
    1066     }
    1067   }
    1068   Pi [0] += tmp[0]*xToJ*F.mvar();
    1069  
    1070   /*// update Pi [l]
    1071   int degPi, degBuf;
    1072   for (int l= 1; l < factors.length() - 1; l++)
    1073   {
    1074     degPi= degree (Pi [l - 1], x);
    1075     degBuf= degree (bufFactors[l + 1], x);
    1076     if (degPi > 0 && degBuf > 0)
    1077       M (j + 1, l + 1)= mulNTL (Pi [l - 1] [j], bufFactors[l + 1] [j]); 
    1078     if (j == 1)
    1079     {
    1080       if (degPi > 0 && degBuf > 0)
    1081         Pi [l] += xToJ*(mulNTL (Pi [l - 1] [0] + Pi [l - 1] [j],
    1082                   bufFactors[l + 1] [0] + buf[l + 1]) - M (j + 1, l +1) -
    1083                   M (1, l + 1));
    1084       else if (degPi > 0)
    1085         Pi [l] += xToJ*(mulNTL (Pi [l - 1] [j], bufFactors[l + 1]));
    1086       else if (degBuf > 0)
    1087         Pi [l] += xToJ*(mulNTL (Pi [l - 1], buf[l + 1]));
    1088     }
    1089     else
    1090     {
    1091       if (degPi > 0 && degBuf > 0)
    1092       {
    1093         uIZeroJ= mulNTL (uIZeroJ, bufFactors [l + 1] [0]);
    1094         uIZeroJ += mulNTL (Pi [l - 1] [0], buf [l + 1]);
    1095       }
    1096       else if (degPi > 0)
    1097         uIZeroJ= mulNTL (uIZeroJ, bufFactors [l + 1]);
    1098       else if (degBuf > 0)
    1099       {
    1100         uIZeroJ= mulNTL (uIZeroJ, bufFactors [l + 1] [0]);
    1101         uIZeroJ += mulNTL (Pi [l - 1], buf[l + 1]);
    1102       }
    1103       Pi[l] += xToJ*uIZeroJ;
    1104     }
    1105     one= bufFactors [l + 1];
    1106     two= Pi [l - 1];
    1107     if (two.exp() == j + 1)
    1108     {
    1109       if (degBuf > 0 && degPi > 0)
    1110       {
    1111           tmp[l] += mulNTL (two.coeff(), bufFactors[l + 1][0]);
    1112           two++;
    1113       }
    1114       else if (degPi > 0)
    1115       {     
    1116           tmp[l] += mulNTL (two.coeff(), bufFactors[l + 1]);
    1117           two++;
    1118       }
    1119     }
    1120     if (degBuf > 0 && degPi > 0)
    1121     {
    1122       for (k= 1; k <= (int) ceil (j/2.0); k++)
    1123       { 
    1124         if (k != j - k + 1)
    1125         {
    1126           if (one.exp() == j - k + 1 && two.exp() == j - k + 1)
    1127           {
    1128             tmp[l] += mulNTL ((bufFactors[l + 1] [k] + one.coeff()), (Pi[l - 1] [k] +
    1129                        two.coeff())) - M (k + 1, l + 1) - M (j - k + 2, l + 1);
    1130             one++;
    1131             two++;
    1132           }
    1133           else if (one.exp() == j - k + 1)
    1134           {
    1135             tmp[l] += mulNTL ((bufFactors[l + 1] [k] + one.coeff()), Pi[l - 1] [k]) -
    1136                        M (k + 1, l + 1);     
    1137             one++;
    1138           }
    1139           else if (two.exp() == j - k + 1)
    1140           {
    1141             tmp[l] += mulNTL (bufFactors[l + 1] [k], (Pi[l - 1] [k] + two.coeff())) -
    1142                       M (k + 1, l + 1);
    1143             two++;
    1144           }
    1145         }
    1146         else
    1147           tmp[l] += M (k + 1, l + 1);
    1148       } 
    1149     } 
    1150     Pi[l] += tmp[l]*xToJ*F.mvar();
    1151   }*/
    1152   return;
    1153 }
    1154 
    1155 void
    1156 henselLift122 (const CanonicalForm& F, CFList& factors, int l, CFArray& Pi,
    1157               CFList& diophant, CFMatrix& M, const CFArray& LCs, bool sort)
    1158 {
    1159   if (sort)
    1160     sortList (factors, Variable (1));
    1161   Pi= CFArray (factors.length() - 2);
    1162   CFList bufFactors2= factors;
    1163   bufFactors2.removeFirst();
    1164   diophant.removeFirst();
    1165   CFListIterator iter= diophant;
    1166   CanonicalForm test= 0, test2= 0;
    1167   CanonicalForm s,t;
    1168   extgcd (bufFactors2.getFirst(), bufFactors2.getLast(), s, t);
    1169   diophant= CFList();
    1170   diophant.append (t);
    1171   diophant.append (s);
    1172   DEBOUTLN (cerr, "diophant= " << diophant);
    1173 
    1174   CFArray bufFactors= CFArray (bufFactors2.length());
    1175   int i= 0;
    1176   for (CFListIterator k= bufFactors2; k.hasItem(); i++, k++)
    1177     bufFactors[i]= replaceLc (k.getItem(), LCs [i]);
    1178 
    1179   Variable x= F.mvar();
    1180   if (degree (bufFactors[0], x) > 0 && degree (bufFactors [1], x) > 0)
    1181   {
    1182     M (1, 1)= mulNTL (bufFactors [0] [0], bufFactors[1] [0]);
    1183     Pi [0]= M (1, 1) + (mulNTL (bufFactors [0] [1], bufFactors[1] [0]) +
    1184                         mulNTL (bufFactors [0] [0], bufFactors [1] [1]))*x;
    1185   }
    1186   else if (degree (bufFactors[0], x) > 0)
    1187   {
    1188     M (1, 1)= mulNTL (bufFactors [0] [0], bufFactors[1]);
    1189     Pi [0]= M (1, 1) +
    1190             mulNTL (bufFactors [0] [1], bufFactors[1])*x;
    1191   }
    1192   else if (degree (bufFactors[1], x) > 0)
    1193   {
    1194     M (1, 1)= mulNTL (bufFactors [0], bufFactors[1] [0]);
    1195     Pi [0]= M (1, 1) +
    1196             mulNTL (bufFactors [0], bufFactors[1] [1])*x;
    1197   }
    1198   else
    1199   {
    1200     M (1, 1)= mulNTL (bufFactors [0], bufFactors[1]);
    1201     Pi [0]= M (1, 1);
    1202   }         
    1203 
    1204   for (i= 1; i < l; i++)
    1205     henselStep122 (F, bufFactors2, bufFactors, diophant, M, Pi, i, LCs);
    1206 
    1207   factors= CFList();
    1208   for (i= 0; i < bufFactors.size(); i++)
    1209     factors.append (bufFactors[i]);
    1210956  return;
    1211957}
    … …  
    13871133  for (int i= 1; i < d; i++)
    13881134  {
    1389     CanonicalForm kk= 0;
    13901135    if (degree (e, y) > 0)
    13911136      coeffE= e[i];
    … …  
    14111156      }
    14121157    }
    1413     if (e.isZero())
    1414       break;
    1415   }
    1416 
    1417 #ifdef DEBUGOUTPUT
    1418     CanonicalForm test= 0;
    1419     j= p;
    1420     for (CFListIterator i= result; i.hasItem(); i++, j++)
    1421       test += mod (i.getItem()*j.getItem(), power (y, d));
    1422     DEBOUTLN (cerr, "test= " << test);
    1423 #endif
    1424   return result;
    1425 }
    1426 
    1427 /// solve \f$ E=sum_{i= 1}^{r}{\sigma_{i}prod_{j=1, j\neq i}^{r}{f_{i}}}\f$
    1428 /// mod M, products contains \f$ prod_{j=1, j\neq i}^{r}{f_{i}}} \f$
    1429 static inline
    1430 CFList
    1431 diophantine (const CFList& recResult, const CFList& factors,
    1432              const CFList& products, const CFList& M, const CanonicalForm& E)
    1433 {
    1434   if (M.isEmpty())
    1435   {
    1436     CFList result;
    1437     CFListIterator j= factors;
    1438     CanonicalForm buf;
    1439     for (CFListIterator i= recResult; i.hasItem(); i++, j++)
    1440     {
    1441       ASSERT (E.isUnivariate() || E.inCoeffDomain(),
    1442               "constant or univariate poly expected");
    1443       ASSERT (i.getItem().isUnivariate() || i.getItem().inCoeffDomain(),
    1444               "constant or univariate poly expected");
    1445       ASSERT (j.getItem().isUnivariate() || j.getItem().inCoeffDomain(),
    1446               "constant or univariate poly expected");
    1447       buf= mulNTL (E, i.getItem());
    1448       result.append (modNTL (buf, j.getItem()));
    1449     }
    1450     return result;
    1451   }
    1452   Variable y= M.getLast().mvar();
    1453   CFList bufFactors= factors;
    1454   for (CFListIterator i= bufFactors; i.hasItem(); i++)
    1455     i.getItem()= mod (i.getItem(), y);
    1456   CFList bufProducts= products;
    1457   for (CFListIterator i= bufProducts; i.hasItem(); i++)
    1458     i.getItem()= mod (i.getItem(), y);
    1459   CFList buf= M;
    1460   buf.removeLast();
    1461   CanonicalForm bufE= mod (E, y);
    1462   CFList recDiophantine= diophantine (recResult, bufFactors, bufProducts, buf,
    1463                                       bufE);
    1464 
    1465   CanonicalForm e= E;
    1466   CFListIterator j= products;
    1467   for (CFListIterator i= recDiophantine; i.hasItem(); i++, j++)
    1468     e -= mulMod (i.getItem(), j.getItem(), M);
    1469 
    1470   CFList result= recDiophantine;
    1471   int d= degree (M.getLast());
    1472   CanonicalForm coeffE;
    1473   for (int i= 1; i < d; i++)
    1474   {
    1475     if (degree (e, y) > 0)
    1476       coeffE= e[i];
    1477     else
    1478       coeffE= 0;
    1479     if (!coeffE.isZero())
    1480     {
    1481       CFListIterator k= result;
    1482       recDiophantine= diophantine (recResult, bufFactors, bufProducts, buf,
    1483                                    coeffE);
    1484       CFListIterator l= products;
    1485       for (j= recDiophantine; j.hasItem(); j++, k++, l++)
    1486       {
    1487         k.getItem() += j.getItem()*power (y, i);
    1488         e -= mulMod (j.getItem()*power (y, i), l.getItem(), M);
    1489       }
    1490     }
    1491     if (e.isZero())
    1492       break;
    1493   }
    1494   return result;
    1495 }
    1496 
    1497 // non monic case
    1498 static inline
    1499 CFList
    1500 multiRecDiophantine2 (const CanonicalForm& F, const CFList& factors,
    1501                      const CFList& recResult, const CFList& M, const int d)
    1502 {
    1503   Variable y= F.mvar();
    1504   CFList result;
    1505   CFListIterator i;
    1506   CanonicalForm e= 1;
    1507   CFListIterator j= factors;
    1508   CFList p;
    1509   CFArray bufFactors= CFArray (factors.length());
    1510   CanonicalForm yToD= power (y, d);
    1511   int k= 0;
    1512   for (CFListIterator i= factors; i.hasItem(); i++, k++)
    1513     bufFactors [k]= i.getItem();
    1514   CanonicalForm b;
    1515   CFList buf= M;
    1516   buf.removeLast();
    1517   buf.append (yToD);
    1518   for (k= 0; k < factors.length(); k++) //TODO compute b's faster
    1519   {
    1520     b= 1;
    1521     if (fdivides (bufFactors[k], F))
    1522       b= F/bufFactors[k];
    1523     else
    1524     {
    1525       for (int l= 0; l < factors.length(); l++)
    1526       {
    1527         if (l == k)
    1528           continue;
    1529         else
    1530         {
    1531           b= mulMod (b, bufFactors[l], buf);
    1532         }
    1533       }
    1534     }
    1535     p.append (b);
    1536   } 
    1537   j= p;
    1538   for (CFListIterator i= recResult; i.hasItem(); i++, j++)
    1539     e -= mulMod (i.getItem(), j.getItem(), M);
    1540   if (e.isZero())
    1541     return recResult;
    1542   CanonicalForm coeffE;
    1543   result= recResult;
    1544   CanonicalForm g;
    1545   buf.removeLast();
    1546   for (int i= 1; i < d; i++)
    1547   {
    1548     if (degree (e, y) > 0)
    1549       coeffE= e[i];
    1550     else
    1551       coeffE= 0;
    1552     if (!coeffE.isZero())
    1553     {
    1554       CFListIterator k= result;
    1555       CFListIterator l= p;
    1556       j= recResult;
    1557       int ii= 0;
    1558       CanonicalForm dummy;
    1559       CFList recDiophantine;
    1560       CFList buf2, buf3;
    1561       buf2= factors;
    1562       buf3= p;
    1563       for (CFListIterator iii= buf2; iii.hasItem(); iii++)
    1564         iii.getItem()= mod (iii.getItem(), y);
    1565       for (CFListIterator iii= buf3; iii.hasItem(); iii++)
    1566         iii.getItem()= mod (iii.getItem(), y);
    1567       recDiophantine= diophantine (recResult, buf2, buf3, buf, coeffE);
    1568       CFListIterator iter= recDiophantine;
    1569       for (; j.hasItem(); j++, k++, l++, ii++, iter++)
    1570       {
    1571         k.getItem() += iter.getItem()*power (y, i);
    1572         e -= mulMod (iter.getItem()*power (y, i), l.getItem(), M);
    1573       }
    1574     }
     1158 
    15751159    if (e.isZero())
    15761160      break;
    … …  
    18261410}
    18271411
    1828 // so far only tested for two! factor Hensel lifting
    1829 static inline
    1830 void
    1831 henselStep2 (const CanonicalForm& F, const CFList& factors, CFArray& bufFactors,
    1832             const CFList& diophant, CFMatrix& M, CFArray& Pi,
    1833             const CFList& products, int j, const CFList& MOD)
    1834 {
    1835   CanonicalForm E;
    1836   CanonicalForm xToJ= power (F.mvar(), j);
    1837   Variable x= F.mvar();
    1838 
    1839   // compute the error
    1840 #ifdef DEBUGOUTPUT
    1841     CanonicalForm test= 1;
    1842     for (int i= 0; i < factors.length(); i++)
    1843     {
    1844       if (i == 0)
    1845         test *= mod (bufFactors [i], power (x, j));
    1846       else
    1847         test *= bufFactors[i];
    1848     }
    1849     test= mod (test, power (x, j));
    1850     test= mod (test, MOD);
    1851     CanonicalForm test2= mod (F, power (x, j - 1)) - mod (test, power (x, j-1));
    1852     DEBOUTLN (cerr, "test= " << test2);
    1853 #endif
    1854 
    1855   if (degree (Pi [factors.length() - 2], x) > 0)
    1856     E= F[j] - Pi [factors.length() - 2] [j];
    1857   else
    1858     E= F[j];
    1859 
    1860   CFArray buf= CFArray (diophant.length());
    1861 
    1862   // actual lifting
    1863   CFList diophantine2= diophantine (diophant, factors, products, MOD, E);
    1864 
    1865   int k= 0;
    1866   for (CFListIterator i= diophantine2; k < factors.length(); k++, i++)
    1867   {
    1868     buf[k]= i.getItem();
    1869     bufFactors[k] += xToJ*i.getItem();
    1870   }
    1871 
    1872   // update Pi [0]
    1873   int degBuf0= degree (bufFactors[0], x);
    1874   int degBuf1= degree (bufFactors[1], x);
    1875   if (degBuf0 > 0 && degBuf1 > 0)
    1876   {
    1877     M (j + 1, 1)= mulMod (bufFactors[0] [j], bufFactors[1] [j], MOD);
    1878     if (j + 2 <= M.rows())
    1879       M (j + 2, 1)= mulMod (bufFactors[0] [j + 1], bufFactors[1] [j + 1], MOD);
    1880   }
    1881   CanonicalForm uIZeroJ;
    1882 
    1883     if (degBuf0 > 0 && degBuf1 > 0)
    1884       uIZeroJ= mulMod (bufFactors[0] [0], buf[1], MOD) +
    1885                mulMod (bufFactors[1] [0], buf[0], MOD);
    1886     else if (degBuf0 > 0)
    1887       uIZeroJ= mulMod (buf[0], bufFactors[1], MOD);
    1888     else if (degBuf1 > 0)
    1889       uIZeroJ= mulMod (bufFactors[0], buf[1], MOD);
    1890     else
    1891       uIZeroJ= 0;
    1892     Pi [0] += xToJ*uIZeroJ;
    1893 
    1894   CFArray tmp= CFArray (factors.length() - 1);
    1895   for (k= 0; k < factors.length() - 1; k++)
    1896     tmp[k]= 0;
    1897   CFIterator one, two;
    1898   one= bufFactors [0];
    1899   two= bufFactors [1];
    1900   if (degBuf0 > 0 && degBuf1 > 0)
    1901   {
    1902     while (one.exp() > j) one++;
    1903     while (two.exp() > j) two++;
    1904     for (k= 1; k <= (int) ceil (j/2.0); k++)
    1905     { 
    1906       if (k != j - k + 1)
    1907       {
    1908         if (one.exp() == j - k + 1 && two.exp() == j - k + 1)
    1909         {
    1910           tmp[0] += mulMod ((bufFactors[0] [k] + one.coeff()),
    1911                     (bufFactors[1] [k] + two.coeff()), MOD) - M (k + 1, 1) -
    1912                     M (j - k + 2, 1);
    1913           one++;
    1914           two++;
    1915         }
    1916         else if (one.exp() == j - k + 1)
    1917         {
    1918           tmp[0] += mulMod ((bufFactors[0] [k] + one.coeff()),
    1919                     bufFactors[1] [k], MOD) - M (k + 1, 1);
    1920           one++;
    1921         }
    1922         else if (two.exp() == j - k + 1)
    1923         {
    1924           tmp[0] += mulMod (bufFactors[0] [k], (bufFactors[1] [k] +
    1925                     two.coeff()), MOD) - M (k + 1, 1);
    1926           two++;
    1927         }
    1928       }
    1929       else
    1930       {
    1931         tmp[0] += M (k + 1, 1);
    1932       }
    1933     }
    1934    
    1935     if (j + 2 <= M.rows())
    1936     {
    1937       if (degBuf0 >= j + 1 && degBuf1 >= j + 1)
    1938         tmp [0] += mulMod ((bufFactors [0] [j + 1]+ bufFactors [0] [0]),
    1939                            (bufFactors [1] [j + 1] + bufFactors [1] [0]), MOD)
    1940                    - M(1,1) - M (j + 2,1);
    1941       else if (degBuf0 >= j + 1)
    1942         tmp[0] += mulMod (bufFactors [0] [j+1], bufFactors [1] [0], MOD);
    1943       else if (degBuf1 >= j + 1)
    1944         tmp[0] += mulMod (bufFactors [0] [0], bufFactors [1] [j + 1], MOD);
    1945     }
    1946   }
    1947   Pi [0] += tmp[0]*xToJ*F.mvar();
    1948 
    1949   // update Pi [l]
    1950   int degPi, degBuf;
    1951   for (int l= 1; l < factors.length() - 1; l++)
    1952   {
    1953     degPi= degree (Pi [l - 1], x);
    1954     degBuf= degree (bufFactors[l + 1], x);
    1955     if (degPi > 0 && degBuf > 0)
    1956       M (j + 1, l + 1)= mulMod (Pi [l - 1] [j], bufFactors[l + 1] [j], MOD); 
    1957     if (j == 1)
    1958     {
    1959       if (degPi > 0 && degBuf > 0)
    1960         Pi [l] += xToJ*(mulMod ((Pi [l - 1] [0] + Pi [l - 1] [j]),
    1961                   (bufFactors[l + 1] [0] + buf[l + 1]), MOD) - M (j + 1, l +1)-
    1962                   M (1, l + 1));
    1963       else if (degPi > 0)
    1964         Pi [l] += xToJ*(mulMod (Pi [l - 1] [j], bufFactors[l + 1], MOD));
    1965       else if (degBuf > 0)
    1966         Pi [l] += xToJ*(mulMod (Pi [l - 1], buf[l + 1], MOD));
    1967     }
    1968     else
    1969     {
    1970       if (degPi > 0 && degBuf > 0)
    1971       {
    1972         uIZeroJ= mulMod (uIZeroJ, bufFactors [l + 1] [0], MOD);
    1973         uIZeroJ += mulMod (Pi [l - 1] [0], buf [l + 1], MOD);
    1974       }
    1975       else if (degPi > 0)
    1976         uIZeroJ= mulMod (uIZeroJ, bufFactors [l + 1], MOD);
    1977       else if (degBuf > 0)
    1978       {
    1979         uIZeroJ= mulMod (uIZeroJ, bufFactors [l + 1] [0], MOD);
    1980         uIZeroJ += mulMod (Pi [l - 1], buf[l + 1], MOD);
    1981       }
    1982       Pi[l] += xToJ*uIZeroJ;
    1983     }
    1984     one= bufFactors [l + 1];
    1985     two= Pi [l - 1];
    1986     if (two.exp() == j + 1)
    1987     {
    1988       if (degBuf > 0 && degPi > 0)
    1989       {
    1990           tmp[l] += mulMod (two.coeff(), bufFactors[l + 1][0], MOD);
    1991           two++;
    1992       }
    1993       else if (degPi > 0)
    1994       {     
    1995           tmp[l] += mulMod (two.coeff(), bufFactors[l + 1], MOD);
    1996           two++;
    1997       }
    1998     }
    1999     if (degBuf > 0 && degPi > 0)
    2000     {
    2001       for (k= 1; k <= (int) ceil (j/2.0); k++)
    2002       { 
    2003         if (k != j - k + 1)
    2004         {
    2005           if (one.exp() == j - k + 1 && two.exp() == j - k + 1)
    2006           {
    2007             tmp[l] += mulMod ((bufFactors[l + 1] [k] + one.coeff()),
    2008                       (Pi[l - 1] [k] + two.coeff()), MOD) - M (k + 1, l + 1) -
    2009                       M (j - k + 2, l + 1);
    2010             one++;
    2011             two++;
    2012           }
    2013           else if (one.exp() == j - k + 1)
    2014           {
    2015             tmp[l] += mulMod ((bufFactors[l + 1] [k] + one.coeff()),
    2016                       Pi[l - 1] [k], MOD) - M (k + 1, l + 1);     
    2017             one++;
    2018           }
    2019           else if (two.exp() == j - k + 1)
    2020           {
    2021             tmp[l] += mulMod (bufFactors[l + 1] [k],
    2022                       (Pi[l - 1] [k] + two.coeff()), MOD) - M (k + 1, l + 1);
    2023             two++;
    2024           }
    2025         }
    2026         else
    2027           tmp[l] += M (k + 1, l + 1);
    2028       } 
    2029     } 
    2030     Pi[l] += tmp[l]*xToJ*F.mvar();
    2031   }
    2032 
    2033   return;
    2034 }
    2035 
    20361412CFList
    20371413henselLift23 (const CFList& eval, const CFList& factors, const int* l, CFList&
    … …  
    20401416  CFList buf= factors;
    20411417  int k= 0;
     1418  int liftBound;
    20421419  int liftBoundBivar= l[k];
    20431420  diophant= biDiophantine (eval.getFirst(), buf, liftBoundBivar);
    … …  
    20611438  if (i.hasItem())
    20621439    i++;
    2063   for (; i.hasItem(); i++, k++)
     1440  for (i; i.hasItem(); i++, k++)
    20641441  {
    20651442    Pi [k]= mulMod (Pi [k - 1], i.getItem(), MOD);
    … …  
    21281505  if (i.hasItem())
    21291506    i++;
    2130   for (; i.hasItem(); i++, k++)
     1507  for (i; i.hasItem(); i++, k++)
    21311508  {
    21321509    Pi [k]= mod (Pi [k], xToLOld);
    … …  
    21441521CFList
    21451522henselLift (const CFList& eval, const CFList& factors, const int* l, const int
    2146             lLength, bool sort)
     1523            lLength)
    21471524{
    21481525  CFList diophant;
    21491526  CFList buf= factors;
    21501527  buf.insert (LC (eval.getFirst(), 1));
    2151   if (sort)
    2152     sortList (buf, Variable (1));
     1528  sortList (buf, Variable (1));
    21531529  CFArray Pi;
    21541530  CFMatrix M= CFMatrix (l[1], factors.length());
    … …  
    21771553}
    21781554
    2179 // wrt. Variable (1)
    2180 CanonicalForm replaceLC (const CanonicalForm& F, const CanonicalForm& c)
    2181 {
    2182   if (degree (F, 1) <= 0)
    2183    return c;
    2184   else
    2185   {
    2186     CanonicalForm result= swapvar (F, Variable (F.level() + 1), Variable (1));
    2187     result += (swapvar (c, Variable (F.level() + 1), Variable (1))
    2188               - LC (result))*power (result.mvar(), degree (result));
    2189     return swapvar (result, Variable (F.level() + 1), Variable (1));
    2190   }
    2191 }
    2192 
    2193 // so far only for two factor Hensel lifting
    2194 CFList
    2195 henselLift232 (const CFList& eval, const CFList& factors, const int* l, CFList&
    2196               diophant, CFArray& Pi, CFMatrix& M, const CFList& LCs1,
    2197               const CFList& LCs2)
    2198 {
    2199   CFList buf= factors;
    2200   int k= 0;
    2201   int liftBoundBivar= l[k];
    2202   CFList bufbuf= factors;
    2203   Variable v= Variable (2);
    2204 
    2205   CFList MOD;
    2206   MOD.append (power (Variable (2), liftBoundBivar));
    2207   CFArray bufFactors= CFArray (factors.length());
    2208   k= 0;
    2209   CFListIterator j= eval;
    2210   j++;
    2211   CFListIterator iter1= LCs1;
    2212   CFListIterator iter2= LCs2;
    2213   iter1++;
    2214   iter2++;
    2215   bufFactors[0]= replaceLC (buf.getFirst(), iter1.getItem());
    2216   bufFactors[1]= replaceLC (buf.getLast(), iter2.getItem());
    2217 
    2218   CFListIterator i= buf;
    2219   i++;
    2220   Variable y= j.getItem().mvar();
    2221   if (y.level() != 3)
    2222     y= Variable (3);
    2223 
    2224   Pi[0]= mod (Pi[0], power (v, liftBoundBivar));
    2225   M (1, 1)= Pi[0];
    2226   if (degree (bufFactors[0], y) > 0 && degree (bufFactors [1], y) > 0)
    2227     Pi [0] += (mulMod (bufFactors [0] [1], bufFactors[1] [0], MOD) +
    2228                         mulMod (bufFactors [0] [0], bufFactors [1] [1], MOD))*y;
    2229   else if (degree (bufFactors[0], y) > 0)
    2230     Pi [0] += mulMod (bufFactors [0] [1], bufFactors[1], MOD)*y;
    2231   else if (degree (bufFactors[1], y) > 0)
    2232     Pi [0] += mulMod (bufFactors [0], bufFactors[1] [1], MOD)*y;
    2233  
    2234   CFList products;
    2235   for (int i= 0; i < bufFactors.size(); i++)
    2236   {
    2237     if (degree (bufFactors[i], y) > 0)
    2238       products.append (M (1, 1)/bufFactors[i] [0]);
    2239     else
    2240       products.append (M (1, 1)/bufFactors[i]);
    2241   }
    2242  
    2243   for (int d= 1; d < l[1]; d++)
    2244     henselStep2 (j.getItem(), buf, bufFactors, diophant, M, Pi, products, d, MOD);
    2245   CFList result;
    2246   for (k= 0; k < factors.length(); k++)
    2247     result.append (bufFactors[k]);
    2248   return result;
    2249 }
    2250 
    2251 
    2252 CFList
    2253 henselLift2 (const CFList& F, const CFList& factors, const CFList& MOD, CFList&
    2254             diophant, CFArray& Pi, CFMatrix& M, const int lOld, const int
    2255             lNew, const CFList& LCs1, const CFList& LCs2)
    2256 {
    2257   int k= 0;
    2258   CFArray bufFactors= CFArray (factors.length());
    2259   bufFactors[0]= replaceLC (factors.getFirst(), LCs1.getLast());
    2260   bufFactors[1]= replaceLC (factors.getLast(), LCs2.getLast());
    2261   CFList buf= factors;
    2262   Variable y= F.getLast().mvar();
    2263   Variable x= F.getFirst().mvar();
    2264   CanonicalForm xToLOld= power (x, lOld);
    2265   Pi [0]= mod (Pi[0], xToLOld);
    2266   M (1, 1)= Pi [0];
    2267 
    2268   if (degree (bufFactors[0], y) > 0 && degree (bufFactors [1], y) > 0)
    2269     Pi [0] += (mulMod (bufFactors [0] [1], bufFactors[1] [0], MOD) +
    2270                         mulMod (bufFactors [0] [0], bufFactors [1] [1], MOD))*y;
    2271   else if (degree (bufFactors[0], y) > 0)
    2272     Pi [0] += mulMod (bufFactors [0] [1], bufFactors[1], MOD)*y;
    2273   else if (degree (bufFactors[1], y) > 0)
    2274     Pi [0] += mulMod (bufFactors [0], bufFactors[1] [1], MOD)*y;
    2275  
    2276   CFList products;
    2277   for (int i= 0; i < bufFactors.size(); i++)
    2278   {
    2279     if (degree (bufFactors[i], y) > 0)
    2280     {
    2281       ASSERT (fdivides (bufFactors[i] [0], M (1, 1)), "expected exact division");
    2282       products.append (M (1, 1)/bufFactors[i] [0]);
    2283     }
    2284     else
    2285     {
    2286       ASSERT (fdivides (bufFactors[i], M (1, 1)), "expected exact division");
    2287       products.append (M (1, 1)/bufFactors[i]);
    2288     }
    2289   }
    2290  
    2291   for (int d= 1; d < lNew; d++)
    2292     henselStep2 (F.getLast(), buf, bufFactors, diophant, M, Pi, products, d, MOD);
    2293   CFList result;
    2294   for (k= 0; k < factors.length(); k++)
    2295     result.append (bufFactors[k]);
    2296   return result;
    2297 }
    2298 
    2299 // so far only for two! factor Hensel lifting
    2300 CFList
    2301 henselLift2 (const CFList& eval, const CFList& factors, const int* l, const int
    2302              lLength, bool sort, const CFList& LCs1, const CFList& LCs2,
    2303              const CFArray& Pi, const CFList& diophant)
    2304 {
    2305   CFList bufDiophant= diophant;
    2306   CFList buf= factors;
    2307   if (sort)
    2308     sortList (buf, Variable (1));
    2309   CFArray bufPi= Pi;
    2310   CFMatrix M= CFMatrix (l[1], factors.length());
    2311   CFList result= henselLift232 (eval, buf, l, bufDiophant, bufPi, M, LCs1, LCs2);
    2312   if (eval.length() == 2)
    2313     return result;
    2314   CFList MOD;
    2315   for (int i= 0; i < 2; i++)
    2316     MOD.append (power (Variable (i + 2), l[i]));
    2317   CFListIterator j= eval;
    2318   j++;
    2319   CFList bufEval;
    2320   bufEval.append (j.getItem());
    2321   j++;
    2322   CFListIterator jj= LCs1;
    2323   CFListIterator jjj= LCs2;
    2324   CFList bufLCs1, bufLCs2;
    2325   jj++, jjj++;
    2326   bufLCs1.append (jj.getItem());
    2327   bufLCs2.append (jjj.getItem());
    2328   jj++, jjj++;
    2329  
    2330   for (int i= 2; i <= lLength && j.hasItem(); i++, j++, jj++, jjj++)
    2331   {
    2332     bufEval.append (j.getItem());
    2333     bufLCs1.append (jj.getItem());
    2334     bufLCs2.append (jjj.getItem());
    2335     M= CFMatrix (l[i], factors.length());
    2336     result= henselLift2 (bufEval, result, MOD, bufDiophant, bufPi, M, l[i - 1],
    2337                          l[i], bufLCs1, bufLCs2);
    2338     MOD.append (power (Variable (i + 2), l[i]));
    2339     bufEval.removeFirst();
    2340     bufLCs1.removeFirst();
    2341     bufLCs2.removeFirst();
    2342   }
    2343   return result;
    2344 }
    2345 
    23461555#endif
    23471556/* HAVE_NTL */

  • factory/facHensel.h

    r7b8a416 rc4f4fd  
    364364           );
    365365
    366 /// two factor Hensel lifting from univariate to bivariate, factors need not to
    367 /// be monic
    368 void
    369 henselLift122 (const CanonicalForm& F,///< [in] a bivariate poly
    370                CFList& factors,       ///< [in, out] a list of univariate polys
    371                                       ///< lifted factors
    372                int l,                 ///< [in] lift bound
    373                CFArray& Pi,           ///< [in, out] stores intermediate results
    374                CFList& diophant,      ///< [in, out] result of diophantine
    375                CFMatrix& M,           ///< [in, out] stores intermediate results
    376                const CFArray& LCs,    ///< [in] leading coefficients
    377                bool sort              ///< [in] if true factors are sorted by
    378                                       ///< their degree
    379               );
    380      
    381 
    382 /// two factor Hensel lifting from bivariate to multivariate, factors need not
    383 /// to be monic
    384 ///
    385 /// @return @a henselLift122 returns a list of lifted factors   
    386 CFList
    387 henselLift2 (const CFList& eval,   ///< [in] a list of polynomials the last
    388                                    ///< element is a compressed multivariate
    389                                    ///< poly, last but one element equals the
    390                                    ///< last elements modulo its main variable
    391                                    ///< and so on. The first element is a
    392                                    ///< compressed bivariate poly.
    393              const CFList& factors,///< [in] bivariate factors
    394              const int* l,         ///< [in] lift bounds
    395              const int lLength,    ///< [in] length of l
    396              bool sort,            ///< [in] if true factors are sorted by
    397                                    ///< their degree in Variable(1)
    398              const CFList& LCs1,   ///< [in] a list of evaluated LC of first
    399                                    ///< factor
    400              const CFList& LCs2,   ///< [in] a list of evaluated LC of second
    401                                    ///< factor
    402              const CFArray& Pi,    ///< [in] intermediate result
    403              const CFList& diophant///< [in] result of diophantine
    404             );
    405366#endif
    406367/* FAC_HENSEL_H */

  • factory/ffreval.h

    r7b8a416 rc4f4fd  
    66#include "canonicalform.h"
    77#include "cf_eval.h"
    8 #include "cf_reval.h"
    98
    109
    … …  
    1817    {
    1918      for( int i=min; i<=max; i++ )
    20         values[i] = start[i] = 0;   
    21     }
    22     FFREvaluation( int min, int max, const AlgExtRandomF & sample ) : REvaluation( min, max, sample ), start( min, max )
    23     {
    24       for( int i=min; i<=max; i++ )
    25         values[i] = start[i] = 0;
    26     }
    27     FFREvaluation( int min, int max, const GFRandom & sample ) : REvaluation( min, max, sample ), start( min, max )
    28     {
    29       for( int i=min; i<=max; i++ )
    30         values[i] = start[i] = 0;
     19        values[i] = start[i] = gen->generate();  //generate random point
     20
     21      nextpoint();
    3122    }
    3223    FFREvaluation& operator= ( const FFREvaluation & e );
Note: See TracChangeset for help on using the changeset viewer.