source: git/factory/cf_gcd.cc @ 80b8fe

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