source: git/factory/cf_gcd.cc @ 12f992

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