source: git/factory/cf_gcd.cc @ e4e36c

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