source: git/factory/cf_gcd.cc @ 2abbc76

spielwiese
Last change on this file since 2abbc76 was bffe62d, checked in by Martin Lee <martinlee84@…>, 11 years ago
chg: skip zz_p::init() if it is already correctly initialized
  • Property mode set to 100644
File size: 34.7 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    CanonicalForm pi, pi1;
610    CanonicalForm C, Ci, Ci1, Hi, bi, pi2;
611    bool bpure;
612    int delta = degree( f ) - degree( g );
613
614    if ( delta >= 0 )
615    {
616        pi = f; pi1 = g;
617    }
618    else
619    {
620        pi = g; pi1 = f; delta = -delta;
621    }
622    Ci = content( pi ); Ci1 = content( pi1 );
623    pi1 = pi1 / Ci1; pi = pi / Ci;
624    C = gcd( Ci, Ci1 );
625    int d= 0;
626    if ( !( pi.isUnivariate() && pi1.isUnivariate() ) )
627    {
628        if ( gcd_test_one( pi1, pi, true, d ) )
629        {
630          C=abs(C);
631          //out_cf("GCD:",C,"\n");
632          return C;
633        }
634        bpure = false;
635    }
636    else
637    {
638        bpure = isPurePoly(pi) && isPurePoly(pi1);
639#ifdef HAVE_FLINT
640        if (bpure && (CFFactory::gettype() != GaloisFieldDomain))
641          return gcd_univar_flintp(pi,pi1)*C;
642#else
643#ifdef HAVE_NTL
644        if ( isOn(SW_USE_NTL_GCD_P) && bpure && (CFFactory::gettype() != GaloisFieldDomain))
645            return gcd_univar_ntlp(pi, pi1 ) * C;
646#endif
647#endif
648    }
649    Variable v = f.mvar();
650    Hi = power( LC( pi1, v ), delta );
651    int maxNumVars= tmax (getNumVars (pi), getNumVars (pi1));
652
653    if (!(pi.isUnivariate() && pi1.isUnivariate()))
654    {
655      if (size (Hi)*size (pi)/(maxNumVars*3) > 500) //maybe this needs more tuning
656      {
657        On (SW_USE_FF_MOD_GCD);
658        C *= gcd (pi, pi1);
659        Off (SW_USE_FF_MOD_GCD);
660        return C;
661      }
662    }
663
664    if ( (delta+1) % 2 )
665        bi = 1;
666    else
667        bi = -1;
668    CanonicalForm oldPi= pi, oldPi1= pi1, powHi;
669    while ( degree( pi1, v ) > 0 )
670    {
671        if (!(pi.isUnivariate() && pi1.isUnivariate()))
672        {
673          if (size (pi)/maxNumVars > 500 || size (pi1)/maxNumVars > 500)
674          {
675            On (SW_USE_FF_MOD_GCD);
676            C *= gcd (oldPi, oldPi1);
677            Off (SW_USE_FF_MOD_GCD);
678            return C;
679          }
680        }
681        pi2 = psr( pi, pi1, v );
682        pi2 = pi2 / bi;
683        pi = pi1; pi1 = pi2;
684        maxNumVars= tmax (getNumVars (pi), getNumVars (pi1));
685        if (!(pi1.isUnivariate()) && (size (pi1)/maxNumVars > 500))
686        {
687            On (SW_USE_FF_MOD_GCD);
688            C *= gcd (oldPi, oldPi1);
689            Off (SW_USE_FF_MOD_GCD);
690            return C;
691        }
692        if ( degree( pi1, v ) > 0 )
693        {
694            delta = degree( pi, v ) - degree( pi1, v );
695            powHi= power (Hi, delta-1);
696            if ( (delta+1) % 2 )
697                bi = LC( pi, v ) * powHi*Hi;
698            else
699                bi = -LC( pi, v ) * powHi*Hi;
700            Hi = power( LC( pi1, v ), delta ) / powHi;
701            if (!(pi.isUnivariate() && pi1.isUnivariate()))
702            {
703              if (size (Hi)*size (pi)/(maxNumVars*3) > 500) //maybe this needs more tuning
704              {
705                On (SW_USE_FF_MOD_GCD);
706                C *= gcd (oldPi, oldPi1);
707                Off (SW_USE_FF_MOD_GCD);
708                return C;
709              }
710            }
711        }
712    }
713    if ( degree( pi1, v ) == 0 )
714    {
715      C=abs(C);
716      //out_cf("GCD:",C,"\n");
717      return C;
718    }
719    pi /= content( pi );
720    if ( bpure )
721        pi /= pi.lc();
722    C=abs(C*pi);
723    //out_cf("GCD:",C,"\n");
724    return C;
725}
726
727static CanonicalForm
728gcd_poly_0( const CanonicalForm & f, const CanonicalForm & g )
729{
730    CanonicalForm pi, pi1;
731    CanonicalForm C, Ci, Ci1, Hi, bi, pi2;
732    int delta = degree( f ) - degree( g );
733
734    if ( delta >= 0 )
735    {
736        pi = f; pi1 = g;
737    }
738    else
739    {
740        pi = g; pi1 = f; delta = -delta;
741    }
742    Ci = content( pi ); Ci1 = content( pi1 );
743    pi1 = pi1 / Ci1; pi = pi / Ci;
744    C = gcd( Ci, Ci1 );
745    int d= 0;
746    if ( pi.isUnivariate() && pi1.isUnivariate() )
747    {
748#ifdef HAVE_FLINT
749        if (isPurePoly(pi) && isPurePoly(pi1) )
750            return gcd_univar_flint0(pi, pi1 ) * C;
751#else
752#ifdef HAVE_NTL
753        if ( isOn(SW_USE_NTL_GCD_0) && isPurePoly(pi) && isPurePoly(pi1) )
754            return gcd_univar_ntl0(pi, pi1 ) * C;
755#endif
756#endif
757        return gcd_poly_univar0( pi, pi1, true ) * C;
758    }
759    else if ( gcd_test_one( pi1, pi, true, d ) )
760      return C;
761    Variable v = f.mvar();
762    Hi = power( LC( pi1, v ), delta );
763    if ( (delta+1) % 2 )
764        bi = 1;
765    else
766        bi = -1;
767    while ( degree( pi1, v ) > 0 )
768    {
769        pi2 = psr( pi, pi1, v );
770        pi2 = pi2 / bi;
771        pi = pi1; pi1 = pi2;
772        if ( degree( pi1, v ) > 0 )
773        {
774            delta = degree( pi, v ) - degree( pi1, v );
775            if ( (delta+1) % 2 )
776                bi = LC( pi, v ) * power( Hi, delta );
777            else
778                bi = -LC( pi, v ) * power( Hi, delta );
779            Hi = power( LC( pi1, v ), delta ) / power( Hi, delta-1 );
780        }
781    }
782    if ( degree( pi1, v ) == 0 )
783        return C;
784    else
785        return C * pp( pi );
786}
787
788//{{{ CanonicalForm gcd_poly ( const CanonicalForm & f, const CanonicalForm & g )
789//{{{ docu
790//
791// gcd_poly() - calculate polynomial gcd.
792//
793// This is the dispatcher for polynomial gcd calculation.  We call either
794// ezgcd(), sparsemod() or gcd_poly1() in dependecy on the current
795// characteristic and settings of SW_USE_EZGCD.
796//
797// Used by gcd() and gcd_poly_univar0().
798//
799//}}}
800#if 0
801int si_factor_reminder=1;
802#endif
803CanonicalForm gcd_poly ( const CanonicalForm & f, const CanonicalForm & g )
804{
805  CanonicalForm fc, gc, d1;
806  bool fc_isUnivariate=f.isUnivariate();
807  bool gc_isUnivariate=g.isUnivariate();
808  bool fc_and_gc_Univariate=fc_isUnivariate && gc_isUnivariate;
809  fc = f;
810  gc = g;
811  if ( getCharacteristic() != 0 )
812  {
813    #ifdef HAVE_NTL
814    if ((!fc_and_gc_Univariate) && (isOn( SW_USE_EZGCD_P )))
815    {
816      fc= EZGCD_P (fc, gc);
817    }
818    else if (isOn(SW_USE_FF_MOD_GCD) && !fc_and_gc_Univariate)
819    {
820      Variable a;
821      if (hasFirstAlgVar (fc, a) || hasFirstAlgVar (gc, a))
822        fc=GCD_Fp_extension (fc, gc, a);
823      else if (CFFactory::gettype() == GaloisFieldDomain)
824        fc=GCD_GF (fc, gc);
825      else
826        fc=GCD_small_p (fc, gc);
827    }
828    else
829    #endif
830    fc = gcd_poly_p( fc, gc );
831  }
832  else if (!fc_and_gc_Univariate)
833  {
834    if ( isOn( SW_USE_EZGCD ) )
835      fc= ezgcd (fc, gc);
836    else if (isOn(SW_USE_CHINREM_GCD))
837      fc = chinrem_gcd( fc, gc);
838    else
839    {
840       fc = gcd_poly_0( fc, gc );
841    }
842  }
843  else
844  {
845    fc = gcd_poly_0( fc, gc );
846  }
847  if ( d1.degree() > 0 )
848    fc *= d1;
849  return fc;
850}
851//}}}
852
853//{{{ static CanonicalForm cf_content ( const CanonicalForm & f, const CanonicalForm & g )
854//{{{ docu
855//
856// cf_content() - return gcd(g, content(f)).
857//
858// content(f) is calculated with respect to f's main variable.
859//
860// Used by gcd(), content(), content( CF, Variable ).
861//
862//}}}
863static CanonicalForm
864cf_content ( const CanonicalForm & f, const CanonicalForm & g )
865{
866    if ( f.inPolyDomain() || ( f.inExtension() && ! getReduce( f.mvar() ) ) )
867    {
868        CFIterator i = f;
869        CanonicalForm result = g;
870        while ( i.hasTerms() && ! result.isOne() )
871        {
872            result = gcd( i.coeff(), result );
873            i++;
874        }
875        return result;
876    }
877    else
878        return abs( f );
879}
880//}}}
881
882//{{{ CanonicalForm content ( const CanonicalForm & f )
883//{{{ docu
884//
885// content() - return content(f) with respect to main variable.
886//
887// Normalizes result.
888//
889//}}}
890CanonicalForm
891content ( const CanonicalForm & f )
892{
893    if ( f.inPolyDomain() || ( f.inExtension() && ! getReduce( f.mvar() ) ) )
894    {
895        CFIterator i = f;
896        CanonicalForm result = abs( i.coeff() );
897        i++;
898        while ( i.hasTerms() && ! result.isOne() )
899        {
900            result = gcd( i.coeff(), result );
901            i++;
902        }
903        return result;
904    }
905    else
906        return abs( f );
907}
908//}}}
909
910//{{{ CanonicalForm content ( const CanonicalForm & f, const Variable & x )
911//{{{ docu
912//
913// content() - return content(f) with respect to x.
914//
915// x should be a polynomial variable.
916//
917//}}}
918CanonicalForm
919content ( const CanonicalForm & f, const Variable & x )
920{
921    ASSERT( x.level() > 0, "cannot calculate content with respect to algebraic variable" );
922    Variable y = f.mvar();
923
924    if ( y == x )
925        return cf_content( f, 0 );
926    else  if ( y < x )
927        return f;
928    else
929        return swapvar( content( swapvar( f, y, x ), y ), y, x );
930}
931//}}}
932
933//{{{ CanonicalForm vcontent ( const CanonicalForm & f, const Variable & x )
934//{{{ docu
935//
936// vcontent() - return content of f with repect to variables >= x.
937//
938// The content is recursively calculated over all coefficients in
939// f having level less than x.  x should be a polynomial
940// variable.
941//
942//}}}
943CanonicalForm
944vcontent ( const CanonicalForm & f, const Variable & x )
945{
946    ASSERT( x.level() > 0, "cannot calculate vcontent with respect to algebraic variable" );
947
948    if ( f.mvar() <= x )
949        return content( f, x );
950    else {
951        CFIterator i;
952        CanonicalForm d = 0;
953        for ( i = f; i.hasTerms() && ! d.isOne(); i++ )
954            d = gcd( d, vcontent( i.coeff(), x ) );
955        return d;
956    }
957}
958//}}}
959
960//{{{ CanonicalForm pp ( const CanonicalForm & f )
961//{{{ docu
962//
963// pp() - return primitive part of f.
964//
965// Returns zero if f equals zero, otherwise f / content(f).
966//
967//}}}
968CanonicalForm
969pp ( const CanonicalForm & f )
970{
971    if ( f.isZero() )
972        return f;
973    else
974        return f / content( f );
975}
976//}}}
977
978CanonicalForm
979gcd ( const CanonicalForm & f, const CanonicalForm & g )
980{
981    bool b = f.isZero();
982    if ( b || g.isZero() )
983    {
984        if ( b )
985            return abs( g );
986        else
987            return abs( f );
988    }
989    if ( f.inPolyDomain() || g.inPolyDomain() )
990    {
991        if ( f.mvar() != g.mvar() )
992        {
993            if ( f.mvar() > g.mvar() )
994                return cf_content( f, g );
995            else
996                return cf_content( g, f );
997        }
998        if (isOn(SW_USE_QGCD))
999        {
1000          Variable m;
1001          if (
1002          (getCharacteristic() == 0) &&
1003          (hasFirstAlgVar(f,m) || hasFirstAlgVar(g,m))
1004          )
1005          {
1006            bool on_rational = isOn(SW_RATIONAL);
1007            CanonicalForm r=QGCD(f,g);
1008            On(SW_RATIONAL);
1009            CanonicalForm cdF = bCommonDen( r );
1010            if (!on_rational) Off(SW_RATIONAL);
1011            return cdF*r;
1012          }
1013        }
1014
1015        if ( f.inExtension() && getReduce( f.mvar() ) )
1016            return CanonicalForm(1);
1017        else
1018        {
1019            if ( fdivides( f, g ) )
1020                return abs( f );
1021            else  if ( fdivides( g, f ) )
1022                return abs( g );
1023            if ( !( getCharacteristic() == 0 && isOn( SW_RATIONAL ) ) )
1024            {
1025                CanonicalForm d;
1026                d = gcd_poly( f, g );
1027                return abs( d );
1028            }
1029            else
1030            {
1031                //printf ("here\n");
1032                CanonicalForm cdF = bCommonDen( f );
1033                CanonicalForm cdG = bCommonDen( g );
1034                Off( SW_RATIONAL );
1035                CanonicalForm l = lcm( cdF, cdG );
1036                On( SW_RATIONAL );
1037                CanonicalForm F = f * l, G = g * l;
1038                Off( SW_RATIONAL );
1039                l = gcd_poly( F, G );
1040                On( SW_RATIONAL );
1041                return abs( l );
1042            }
1043        }
1044    }
1045    if ( f.inBaseDomain() && g.inBaseDomain() )
1046        return bgcd( f, g );
1047    else
1048        return 1;
1049}
1050
1051//{{{ CanonicalForm lcm ( const CanonicalForm & f, const CanonicalForm & g )
1052//{{{ docu
1053//
1054// lcm() - return least common multiple of f and g.
1055//
1056// The lcm is calculated using the formula lcm(f, g) = f * g / gcd(f, g).
1057//
1058// Returns zero if one of f or g equals zero.
1059//
1060//}}}
1061CanonicalForm
1062lcm ( const CanonicalForm & f, const CanonicalForm & g )
1063{
1064    if ( f.isZero() || g.isZero() )
1065        return 0;
1066    else
1067        return ( f / gcd( f, g ) ) * g;
1068}
1069//}}}
1070
1071#ifdef HAVE_NTL
1072#ifndef HAVE_FLINT
1073static CanonicalForm
1074gcd_univar_ntl0( const CanonicalForm & F, const CanonicalForm & G )
1075{
1076    ZZX F1=convertFacCF2NTLZZX(F);
1077    ZZX G1=convertFacCF2NTLZZX(G);
1078    ZZX R=GCD(F1,G1);
1079    return convertNTLZZX2CF(R,F.mvar());
1080}
1081
1082static CanonicalForm
1083gcd_univar_ntlp( const CanonicalForm & F, const CanonicalForm & G )
1084{
1085  if (fac_NTL_char!=getCharacteristic())
1086  {
1087    fac_NTL_char=getCharacteristic();
1088    #ifdef NTL_ZZ
1089    ZZ r;
1090    r=getCharacteristic();
1091    ZZ_pContext ccc(r);
1092    #else
1093    zz_pContext ccc(getCharacteristic());
1094    #endif
1095    ccc.restore();
1096    #ifdef NTL_ZZ
1097    ZZ_p::init(r);
1098    #else
1099    zz_p::init(getCharacteristic());
1100    #endif
1101  }
1102  #ifdef NTL_ZZ
1103  ZZ_pX F1=convertFacCF2NTLZZpX(F);
1104  ZZ_pX G1=convertFacCF2NTLZZpX(G);
1105  ZZ_pX R=GCD(F1,G1);
1106  return  convertNTLZZpX2CF(R,F.mvar());
1107  #else
1108  zz_pX F1=convertFacCF2NTLzzpX(F);
1109  zz_pX G1=convertFacCF2NTLzzpX(G);
1110  zz_pX R=GCD(F1,G1);
1111  return  convertNTLzzpX2CF(R,F.mvar());
1112  #endif
1113}
1114#endif
1115#endif
1116
1117#ifdef HAVE_FLINT
1118static CanonicalForm
1119gcd_univar_flintp (const CanonicalForm& F, const CanonicalForm& G)
1120{
1121  nmod_poly_t F1, G1;
1122  convertFacCF2nmod_poly_t (F1, F);
1123  convertFacCF2nmod_poly_t (G1, G);
1124  nmod_poly_gcd (F1, F1, G1);
1125  CanonicalForm result= convertnmod_poly_t2FacCF (F1, F.mvar());
1126  nmod_poly_clear (F1);
1127  nmod_poly_clear (G1);
1128  return result;
1129}
1130
1131static CanonicalForm
1132gcd_univar_flint0( const CanonicalForm & F, const CanonicalForm & G )
1133{
1134  fmpz_poly_t F1, G1;
1135  convertFacCF2Fmpz_poly_t(F1, F);
1136  convertFacCF2Fmpz_poly_t(G1, G);
1137  fmpz_poly_gcd (F1, F1, G1);
1138  CanonicalForm result= convertFmpz_poly_t2FacCF (F1, F.mvar());
1139  fmpz_poly_clear (F1);
1140  fmpz_poly_clear (G1);
1141  return result;
1142}
1143#endif
1144
1145
1146/*
1147*  compute positions p1 and pe of optimal variables:
1148*    pe is used in "ezgcd" and
1149*    p1 in "gcd_poly1"
1150*/
1151/*static
1152void optvalues ( const int * df, const int * dg, const int n, int & p1, int &pe )
1153{
1154    int i, o1, oe;
1155    if ( df[n] > dg[n] )
1156    {
1157        o1 = df[n]; oe = dg[n];
1158    }
1159    else
1160    {
1161        o1 = dg[n]; oe = df[n];
1162    }
1163    i = n-1;
1164    while ( i > 0 )
1165    {
1166        if ( df[i] != 0 )
1167        {
1168            if ( df[i] > dg[i] )
1169            {
1170                if ( o1 > df[i]) { o1 = df[i]; p1 = i; }
1171                if ( oe <= dg[i]) { oe = dg[i]; pe = i; }
1172            }
1173            else
1174            {
1175                if ( o1 > dg[i]) { o1 = dg[i]; p1 = i; }
1176                if ( oe <= df[i]) { oe = df[i]; pe = i; }
1177            }
1178        }
1179        i--;
1180    }
1181}*/
1182
1183/*
1184*  make some changes of variables, see optvalues
1185*/
1186/*static void
1187cf_prepgcd( const CanonicalForm & f, const CanonicalForm & g, int & cc, int & p1, int &pe )
1188{
1189    int i, k, n;
1190    n = f.level();
1191    cc = 0;
1192    p1 = pe = n;
1193    if ( n == 1 )
1194        return;
1195    int * degsf = new int[n+1];
1196    int * degsg = new int[n+1];
1197    for ( i = n; i > 0; i-- )
1198    {
1199        degsf[i] = degsg[i] = 0;
1200    }
1201    degsf = degrees( f, degsf );
1202    degsg = degrees( g, degsg );
1203
1204    k = 0;
1205    for ( i = n-1; i > 0; i-- )
1206    {
1207        if ( degsf[i] == 0 )
1208        {
1209            if ( degsg[i] != 0 )
1210            {
1211                cc = -i;
1212                break;
1213            }
1214        }
1215        else
1216        {
1217            if ( degsg[i] == 0 )
1218            {
1219                cc = i;
1220                break;
1221            }
1222            else k++;
1223        }
1224    }
1225
1226    if ( ( cc == 0 ) && ( k != 0 ) )
1227        optvalues( degsf, degsg, n, p1, pe );
1228    if ( ( pe != 1 ) && ( degsf[1] != 0 ) )
1229        pe = -pe;
1230
1231    delete [] degsf;
1232    delete [] degsg;
1233}*/
1234
1235TIMING_DEFINE_PRINT(chinrem_termination)
1236TIMING_DEFINE_PRINT(chinrem_recursion)
1237TIMING_DEFINE_PRINT(chinrem_reconstruction)
1238
1239CanonicalForm chinrem_gcd ( const CanonicalForm & FF, const CanonicalForm & GG )
1240{
1241  CanonicalForm f, g, cl, q(0), Dp, newD, D, newq;
1242  int p, i, dp_deg, d_deg=-1;
1243
1244  CanonicalForm cd ( bCommonDen( FF ));
1245  f=cd*FF;
1246  Variable x= Variable (1);
1247  CanonicalForm cf, cg;
1248  cf= icontent (f);
1249  f /= cf;
1250  //cd = bCommonDen( f ); f *=cd;
1251  //f /=vcontent(f,Variable(1));
1252
1253  cd = bCommonDen( GG );
1254  g=cd*GG;
1255  cg= icontent (g);
1256  g /= cg;
1257  //cd = bCommonDen( g ); g *=cd;
1258  //g /=vcontent(g,Variable(1));
1259
1260  CanonicalForm Dn, test= 0;
1261  cl =  gcd (f.lc(),g.lc());
1262  CanonicalForm gcdcfcg= gcd (cf, cg);
1263  CanonicalForm fp, gp;
1264  CanonicalForm b= 1;
1265  int minCommonDeg= 0;
1266  for (i= tmax (f.level(), g.level()); i > 0; i--)
1267  {
1268    if (degree (f, i) <= 0 || degree (g, i) <= 0)
1269      continue;
1270    else
1271    {
1272      minCommonDeg= tmin (degree (g, i), degree (f, i));
1273      break;
1274    }
1275  }
1276  if (i == 0)
1277    return gcdcfcg;
1278  for (; i > 0; i--)
1279  {
1280    if (degree (f, i) <= 0 || degree (g, i) <= 0)
1281      continue;
1282    else
1283      minCommonDeg= tmin (minCommonDeg, tmin (degree (g, i), degree (f, i)));
1284  }
1285  b= 2*tmin (maxNorm (f), maxNorm (g))*abs (cl)*
1286     power (CanonicalForm (2), minCommonDeg);
1287  bool equal= false;
1288  i = cf_getNumBigPrimes() - 1;
1289
1290  CanonicalForm cof, cog, cofp, cogp, newCof, newCog, cofn, cogn;
1291  int maxNumVars= tmax (getNumVars (f), getNumVars (g));
1292  //Off (SW_RATIONAL);
1293  while ( true )
1294  {
1295    p = cf_getBigPrime( i );
1296    i--;
1297    while ( i >= 0 && mod( cl*(lc(f)/cl)*(lc(g)/cl), p ) == 0 )
1298    {
1299      p = cf_getBigPrime( i );
1300      i--;
1301    }
1302    //printf("try p=%d\n",p);
1303    setCharacteristic( p );
1304    fp= mapinto (f);
1305    gp= mapinto (g);
1306    TIMING_START (chinrem_recursion)
1307#ifdef HAVE_NTL
1308    if (size (fp)/maxNumVars > 500 && size (gp)/maxNumVars > 500)
1309      Dp = GCD_small_p (fp, gp, cofp, cogp);
1310    else
1311    {
1312      Dp= gcd_poly (fp, gp);
1313      cofp= fp/Dp;
1314      cogp= gp/Dp;
1315    }
1316#else
1317    Dp= gcd_poly (fp, gp);
1318    cofp= fp/Dp;
1319    cogp= gp/Dp;
1320#endif
1321    TIMING_END_AND_PRINT (chinrem_recursion,
1322                          "time for gcd mod p in modular gcd: ");
1323    Dp /=Dp.lc();
1324    cofp /= lc (cofp);
1325    cogp /= lc (cogp);
1326    setCharacteristic( 0 );
1327    dp_deg=totaldegree(Dp);
1328    if ( dp_deg == 0 )
1329    {
1330      //printf(" -> 1\n");
1331      return CanonicalForm(gcdcfcg);
1332    }
1333    if ( q.isZero() )
1334    {
1335      D = mapinto( Dp );
1336      cof= mapinto (cofp);
1337      cog= mapinto (cogp);
1338      d_deg=dp_deg;
1339      q = p;
1340    }
1341    else
1342    {
1343      if ( dp_deg == d_deg )
1344      {
1345        chineseRemainder( D, q, mapinto( Dp ), p, newD, newq );
1346        chineseRemainder( cof, q, mapinto (cofp), p, newCof, newq);
1347        chineseRemainder( cog, q, mapinto (cogp), p, newCog, newq);
1348        cof= newCof;
1349        cog= newCog;
1350        q = newq;
1351        D = newD;
1352      }
1353      else if ( dp_deg < d_deg )
1354      {
1355        //printf(" were all bad, try more\n");
1356        // all previous p's are bad primes
1357        q = p;
1358        D = mapinto( Dp );
1359        cof= mapinto (cof);
1360        cog= mapinto (cog);
1361        d_deg=dp_deg;
1362        test= 0;
1363        equal= false;
1364      }
1365      else
1366      {
1367        //printf(" was bad, try more\n");
1368      }
1369      //else dp_deg > d_deg: bad prime
1370    }
1371    if ( i >= 0 )
1372    {
1373      TIMING_START (chinrem_reconstruction);
1374      Dn= Farey(D,q);
1375      cofn= Farey(cof,q);
1376      cogn= Farey(cog,q);
1377      TIMING_END_AND_PRINT (chinrem_reconstruction,
1378                           "time for rational reconstruction in modular gcd: ");
1379      int is_rat= isOn (SW_RATIONAL);
1380      On (SW_RATIONAL);
1381      cd = bCommonDen( Dn ); // we need On(SW_RATIONAL)
1382      cofn *= bCommonDen (cofn);
1383      cogn *= bCommonDen (cogn);
1384      if (!is_rat)
1385        Off (SW_RATIONAL);
1386      Dn *=cd;
1387      if (test != Dn)
1388        test= Dn;
1389      else
1390        equal= true;
1391      //Dn /=vcontent(Dn,Variable(1));
1392      TIMING_START (chinrem_termination);
1393      if ((terminationTest (f,g, cofn, cogn, Dn)) ||
1394          ((equal || q > b) && fdivides (Dn, f) && fdivides (Dn, g)))
1395      {
1396        TIMING_END_AND_PRINT (chinrem_termination,
1397                            "time for successful termination in modular gcd: ");
1398        //printf(" -> success\n");
1399        return Dn*gcdcfcg;
1400      }
1401      TIMING_END_AND_PRINT (chinrem_termination,
1402                          "time for unsuccessful termination in modular gcd: ");
1403      equal= false;
1404      //else: try more primes
1405    }
1406    else
1407    { // try other method
1408      //printf("try other gcd\n");
1409      Off(SW_USE_CHINREM_GCD);
1410      D=gcd_poly( f, g );
1411      On(SW_USE_CHINREM_GCD);
1412      return D*gcdcfcg;
1413    }
1414  }
1415}
1416
1417
Note: See TracBrowser for help on using the repository browser.