source: git/factory/cf_gcd.cc @ b52d27

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