source: git/factory/cf_gcd.cc @ e2e4be

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