source: git/factory/cf_gcd.cc @ 2488dc3

spielwiese
Last change on this file since 2488dc3 was 2488dc3, checked in by Martin Lee <martinlee84@…>, 12 years ago
chg: use other GCD if input is not dense in chinrem_gcd
  • Property mode set to 100644
File size: 33.2 KB
RevLine 
[493c477]1/* emacs edit mode for this file is -*- C++ -*- */
[341696]2/* $Id$ */
[9bab9f]3
[e4fe2b]4#include "config.h"
[ab4548f]5
[650f2d8]6#include "cf_assert.h"
[93b061]7#include "debug.h"
8
[9bab9f]9#include "cf_defs.h"
10#include "canonicalform.h"
11#include "cf_iter.h"
12#include "cf_reval.h"
[edb4893]13#include "cf_primes.h"
[fbefc9]14#include "cf_algorithm.h"
[2072126]15#include "cf_factory.h"
[f63dbca]16#include "fac_util.h"
[6db552]17#include "templates/ftmpl_functions.h"
[bb82f0]18#include "algext.h"
[10af64]19#include "cf_gcd_smallp.h"
[6e2ef0e]20#include "cf_map_ext.h"
21#include "cf_util.h"
[d990001]22#include "gfops.h"
[edb4893]23
[f11d7b]24#ifdef HAVE_NTL
[034eec]25#include <NTL/ZZX.h>
[f11d7b]26#include "NTLconvert.h"
[a7ec94]27bool isPurePoly(const CanonicalForm & );
28static CanonicalForm gcd_univar_ntl0( const CanonicalForm &, const CanonicalForm & );
[c4f4fd]29static CanonicalForm gcd_univar_ntlp( const CanonicalForm &, const CanonicalForm & );
[f11d7b]30#endif
31
[7e8c9e]32#ifdef HAVE_FLINT
33#include "FLINTconvert.h"
34static CanonicalForm gcd_univar_flint0 (const CanonicalForm &, const CanonicalForm &);
35static CanonicalForm gcd_univar_flintp (const CanonicalForm &, const CanonicalForm &);
36#endif
37
[a7ec94]38static CanonicalForm cf_content ( const CanonicalForm &, const CanonicalForm & );
39static void cf_prepgcd( const CanonicalForm &, const CanonicalForm &, int &, int &, int & );
[edb4893]40
[27bb97f]41void out_cf(const char *s1,const CanonicalForm &f,const char *s2);
[6f62c3]42
[110718]43CanonicalForm chinrem_gcd(const CanonicalForm & FF,const CanonicalForm & GG);
[f4b180]44
[f63dbca]45bool
46gcd_test_one ( const CanonicalForm & f, const CanonicalForm & g, bool swap )
[9bab9f]47{
48    int count = 0;
49    // assume polys have same level;
[6e2ef0e]50
51    Variable v= Variable (1);
52    bool algExtension= (hasFirstAlgVar (f, v) || hasFirstAlgVar (g, v));
[f63dbca]53    CanonicalForm lcf, lcg;
[6f62c3]54    if ( swap )
55    {
[150dc8]56        lcf = swapvar( LC( f ), Variable(1), f.mvar() );
57        lcg = swapvar( LC( g ), Variable(1), f.mvar() );
[f63dbca]58    }
[6f62c3]59    else
60    {
[150dc8]61        lcf = LC( f, Variable(1) );
62        lcg = LC( g, Variable(1) );
[f63dbca]63    }
[6e2ef0e]64
[f63dbca]65    CanonicalForm F, G;
[6f62c3]66    if ( swap )
67    {
[150dc8]68        F=swapvar( f, Variable(1), f.mvar() );
69        G=swapvar( g, Variable(1), g.mvar() );
[f63dbca]70    }
[6f62c3]71    else
72    {
[150dc8]73        F = f;
74        G = g;
[f63dbca]75    }
[6e2ef0e]76
77    #define TEST_ONE_MAX 50
78    int p= getCharacteristic();
79    bool passToGF= false;
80    int k= 1;
81    if (p > 0 && p < TEST_ONE_MAX && CFFactory::gettype() != GaloisFieldDomain && !algExtension)
82    {
83      if (p == 2)
84        setCharacteristic (2, 6, 'Z');
85      else if (p == 3)
86        setCharacteristic (3, 4, 'Z');
87      else if (p == 5 || p == 7)
88        setCharacteristic (p, 3, 'Z');
89      else
90        setCharacteristic (p, 2, 'Z');
91      passToGF= true;
92    }
93    else if (p > 0 && CFFactory::gettype() == GaloisFieldDomain && ipower (p , getGFDegree()) < TEST_ONE_MAX)
94    {
95      k= getGFDegree();
96      if (ipower (p, 2*k) > TEST_ONE_MAX)
97        setCharacteristic (p, 2*k, gf_name);
98      else
99        setCharacteristic (p, 3*k, gf_name);
100      F= GFMapUp (F, k);
101      G= GFMapUp (G, k);
102      lcf= GFMapUp (lcf, k);
103      lcg= GFMapUp (lcg, k);
104    }
105    else if (p > 0 && p < TEST_ONE_MAX && algExtension)
106    {
107      bool extOfExt= false;
[d990001]108#ifdef HAVE_NTL
[6e2ef0e]109      int d= degree (getMipo (v));
110      CFList source, dest;
111      Variable v2;
112      CanonicalForm primElem, imPrimElem;
113      if (p == 2 && d < 6)
114      {
115        zz_p::init (p);
116        bool primFail= false;
117        Variable vBuf;
118        primElem= primitiveElement (v, vBuf, primFail);
119        ASSERT (!primFail, "failure in integer factorizer");
120        if (d < 3)
121        {
122          zz_pX NTLIrredpoly;
123          BuildIrred (NTLIrredpoly, d*3);
124          CanonicalForm newMipo= convertNTLzzpX2CF (NTLIrredpoly, Variable (1));
125          v2= rootOf (newMipo);
126        }
127        else
128        {
129          zz_pX NTLIrredpoly;
130          BuildIrred (NTLIrredpoly, d*2);
131          CanonicalForm newMipo= convertNTLzzpX2CF (NTLIrredpoly, Variable (1));
132          v2= rootOf (newMipo);
133        }
134        imPrimElem= mapPrimElem (primElem, v, v2);
135        extOfExt= true;
136      }
137      else if ((p == 3 && d < 4) || ((p == 5 || p == 7) && d < 3))
138      {
139        zz_p::init (p);
140        bool primFail= false;
141        Variable vBuf;
142        primElem= primitiveElement (v, vBuf, primFail);
143        ASSERT (!primFail, "failure in integer factorizer");
144        zz_pX NTLIrredpoly;
145        BuildIrred (NTLIrredpoly, d*2);
146        CanonicalForm newMipo= convertNTLzzpX2CF (NTLIrredpoly, Variable (1));
147        v2= rootOf (newMipo);
148        imPrimElem= mapPrimElem (primElem, v, v2);
149        extOfExt= true;
150      }
151      if (extOfExt)
152      {
153        F= mapUp (F, v, v2, primElem, imPrimElem, source, dest);
154        G= mapUp (G, v, v2, primElem, imPrimElem, source, dest);
155        lcf= mapUp (lcf, v, v2, primElem, imPrimElem, source, dest);
156        lcg= mapUp (lcg, v, v2, primElem, imPrimElem, source, dest);
157        v= v2;
158      }
[d990001]159#endif
[6e2ef0e]160    }
161
162    CFRandom * sample;
163    if ((!algExtension && p > 0) || p == 0)
164      sample = CFRandomFactory::generate();
165    else
166      sample = AlgExtRandomF (v).clone();
167
168    REvaluation e( 2, tmax( f.level(), g.level() ), *sample );
169    delete sample;
170
171    if (passToGF)
172    {
173      lcf= lcf.mapinto();
174      lcg= lcg.mapinto();
175    }
176
177    CanonicalForm eval1, eval2;
178    if (passToGF)
179    {
180      eval1= e (lcf);
181      eval2= e (lcg);
182    }
183    else
184    {
185      eval1= e (lcf);
186      eval2= e (lcg);
187    }
188
189    while ( ( eval1.isZero() || eval2.isZero() ) && count < TEST_ONE_MAX )
190    {
191        e.nextpoint();
192        count++;
193        eval1= e (lcf);
194        eval2= e (lcg);
195    }
196    if ( count >= TEST_ONE_MAX )
197    {
198        if (passToGF)
199          setCharacteristic (p);
200        if (k > 1)
201          setCharacteristic (p, k, gf_name);
202        return false;
203    }
204
205
206    if (passToGF)
207    {
208      F= F.mapinto();
209      G= G.mapinto();
210      eval1= e (F);
211      eval2= e (G);
212    }
213    else
214    {
215      eval1= e (F);
216      eval2= e (G);
217    }
218
219    if ( eval1.taildegree() > 0 && eval2.taildegree() > 0 )
220    {
221        if (passToGF)
222          setCharacteristic (p);
223        if (k > 1)
224          setCharacteristic (p, k, gf_name);
[150dc8]225        return false;
[6e2ef0e]226    }
227
228    CanonicalForm c= gcd (eval1, eval2);
229    bool result= c.degree() < 1;
230
231    if (passToGF)
232      setCharacteristic (p);
233    if (k > 1)
234      setCharacteristic (p, k, gf_name);
235    return result;
[9bab9f]236}
237
[dd3e561]238//{{{ static CanonicalForm icontent ( const CanonicalForm & f, const CanonicalForm & c )
239//{{{ docu
240//
241// icontent() - return gcd of c and all coefficients of f which
242//   are in a coefficient domain.
243//
244// Used by icontent().
245//
246//}}}
[9bab9f]247static CanonicalForm
248icontent ( const CanonicalForm & f, const CanonicalForm & c )
249{
[c30347]250    if ( f.inBaseDomain() )
251    {
252      if (c.isZero()) return abs(f);
253      return bgcd( f, c );
254    }
[ef20c7]255    //else if ( f.inCoeffDomain() )
256    //   return gcd(f,c);
[c30347]257    else
258    {
[150dc8]259        CanonicalForm g = c;
260        for ( CFIterator i = f; i.hasTerms() && ! g.isOne(); i++ )
261            g = icontent( i.coeff(), g );
262        return g;
[9bab9f]263    }
264}
[dd3e561]265//}}}
[9bab9f]266
[dd3e561]267//{{{ CanonicalForm icontent ( const CanonicalForm & f )
268//{{{ docu
269//
270// icontent() - return gcd over all coefficients of f which are
271//   in a coefficient domain.
272//
273//}}}
[9bab9f]274CanonicalForm
275icontent ( const CanonicalForm & f )
276{
277    return icontent( f, 0 );
278}
[dd3e561]279//}}}
[9bab9f]280
[dd3e561]281//{{{ CanonicalForm extgcd ( const CanonicalForm & f, const CanonicalForm & g, CanonicalForm & a, CanonicalForm & b )
282//{{{ docu
283//
284// extgcd() - returns polynomial extended gcd of f and g.
285//
286// Returns gcd(f, g) and a and b sucht that f*a+g*b=gcd(f, g).
287// The gcd is calculated using an extended euclidean polynomial
288// remainder sequence, so f and g should be polynomials over an
289// euclidean domain.  Normalizes result.
290//
291// Note: be sure that f and g have the same level!
292//
293//}}}
[9bab9f]294CanonicalForm
295extgcd ( const CanonicalForm & f, const CanonicalForm & g, CanonicalForm & a, CanonicalForm & b )
296{
[e9a5b62]297  if (f.isZero())
298  {
299    a= 0;
300    b= 1;
301    return g;
302  }
303  else if (g.isZero())
304  {
305    a= 1;
306    b= 0;
307    return f;
308  }
[034eec]309#ifdef HAVE_NTL
[963057]310  if (isOn(SW_USE_NTL_GCD_P) && ( getCharacteristic() > 0 ) && (CFFactory::gettype() != GaloisFieldDomain)
[a86cda]311  &&  (f.level()==g.level()) && isPurePoly(f) && isPurePoly(g))
[034eec]312  {
[c6eecb]313    if (fac_NTL_char!=getCharacteristic())
314    {
315      fac_NTL_char=getCharacteristic();
316      #ifdef NTL_ZZ
317      ZZ r;
318      r=getCharacteristic();
319      ZZ_pContext ccc(r);
320      #else
321      zz_pContext ccc(getCharacteristic());
322      #endif
323      ccc.restore();
324      #ifdef NTL_ZZ
325      ZZ_p::init(r);
326      #else
327      zz_p::init(getCharacteristic());
328      #endif
329    }
330    #ifdef NTL_ZZ
331    ZZ_pX F1=convertFacCF2NTLZZpX(f);
332    ZZ_pX G1=convertFacCF2NTLZZpX(g);
333    ZZ_pX R;
334    ZZ_pX A,B;
335    XGCD(R,A,B,F1,G1);
336    a=convertNTLZZpX2CF(A,f.mvar());
337    b=convertNTLZZpX2CF(B,f.mvar());
338    return convertNTLZZpX2CF(R,f.mvar());
339    #else
[034eec]340    zz_pX F1=convertFacCF2NTLzzpX(f);
341    zz_pX G1=convertFacCF2NTLzzpX(g);
342    zz_pX R;
343    zz_pX A,B;
344    XGCD(R,A,B,F1,G1);
345    a=convertNTLzzpX2CF(A,f.mvar());
346    b=convertNTLzzpX2CF(B,f.mvar());
347    return convertNTLzzpX2CF(R,f.mvar());
[c6eecb]348    #endif
[034eec]349  }
[a86cda]350  if (isOn(SW_USE_NTL_GCD_0) && ( getCharacteristic() ==0)
351  && (f.level()==g.level()) && isPurePoly(f) && isPurePoly(g))
352  {
353    CanonicalForm fc=bCommonDen(f);
354    CanonicalForm gc=bCommonDen(g);
355    ZZX F1=convertFacCF2NTLZZX(f*fc);
356    ZZX G1=convertFacCF2NTLZZX(g*gc);
357    ZZX R=GCD(F1,G1);
358    CanonicalForm r=convertNTLZZX2CF(R,f.mvar());
359    ZZ RR;
360    ZZX A,B;
361    if (r.inCoeffDomain())
362    {
363      XGCD(RR,A,B,F1,G1,1);
364      CanonicalForm rr=convertZZ2CF(RR);
365      ASSERT (!rr.isZero(), "NTL:XGCD failed");
366      a=convertNTLZZX2CF(A,f.mvar())*fc/rr;
367      b=convertNTLZZX2CF(B,f.mvar())*gc/rr;
368      return CanonicalForm(1);
369    }
370    else
371    {
372      fc=bCommonDen(f);
373      gc=bCommonDen(g);
374      F1=convertFacCF2NTLZZX(f*fc/r);
375      G1=convertFacCF2NTLZZX(g*gc/r);
376      XGCD(RR,A,B,F1,G1,1);
377      a=convertNTLZZX2CF(A,f.mvar())*fc;
378      b=convertNTLZZX2CF(B,f.mvar())*gc;
379      CanonicalForm rr=convertZZ2CF(RR);
380      ASSERT (!rr.isZero(), "NTL:XGCD failed");
381      r*=rr;
382      if ( r.sign() < 0 ) { r= -r; a= -a; b= -b; }
383      return r;
384    }
385  }
[034eec]386#endif
[a86cda]387  // may contain bug in the co-factors, see track 107
[034eec]388  CanonicalForm contf = content( f );
389  CanonicalForm contg = content( g );
[9bab9f]390
[034eec]391  CanonicalForm p0 = f / contf, p1 = g / contg;
392  CanonicalForm f0 = 1, f1 = 0, g0 = 0, g1 = 1, q, r;
[9bab9f]393
[c6eecb]394  while ( ! p1.isZero() )
395  {
[034eec]396      divrem( p0, p1, q, r );
397      p0 = p1; p1 = r;
398      r = g0 - g1 * q;
399      g0 = g1; g1 = r;
400      r = f0 - f1 * q;
401      f0 = f1; f1 = r;
402  }
403  CanonicalForm contp0 = content( p0 );
404  a = f0 / ( contf * contp0 );
405  b = g0 / ( contg * contp0 );
406  p0 /= contp0;
[c6eecb]407  if ( p0.sign() < 0 )
408  {
[034eec]409      p0 = -p0;
410      a = -a;
411      b = -b;
412  }
413  return p0;
[9bab9f]414}
[dd3e561]415//}}}
[9bab9f]416
[a7ec94]417//{{{ static CanonicalForm balance ( const CanonicalForm & f, const CanonicalForm & q )
418//{{{ docu
419//
420// balance() - map f from positive to symmetric representation
421//   mod q.
422//
423// This makes sense for univariate polynomials over Z only.
424// q should be an integer.
425//
426// Used by gcd_poly_univar0().
427//
428//}}}
[edb4893]429static CanonicalForm
[a7ec94]430balance ( const CanonicalForm & f, const CanonicalForm & q )
[edb4893]431{
[a7ec94]432    Variable x = f.mvar();
433    CanonicalForm result = 0, qh = q / 2;
434    CanonicalForm c;
435    CFIterator i;
436    for ( i = f; i.hasTerms(); i++ ) {
437        c = mod( i.coeff(), q );
438        if ( c > qh )
439            result += power( x, i.exp() ) * (c - q);
440        else
441            result += power( x, i.exp() ) * c;
[edb4893]442    }
[a7ec94]443    return result;
444}
445//}}}
446
[01e8874]447static CanonicalForm gcd_poly_univar0( const CanonicalForm & F, const CanonicalForm & G, bool primitive )
[a7ec94]448{
[f11d7b]449  CanonicalForm f, g, c, cg, cl, BB, B, M, q, Dp, newD, D, newq;
[01e8874]450  int p, i;
[f11d7b]451
452  if ( primitive )
453  {
454    f = F;
455    g = G;
456    c = 1;
457  }
458  else
459  {
460    CanonicalForm cF = content( F ), cG = content( G );
461    f = F / cF;
462    g = G / cG;
463    c = bgcd( cF, cG );
464  }
465  cg = gcd( f.lc(), g.lc() );
466  cl = ( f.lc() / cg ) * g.lc();
[93b061]467//     B = 2 * cg * tmin(
[150dc8]468//         maxnorm(f)*power(CanonicalForm(2),f.degree())*isqrt(f.degree()+1),
469//         maxnorm(g)*power(CanonicalForm(2),g.degree())*isqrt(g.degree()+1)
470//         )+1;
[f11d7b]471  M = tmin( maxNorm(f), maxNorm(g) );
472  BB = power(CanonicalForm(2),tmin(f.degree(),g.degree()))*M;
473  q = 0;
474  i = cf_getNumSmallPrimes() - 1;
475  while ( true )
476  {
477    B = BB;
478    while ( i >= 0 && q < B )
479    {
480      p = cf_getSmallPrime( i );
481      i--;
482      while ( i >= 0 && mod( cl, p ) == 0 )
483      {
484        p = cf_getSmallPrime( i );
485        i--;
486      }
487      setCharacteristic( p );
488      Dp = gcd( mapinto( f ), mapinto( g ) );
489      Dp = ( Dp / Dp.lc() ) * mapinto( cg );
490      setCharacteristic( 0 );
491      if ( Dp.degree() == 0 )
492        return c;
493      if ( q.isZero() )
494      {
495        D = mapinto( Dp );
496        q = p;
497        B = power(CanonicalForm(2),D.degree())*M+1;
498      }
499      else
500      {
501        if ( Dp.degree() == D.degree() )
502        {
503          chineseRemainder( D, q, mapinto( Dp ), p, newD, newq );
504          q = newq;
505          D = newD;
[150dc8]506        }
[f11d7b]507        else if ( Dp.degree() < D.degree() )
508        {
509          // all previous p's are bad primes
510          q = p;
511          D = mapinto( Dp );
512          B = power(CanonicalForm(2),D.degree())*M+1;
[150dc8]513        }
[f11d7b]514        // else p is a bad prime
515      }
516    }
517    if ( i >= 0 )
518    {
519      // now balance D mod q
520      D = pp( balance( D, q ) );
[ebc602]521      if ( fdivides( D, f ) && fdivides( D, g ) )
[f11d7b]522        return D * c;
523      else
524        q = 0;
[edb4893]525    }
[f11d7b]526    else
[a7ec94]527      return gcd_poly( F, G );
[f11d7b]528    DEBOUTLN( cerr, "another try ..." );
529  }
[edb4893]530}
531
[c4f4fd]532static CanonicalForm
533gcd_poly_p( const CanonicalForm & f, const CanonicalForm & g )
534{
535    CanonicalForm pi, pi1;
536    CanonicalForm C, Ci, Ci1, Hi, bi, pi2;
537    bool bpure;
538    int delta = degree( f ) - degree( g );
539
540    if ( delta >= 0 )
541    {
542        pi = f; pi1 = g;
543    }
544    else
545    {
546        pi = g; pi1 = f; delta = -delta;
547    }
548    Ci = content( pi ); Ci1 = content( pi1 );
549    pi1 = pi1 / Ci1; pi = pi / Ci;
550    C = gcd( Ci, Ci1 );
551    if ( !( pi.isUnivariate() && pi1.isUnivariate() ) )
552    {
553        if ( gcd_test_one( pi1, pi, true ) )
554        {
555          C=abs(C);
556          //out_cf("GCD:",C,"\n");
557          return C;
558        }
559        bpure = false;
560    }
561    else
562    {
563        bpure = isPurePoly(pi) && isPurePoly(pi1);
[7e8c9e]564#ifdef HAVE_FLINT
565        if (bpure && (CFFactory::gettype() != GaloisFieldDomain))
566          return gcd_univar_flintp(pi,pi1)*C;
[7cb5590]567#else
568#ifdef HAVE_NTL
[c4f4fd]569        if ( isOn(SW_USE_NTL_GCD_P) && bpure && (CFFactory::gettype() != GaloisFieldDomain))
570            return gcd_univar_ntlp(pi, pi1 ) * C;
[7cb5590]571#endif
[c4f4fd]572#endif
573    }
574    Variable v = f.mvar();
575    Hi = power( LC( pi1, v ), delta );
576    if ( (delta+1) % 2 )
577        bi = 1;
578    else
579        bi = -1;
[6e2ef0e]580    int maxNumVars= tmax (getNumVars (pi), getNumVars (pi1));
581    CanonicalForm oldPi= pi, oldPi1= pi1;
[c4f4fd]582    while ( degree( pi1, v ) > 0 )
583    {
[6e2ef0e]584        if (!(pi.isUnivariate() && pi1.isUnivariate()))
585        {
586          if (size (pi)/maxNumVars > 500 || size (pi1)/maxNumVars > 500)
587          {
588            On (SW_USE_FF_MOD_GCD);
589            C *= gcd (oldPi, oldPi1);
590            Off (SW_USE_FF_MOD_GCD);
591            return C;
592          }
593        }
[c4f4fd]594        pi2 = psr( pi, pi1, v );
595        pi2 = pi2 / bi;
596        pi = pi1; pi1 = pi2;
[6e2ef0e]597        maxNumVars= tmax (getNumVars (pi), getNumVars (pi1));
[c4f4fd]598        if ( degree( pi1, v ) > 0 )
599        {
600            delta = degree( pi, v ) - degree( pi1, v );
601            if ( (delta+1) % 2 )
602                bi = LC( pi, v ) * power( Hi, delta );
603            else
604                bi = -LC( pi, v ) * power( Hi, delta );
605            Hi = power( LC( pi1, v ), delta ) / power( Hi, delta-1 );
606        }
607    }
608    if ( degree( pi1, v ) == 0 )
609    {
610      C=abs(C);
611      //out_cf("GCD:",C,"\n");
612      return C;
613    }
614    pi /= content( pi );
615    if ( bpure )
616        pi /= pi.lc();
617    C=abs(C*pi);
618    //out_cf("GCD:",C,"\n");
619    return C;
620}
621
[a7ec94]622static CanonicalForm
623gcd_poly_0( const CanonicalForm & f, const CanonicalForm & g )
624{
625    CanonicalForm pi, pi1;
[df497a]626    CanonicalForm C, Ci, Ci1, Hi, bi, pi2;
[a7ec94]627    int delta = degree( f ) - degree( g );
628
629    if ( delta >= 0 )
630    {
631        pi = f; pi1 = g;
632    }
633    else
634    {
635        pi = g; pi1 = f; delta = -delta;
636    }
[9bab9f]637    Ci = content( pi ); Ci1 = content( pi1 );
638    pi1 = pi1 / Ci1; pi = pi / Ci;
[df497a]639    C = gcd( Ci, Ci1 );
[034eec]640    if ( pi.isUnivariate() && pi1.isUnivariate() )
641    {
[c53fdc]642/*#ifdef HAVE_FLINT
[7e8c9e]643        if (isPurePoly(pi) && isPurePoly(pi1) )
644            return gcd_univar_flint0(pi, pi1 ) * C;
[c53fdc]645#else*/
[7cb5590]646#ifdef HAVE_NTL
[a7ec94]647        if ( isOn(SW_USE_NTL_GCD_0) && isPurePoly(pi) && isPurePoly(pi1) )
648            return gcd_univar_ntl0(pi, pi1 ) * C;
[7cb5590]649#endif
[c53fdc]650//#endif
[a7ec94]651        return gcd_poly_univar0( pi, pi1, true ) * C;
[edb4893]652    }
[034eec]653    else if ( gcd_test_one( pi1, pi, true ) )
654      return C;
[a7ec94]655    Variable v = f.mvar();
[9bab9f]656    Hi = power( LC( pi1, v ), delta );
657    if ( (delta+1) % 2 )
[150dc8]658        bi = 1;
[9bab9f]659    else
[150dc8]660        bi = -1;
[6f62c3]661    while ( degree( pi1, v ) > 0 )
662    {
[150dc8]663        pi2 = psr( pi, pi1, v );
664        pi2 = pi2 / bi;
665        pi = pi1; pi1 = pi2;
[6f62c3]666        if ( degree( pi1, v ) > 0 )
667        {
[150dc8]668            delta = degree( pi, v ) - degree( pi1, v );
669            if ( (delta+1) % 2 )
670                bi = LC( pi, v ) * power( Hi, delta );
671            else
672                bi = -LC( pi, v ) * power( Hi, delta );
673            Hi = power( LC( pi1, v ), delta ) / power( Hi, delta-1 );
674        }
[9bab9f]675    }
676    if ( degree( pi1, v ) == 0 )
[150dc8]677        return C;
[df497a]678    else
[150dc8]679        return C * pp( pi );
[9bab9f]680}
681
[b809a8]682//{{{ CanonicalForm gcd_poly ( const CanonicalForm & f, const CanonicalForm & g )
[dd3e561]683//{{{ docu
684//
685// gcd_poly() - calculate polynomial gcd.
686//
687// This is the dispatcher for polynomial gcd calculation.  We call either
688// ezgcd(), sparsemod() or gcd_poly1() in dependecy on the current
[e88604]689// characteristic and settings of SW_USE_EZGCD.
[dd3e561]690//
691// Used by gcd() and gcd_poly_univar0().
692//
693//}}}
[0b6919]694#if 0
[bfc606]695int si_factor_reminder=1;
[0b6919]696#endif
[b809a8]697CanonicalForm gcd_poly ( const CanonicalForm & f, const CanonicalForm & g )
[f63dbca]698{
[110718]699  CanonicalForm fc, gc, d1;
700  int mp, cc, p1, pe;
701  mp = f.level()+1;
[ed9927]702  bool fc_isUnivariate=f.isUnivariate();
703  bool gc_isUnivariate=g.isUnivariate();
704  bool fc_and_gc_Univariate=fc_isUnivariate && gc_isUnivariate;
705  cf_prepgcd( f, g, cc, p1, pe);
706  if ( cc != 0 )
[110718]707  {
[ed9927]708    if ( cc > 0 )
[abfc3b]709    {
[ed9927]710      fc = replacevar( f, Variable(cc), Variable(mp) );
711      gc = g;
[e074407]712    }
[ed9927]713    else
[110718]714    {
[ed9927]715      fc = replacevar( g, Variable(-cc), Variable(mp) );
716      gc = f;
[110718]717    }
[ed9927]718    return cf_content( fc, gc );
719  }
720// now each appearing variable is in f and g
721  fc = f;
722  gc = g;
723  if ( getCharacteristic() != 0 )
724  {
[2072126]725    #ifdef HAVE_NTL
[e16f7d]726    if ((!fc_and_gc_Univariate) && (isOn( SW_USE_EZGCD_P )))
[49f1f45]727    {
[08daea]728      fc= EZGCD_P (fc, gc);
[c30347]729    }
[10af64]730    else if (isOn(SW_USE_FF_MOD_GCD) && !fc_and_gc_Univariate)
731    {
732      Variable a;
733      if (hasFirstAlgVar (fc, a) || hasFirstAlgVar (gc, a))
734        fc=GCD_Fp_extension (fc, gc, a);
[b5c084]735      else if (CFFactory::gettype() == GaloisFieldDomain)
[10af64]736        fc=GCD_GF (fc, gc);
[b5c084]737      else
738        fc=GCD_small_p (fc, gc);
[10af64]739    }
[efcd2dc]740    else
[2072126]741    #endif
[efcd2dc]742    if ( p1 == fc.level() )
[ed9927]743      fc = gcd_poly_p( fc, gc );
744    else
745    {
746      fc = replacevar( fc, Variable(p1), Variable(mp) );
747      gc = replacevar( gc, Variable(p1), Variable(mp) );
748      fc = replacevar( gcd_poly_p( fc, gc ), Variable(mp), Variable(p1) );
749    }
[110718]750  }
[c30347]751  else if (!fc_and_gc_Univariate)
[110718]752  {
[f7a4e9]753    if ( isOn( SW_USE_EZGCD ) )
754    {
755      fc= ezgcd (fc, gc);
756      /*if ( pe == 1 )
757        fc = ezgcd( fc, gc );
758      else if ( pe > 0 )// no variable at position 1
759      {
760        fc = replacevar( fc, Variable(pe), Variable(1) );
761        gc = replacevar( gc, Variable(pe), Variable(1) );
762        fc = replacevar( ezgcd( fc, gc ), Variable(1), Variable(pe) );
763      }
764      else
765      {
766        pe = -pe;
767        fc = swapvar( fc, Variable(pe), Variable(1) );
768        gc = swapvar( gc, Variable(pe), Variable(1) );
769        fc = swapvar( ezgcd( fc, gc ), Variable(1), Variable(pe) );
770      }*/
771    }
772    else if (
[c30347]773    isOn(SW_USE_CHINREM_GCD)
[ed9927]774    && (isPurePoly_m(fc)) && (isPurePoly_m(gc))
[c30347]775    )
776    {
[ed9927]777    #if 0
778      if ( p1 == fc.level() )
779        fc = chinrem_gcd( fc, gc );
780      else
781      {
782        fc = replacevar( fc, Variable(p1), Variable(mp) );
783        gc = replacevar( gc, Variable(p1), Variable(mp) );
784        fc = replacevar( chinrem_gcd( fc, gc ), Variable(mp), Variable(p1) );
785      }
786    #else
787      fc = chinrem_gcd( fc, gc);
788    #endif
[c30347]789    }
790    else
791    {
[ed9927]792       fc = gcd_poly_0( fc, gc );
[c30347]793    }
[110718]794  }
795  else
796  {
797    fc = gcd_poly_0( fc, gc );
798  }
799  if ( d1.degree() > 0 )
800    fc *= d1;
801  return fc;
[f63dbca]802}
[dd3e561]803//}}}
[93b061]804
[dd3e561]805//{{{ static CanonicalForm cf_content ( const CanonicalForm & f, const CanonicalForm & g )
806//{{{ docu
807//
808// cf_content() - return gcd(g, content(f)).
809//
810// content(f) is calculated with respect to f's main variable.
811//
812// Used by gcd(), content(), content( CF, Variable ).
813//
814//}}}
[9bab9f]815static CanonicalForm
816cf_content ( const CanonicalForm & f, const CanonicalForm & g )
817{
[6f62c3]818    if ( f.inPolyDomain() || ( f.inExtension() && ! getReduce( f.mvar() ) ) )
819    {
[150dc8]820        CFIterator i = f;
821        CanonicalForm result = g;
[6f62c3]822        while ( i.hasTerms() && ! result.isOne() )
823        {
[a7ec94]824            result = gcd( i.coeff(), result );
[150dc8]825            i++;
826        }
827        return result;
[9bab9f]828    }
829    else
[a7ec94]830        return abs( f );
[9bab9f]831}
[dd3e561]832//}}}
[9bab9f]833
[4ea0ab]834//{{{ CanonicalForm content ( const CanonicalForm & f )
835//{{{ docu
836//
837// content() - return content(f) with respect to main variable.
838//
[dd3e561]839// Normalizes result.
840//
[4ea0ab]841//}}}
[9bab9f]842CanonicalForm
843content ( const CanonicalForm & f )
844{
[6f62c3]845    if ( f.inPolyDomain() || ( f.inExtension() && ! getReduce( f.mvar() ) ) )
846    {
[a7ec94]847        CFIterator i = f;
848        CanonicalForm result = abs( i.coeff() );
849        i++;
[6f62c3]850        while ( i.hasTerms() && ! result.isOne() )
851        {
[a7ec94]852            result = gcd( i.coeff(), result );
853            i++;
854        }
855        return result;
856    }
857    else
858        return abs( f );
[9bab9f]859}
[4ea0ab]860//}}}
[9bab9f]861
[dd3e561]862//{{{ CanonicalForm content ( const CanonicalForm & f, const Variable & x )
863//{{{ docu
864//
865// content() - return content(f) with respect to x.
866//
867// x should be a polynomial variable.
868//
869//}}}
[9bab9f]870CanonicalForm
871content ( const CanonicalForm & f, const Variable & x )
872{
[dd3e561]873    ASSERT( x.level() > 0, "cannot calculate content with respect to algebraic variable" );
874    Variable y = f.mvar();
875
876    if ( y == x )
[150dc8]877        return cf_content( f, 0 );
[dd3e561]878    else  if ( y < x )
[150dc8]879        return f;
[9bab9f]880    else
[150dc8]881        return swapvar( content( swapvar( f, y, x ), y ), y, x );
[9bab9f]882}
[dd3e561]883//}}}
[9bab9f]884
[dd3e561]885//{{{ CanonicalForm vcontent ( const CanonicalForm & f, const Variable & x )
886//{{{ docu
887//
888// vcontent() - return content of f with repect to variables >= x.
889//
890// The content is recursively calculated over all coefficients in
891// f having level less than x.  x should be a polynomial
892// variable.
893//
894//}}}
[9bab9f]895CanonicalForm
896vcontent ( const CanonicalForm & f, const Variable & x )
897{
[dd3e561]898    ASSERT( x.level() > 0, "cannot calculate vcontent with respect to algebraic variable" );
899
[9bab9f]900    if ( f.mvar() <= x )
[150dc8]901        return content( f, x );
[9bab9f]902    else {
[150dc8]903        CFIterator i;
904        CanonicalForm d = 0;
905        for ( i = f; i.hasTerms() && ! d.isOne(); i++ )
906            d = gcd( d, vcontent( i.coeff(), x ) );
907        return d;
[9bab9f]908    }
909}
[dd3e561]910//}}}
[9bab9f]911
[4ea0ab]912//{{{ CanonicalForm pp ( const CanonicalForm & f )
913//{{{ docu
914//
915// pp() - return primitive part of f.
916//
[dd3e561]917// Returns zero if f equals zero, otherwise f / content(f).
918//
[4ea0ab]919//}}}
[9bab9f]920CanonicalForm
921pp ( const CanonicalForm & f )
922{
923    if ( f.isZero() )
[150dc8]924        return f;
[9bab9f]925    else
[150dc8]926        return f / content( f );
[9bab9f]927}
[4ea0ab]928//}}}
[9bab9f]929
[ff6222]930CanonicalForm
[9bab9f]931gcd ( const CanonicalForm & f, const CanonicalForm & g )
932{
[a7ec94]933    bool b = f.isZero();
934    if ( b || g.isZero() )
935    {
936        if ( b )
937            return abs( g );
[abfc3b]938        else
[a7ec94]939            return abs( f );
940    }
941    if ( f.inPolyDomain() || g.inPolyDomain() )
942    {
943        if ( f.mvar() != g.mvar() )
944        {
945            if ( f.mvar() > g.mvar() )
946                return cf_content( f, g );
947            else
948                return cf_content( g, f );
949        }
[bb82f0]950        if (isOn(SW_USE_QGCD))
951        {
952          Variable m;
[fc9f44]953          if (
954          (getCharacteristic() == 0) &&
[e6f7ee1]955          (hasFirstAlgVar(f,m) || hasFirstAlgVar(g,m))
[bb82f0]956          )
[fc31bce]957          {
[713bdb]958            bool on_rational = isOn(SW_RATIONAL);
959            CanonicalForm r=QGCD(f,g);
[f06059]960            On(SW_RATIONAL);
[713bdb]961            CanonicalForm cdF = bCommonDen( r );
962            if (!on_rational) Off(SW_RATIONAL);
963            return cdF*r;
[fc31bce]964          }
[bb82f0]965        }
[713bdb]966
[150dc8]967        if ( f.inExtension() && getReduce( f.mvar() ) )
[bb82f0]968            return CanonicalForm(1);
[a7ec94]969        else
970        {
[ebc602]971            if ( fdivides( f, g ) )
[a7ec94]972                return abs( f );
[ebc602]973            else  if ( fdivides( g, f ) )
[a7ec94]974                return abs( g );
975            if ( !( getCharacteristic() == 0 && isOn( SW_RATIONAL ) ) )
976            {
977                CanonicalForm d;
[64a501]978                d = gcd_poly( f, g );
[a7ec94]979                return abs( d );
980            }
981            else
982            {
[56d3c6]983                //printf ("here\n");
[150dc8]984                CanonicalForm cdF = bCommonDen( f );
985                CanonicalForm cdG = bCommonDen( g );
986                Off( SW_RATIONAL );
987                CanonicalForm l = lcm( cdF, cdG );
988                On( SW_RATIONAL );
989                CanonicalForm F = f * l, G = g * l;
990                Off( SW_RATIONAL );
[64a501]991                l = gcd_poly( F, G );
[150dc8]992                On( SW_RATIONAL );
[a7ec94]993                return abs( l );
[150dc8]994            }
995        }
[a7ec94]996    }
997    if ( f.inBaseDomain() && g.inBaseDomain() )
998        return bgcd( f, g );
[9bab9f]999    else
[a7ec94]1000        return 1;
[9bab9f]1001}
1002
[dd3e561]1003//{{{ CanonicalForm lcm ( const CanonicalForm & f, const CanonicalForm & g )
1004//{{{ docu
1005//
1006// lcm() - return least common multiple of f and g.
1007//
1008// The lcm is calculated using the formula lcm(f, g) = f * g / gcd(f, g).
1009//
1010// Returns zero if one of f or g equals zero.
1011//
1012//}}}
[9bab9f]1013CanonicalForm
1014lcm ( const CanonicalForm & f, const CanonicalForm & g )
1015{
[dd3e561]1016    if ( f.isZero() || g.isZero() )
[a7ec94]1017        return 0;
[dd3e561]1018    else
[150dc8]1019        return ( f / gcd( f, g ) ) * g;
[9bab9f]1020}
[dd3e561]1021//}}}
[a7ec94]1022
1023#ifdef HAVE_NTL
1024
1025static CanonicalForm
1026gcd_univar_ntl0( const CanonicalForm & F, const CanonicalForm & G )
1027{
1028    ZZX F1=convertFacCF2NTLZZX(F);
1029    ZZX G1=convertFacCF2NTLZZX(G);
1030    ZZX R=GCD(F1,G1);
1031    return convertNTLZZX2CF(R,F.mvar());
1032}
1033
[c4f4fd]1034static CanonicalForm
1035gcd_univar_ntlp( const CanonicalForm & F, const CanonicalForm & G )
1036{
1037  if (fac_NTL_char!=getCharacteristic())
1038  {
1039    fac_NTL_char=getCharacteristic();
1040    #ifdef NTL_ZZ
1041    ZZ r;
1042    r=getCharacteristic();
1043    ZZ_pContext ccc(r);
1044    #else
1045    zz_pContext ccc(getCharacteristic());
1046    #endif
1047    ccc.restore();
1048    #ifdef NTL_ZZ
1049    ZZ_p::init(r);
1050    #else
1051    zz_p::init(getCharacteristic());
1052    #endif
1053  }
1054  #ifdef NTL_ZZ
1055  ZZ_pX F1=convertFacCF2NTLZZpX(F);
1056  ZZ_pX G1=convertFacCF2NTLZZpX(G);
1057  ZZ_pX R=GCD(F1,G1);
1058  return  convertNTLZZpX2CF(R,F.mvar());
1059  #else
1060  zz_pX F1=convertFacCF2NTLzzpX(F);
1061  zz_pX G1=convertFacCF2NTLzzpX(G);
1062  zz_pX R=GCD(F1,G1);
1063  return  convertNTLzzpX2CF(R,F.mvar());
1064  #endif
1065}
1066
[a7ec94]1067#endif
1068
[7e8c9e]1069#ifdef HAVE_FLINT
1070static CanonicalForm
1071gcd_univar_flintp (const CanonicalForm& F, const CanonicalForm& G)
1072{
1073  nmod_poly_t F1, G1;
1074  convertFacCF2nmod_poly_t (F1, F);
1075  convertFacCF2nmod_poly_t (G1, G);
1076  nmod_poly_gcd (F1, F1, G1);
1077  CanonicalForm result= convertnmod_poly_t2FacCF (F1, F.mvar());
1078  nmod_poly_clear (F1);
1079  nmod_poly_clear (G1);
1080  return result;
1081}
1082
1083static CanonicalForm
1084gcd_univar_flint0( const CanonicalForm & F, const CanonicalForm & G )
1085{
1086  fmpz_poly_t F1, G1;
1087  convertFacCF2Fmpz_poly_t(F1, F);
1088  convertFacCF2Fmpz_poly_t(G1, G);
1089  fmpz_poly_gcd (F1, F1, G1);
1090  CanonicalForm result= convertFmpz_poly_t2FacCF (F1, F.mvar());
1091  fmpz_poly_clear (F1);
1092  fmpz_poly_clear (G1);
1093  return result;
1094}
1095#endif
1096
1097
[a7ec94]1098/*
1099*  compute positions p1 and pe of optimal variables:
1100*    pe is used in "ezgcd" and
1101*    p1 in "gcd_poly1"
1102*/
1103static
1104void optvalues ( const int * df, const int * dg, const int n, int & p1, int &pe )
1105{
1106    int i, o1, oe;
1107    if ( df[n] > dg[n] )
1108    {
1109        o1 = df[n]; oe = dg[n];
1110    }
1111    else
1112    {
1113        o1 = dg[n]; oe = df[n];
1114    }
1115    i = n-1;
1116    while ( i > 0 )
1117    {
1118        if ( df[i] != 0 )
1119        {
1120            if ( df[i] > dg[i] )
1121            {
1122                if ( o1 > df[i]) { o1 = df[i]; p1 = i; }
1123                if ( oe <= dg[i]) { oe = dg[i]; pe = i; }
1124            }
1125            else
1126            {
1127                if ( o1 > dg[i]) { o1 = dg[i]; p1 = i; }
1128                if ( oe <= df[i]) { oe = df[i]; pe = i; }
1129            }
1130        }
1131        i--;
1132    }
1133}
1134
1135/*
1136*  make some changes of variables, see optvalues
1137*/
1138static void
1139cf_prepgcd( const CanonicalForm & f, const CanonicalForm & g, int & cc, int & p1, int &pe )
1140{
1141    int i, k, n;
1142    n = f.level();
1143    cc = 0;
1144    p1 = pe = n;
1145    if ( n == 1 )
1146        return;
1147    int * degsf = new int[n+1];
1148    int * degsg = new int[n+1];
1149    for ( i = n; i > 0; i-- )
1150    {
1151        degsf[i] = degsg[i] = 0;
1152    }
1153    degsf = degrees( f, degsf );
1154    degsg = degrees( g, degsg );
1155
1156    k = 0;
1157    for ( i = n-1; i > 0; i-- )
1158    {
[f4b180]1159        if ( degsf[i] == 0 )
[a7ec94]1160        {
1161            if ( degsg[i] != 0 )
1162            {
1163                cc = -i;
1164                break;
1165            }
1166        }
1167        else
1168        {
1169            if ( degsg[i] == 0 )
1170            {
1171                cc = i;
1172                break;
1173            }
1174            else k++;
1175        }
1176    }
1177
1178    if ( ( cc == 0 ) && ( k != 0 ) )
1179        optvalues( degsf, degsg, n, p1, pe );
1180    if ( ( pe != 1 ) && ( degsf[1] != 0 ) )
1181        pe = -pe;
[f4b180]1182
[a7ec94]1183    delete [] degsf;
1184    delete [] degsg;
1185}
[6f62c3]1186
1187
1188static CanonicalForm
1189balance_p ( const CanonicalForm & f, const CanonicalForm & q )
1190{
1191    Variable x = f.mvar();
1192    CanonicalForm result = 0, qh = q / 2;
1193    CanonicalForm c;
1194    CFIterator i;
1195    for ( i = f; i.hasTerms(); i++ )
1196    {
1197        c = i.coeff();
1198        if ( c.inCoeffDomain())
1199        {
1200          if ( c > qh )
1201            result += power( x, i.exp() ) * (c - q);
1202          else
1203            result += power( x, i.exp() ) * c;
1204        }
[f4b180]1205        else
[6f62c3]1206          result += power( x, i.exp() ) * balance_p(c,q);
1207    }
1208    return result;
1209}
1210
1211CanonicalForm chinrem_gcd ( const CanonicalForm & FF, const CanonicalForm & GG )
1212{
[297e92]1213  CanonicalForm f, g, cl, q(0), Dp, newD, D, newq;
[c1b9927]1214  int p, i, dp_deg, d_deg=-1;
[6f62c3]1215
[01e8874]1216  CanonicalForm cd ( bCommonDen( FF ));
[6f62c3]1217  f=cd*FF;
[297e92]1218  Variable x= Variable (1);
1219  CanonicalForm cf, cg;
1220  cf= icontent (f);
1221  f /= cf;
[08a6ebb]1222  //cd = bCommonDen( f ); f *=cd;
1223  //f /=vcontent(f,Variable(1));
[6f62c3]1224
1225  cd = bCommonDen( GG );
1226  g=cd*GG;
[297e92]1227  cg= icontent (g);
1228  g /= cg;
[08a6ebb]1229  //cd = bCommonDen( g ); g *=cd;
1230  //g /=vcontent(g,Variable(1));
[6f62c3]1231
[297e92]1232  CanonicalForm Dn, test= 0;
[4704674]1233  cl =  gcd (f.lc(),g.lc());
[597783]1234  CanonicalForm gcdcfcg= gcd (cf, cg);
[2488dc3]1235  CanonicalForm fp, gp;
[4704674]1236  CanonicalForm b= 1;
1237  int minCommonDeg= 0;
1238  for (i= tmax (f.level(), g.level()); i > 0; i--)
1239  {
1240    if (degree (f, i) <= 0 || degree (g, i) <= 0)
1241      continue;
1242    else
1243    {
1244      minCommonDeg= tmin (degree (g, i), degree (f, i));
1245      break;
1246    }
1247  }
1248  if (i == 0)
1249    return gcdcfcg;
1250  for (; i > 0; i--)
1251  {
1252    if (degree (f, i) <= 0 || degree (g, i) <= 0)
1253      continue;
1254    else
1255      minCommonDeg= tmin (minCommonDeg, tmin (degree (g, i), degree (f, i)));
1256  }
[cb7827]1257  b= 2*tmin (maxNorm (f), maxNorm (g))*abs (cl)*
1258     power (CanonicalForm (2), minCommonDeg);
[297e92]1259  bool equal= false;
[6f62c3]1260  i = cf_getNumBigPrimes() - 1;
1261
[cb7827]1262  CanonicalForm cof, cog, cofp, cogp, newCof, newCog, cofn, cogn;
[2488dc3]1263  int maxNumVars= tmax (getNumVars (f), getNumVars (g));
[297e92]1264  //Off (SW_RATIONAL);
[6f62c3]1265  while ( true )
1266  {
1267    p = cf_getBigPrime( i );
1268    i--;
[597783]1269    while ( i >= 0 && mod( cl*(lc(f)/cl)*(lc(g)/cl), p ) == 0 )
[6f62c3]1270    {
1271      p = cf_getBigPrime( i );
1272      i--;
1273    }
[c30347]1274    //printf("try p=%d\n",p);
[6f62c3]1275    setCharacteristic( p );
[2488dc3]1276    fp= mapinto (f);
1277    gp= mapinto (g);
[517530]1278#ifdef HAVE_NTL
[2488dc3]1279    if (size (fp)/maxNumVars > 500 && size (gp)/maxNumVars > 500)
1280      Dp = GCD_small_p (fp, gp, cofp, cogp);
1281    else
1282    {
1283      Dp= gcd_poly (fp, gp);
1284      cofp= fp/Dp;
1285      cogp= gp/Dp;
1286    }
[517530]1287#else
[2488dc3]1288    Dp= gcd_poly (fp, gp);
1289    cofp= fp/Dp;
1290    cogp= gp/Dp;
[517530]1291#endif
[08a6ebb]1292    Dp /=Dp.lc();
[cb7827]1293    cofp /= lc (cofp);
1294    cogp /= lc (cogp);
[6f62c3]1295    setCharacteristic( 0 );
1296    dp_deg=totaldegree(Dp);
1297    if ( dp_deg == 0 )
[c30347]1298    {
1299      //printf(" -> 1\n");
[297e92]1300      return CanonicalForm(gcdcfcg);
[c30347]1301    }
[6f62c3]1302    if ( q.isZero() )
1303    {
1304      D = mapinto( Dp );
[cb7827]1305      cof= mapinto (cofp);
1306      cog= mapinto (cogp);
[6f62c3]1307      d_deg=dp_deg;
1308      q = p;
1309    }
1310    else
1311    {
1312      if ( dp_deg == d_deg )
1313      {
1314        chineseRemainder( D, q, mapinto( Dp ), p, newD, newq );
[cb7827]1315        chineseRemainder( cof, q, mapinto (cofp), p, newCof, newq);
1316        chineseRemainder( cog, q, mapinto (cogp), p, newCog, newq);
1317        cof= newCof;
1318        cog= newCog;
[6f62c3]1319        q = newq;
1320        D = newD;
1321      }
[f4b180]1322      else if ( dp_deg < d_deg )
[6f62c3]1323      {
[c30347]1324        //printf(" were all bad, try more\n");
[6f62c3]1325        // all previous p's are bad primes
1326        q = p;
1327        D = mapinto( Dp );
[cb7827]1328        cof= mapinto (cof);
1329        cog= mapinto (cog);
[6f62c3]1330        d_deg=dp_deg;
[297e92]1331        test= 0;
1332        equal= false;
[6f62c3]1333      }
[c30347]1334      else
1335      {
1336        //printf(" was bad, try more\n");
1337      }
[f4b180]1338      //else dp_deg > d_deg: bad prime
[6f62c3]1339    }
[08a6ebb]1340    if ( i >= 0 )
[6f62c3]1341    {
[297e92]1342      Dn= Farey(D,q);
[cb7827]1343      cofn= Farey(cof,q);
1344      cogn= Farey(cog,q);
[297e92]1345      int is_rat= isOn (SW_RATIONAL);
[56d3c6]1346      On (SW_RATIONAL);
[297e92]1347      cd = bCommonDen( Dn ); // we need On(SW_RATIONAL)
[cb7827]1348      cofn *= bCommonDen (cofn);
1349      cogn *= bCommonDen (cogn);
[297e92]1350      if (!is_rat)
1351        Off (SW_RATIONAL);
[c992ec1]1352      Dn *=cd;
[297e92]1353      if (test != Dn)
1354        test= Dn;
1355      else
1356        equal= true;
[c992ec1]1357      //Dn /=vcontent(Dn,Variable(1));
[1e4b53]1358      if ((terminationTest (f,g, cofn, cogn, Dn)) ||
[597783]1359          ((equal || q > b) && fdivides (Dn, f) && fdivides (Dn, g)))
[6f62c3]1360      {
[c30347]1361        //printf(" -> success\n");
[297e92]1362        return Dn*gcdcfcg;
[6f62c3]1363      }
[297e92]1364      equal= false;
[c992ec1]1365      //else: try more primes
[6f62c3]1366    }
1367    else
[c992ec1]1368    { // try other method
[c30347]1369      //printf("try other gcd\n");
[6f62c3]1370      Off(SW_USE_CHINREM_GCD);
1371      D=gcd_poly( f, g );
1372      On(SW_USE_CHINREM_GCD);
[297e92]1373      return D*gcdcfcg;
[6f62c3]1374    }
1375  }
1376}
[c5d0aed]1377
[cb7827]1378
Note: See TracBrowser for help on using the repository browser.