source: git/factory/cf_gcd.cc @ ba5e9e

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