source: git/factory/cf_gcd.cc @ 4704674

spielwiese
Last change on this file since 4704674 was 4704674, checked in by Martin Lee <martinlee84@…>, 11 years ago
chg: avoid divisibility tests chg: use a bound on coeffs in chinrem_gcd
  • Property mode set to 100644
File size: 32.4 KB
Line 
1/* emacs edit mode for this file is -*- C++ -*- */
2/* $Id$ */
3
4#include "config.h"
5
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 "fieldGCD.h"
20#include "cf_gcd_smallp.h"
21#include "cf_map_ext.h"
22#include "cf_util.h"
23#include "gfops.h"
24
25#ifdef HAVE_NTL
26#include <NTL/ZZX.h>
27#include "NTLconvert.h"
28bool isPurePoly(const CanonicalForm & );
29static CanonicalForm gcd_univar_ntl0( const CanonicalForm &, const CanonicalForm & );
30static CanonicalForm gcd_univar_ntlp( const CanonicalForm &, const CanonicalForm & );
31#endif
32
33#ifdef HAVE_FLINT
34#include "FLINTconvert.h"
35static CanonicalForm gcd_univar_flint0 (const CanonicalForm &, const CanonicalForm &);
36static CanonicalForm gcd_univar_flintp (const CanonicalForm &, const CanonicalForm &);
37#endif
38
39static CanonicalForm cf_content ( const CanonicalForm &, const CanonicalForm & );
40static void cf_prepgcd( const CanonicalForm &, const CanonicalForm &, int &, int &, int & );
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  if (isOn(SW_USE_NTL_GCD_P) && ( getCharacteristic() > 0 ) && (CFFactory::gettype() != GaloisFieldDomain)
312  &&  (f.level()==g.level()) && isPurePoly(f) && isPurePoly(g))
313  {
314    if (fac_NTL_char!=getCharacteristic())
315    {
316      fac_NTL_char=getCharacteristic();
317      #ifdef NTL_ZZ
318      ZZ r;
319      r=getCharacteristic();
320      ZZ_pContext ccc(r);
321      #else
322      zz_pContext ccc(getCharacteristic());
323      #endif
324      ccc.restore();
325      #ifdef NTL_ZZ
326      ZZ_p::init(r);
327      #else
328      zz_p::init(getCharacteristic());
329      #endif
330    }
331    #ifdef NTL_ZZ
332    ZZ_pX F1=convertFacCF2NTLZZpX(f);
333    ZZ_pX G1=convertFacCF2NTLZZpX(g);
334    ZZ_pX R;
335    ZZ_pX A,B;
336    XGCD(R,A,B,F1,G1);
337    a=convertNTLZZpX2CF(A,f.mvar());
338    b=convertNTLZZpX2CF(B,f.mvar());
339    return convertNTLZZpX2CF(R,f.mvar());
340    #else
341    zz_pX F1=convertFacCF2NTLzzpX(f);
342    zz_pX G1=convertFacCF2NTLzzpX(g);
343    zz_pX R;
344    zz_pX A,B;
345    XGCD(R,A,B,F1,G1);
346    a=convertNTLzzpX2CF(A,f.mvar());
347    b=convertNTLzzpX2CF(B,f.mvar());
348    return convertNTLzzpX2CF(R,f.mvar());
349    #endif
350  }
351  if (isOn(SW_USE_NTL_GCD_0) && ( getCharacteristic() ==0)
352  && (f.level()==g.level()) && isPurePoly(f) && isPurePoly(g))
353  {
354    CanonicalForm fc=bCommonDen(f);
355    CanonicalForm gc=bCommonDen(g);
356    ZZX F1=convertFacCF2NTLZZX(f*fc);
357    ZZX G1=convertFacCF2NTLZZX(g*gc);
358    ZZX R=GCD(F1,G1);
359    CanonicalForm r=convertNTLZZX2CF(R,f.mvar());
360    ZZ RR;
361    ZZX A,B;
362    if (r.inCoeffDomain())
363    {
364      XGCD(RR,A,B,F1,G1,1);
365      CanonicalForm rr=convertZZ2CF(RR);
366      ASSERT (!rr.isZero(), "NTL:XGCD failed");
367      a=convertNTLZZX2CF(A,f.mvar())*fc/rr;
368      b=convertNTLZZX2CF(B,f.mvar())*gc/rr;
369      return CanonicalForm(1);
370    }
371    else
372    {
373      fc=bCommonDen(f);
374      gc=bCommonDen(g);
375      F1=convertFacCF2NTLZZX(f*fc/r);
376      G1=convertFacCF2NTLZZX(g*gc/r);
377      XGCD(RR,A,B,F1,G1,1);
378      a=convertNTLZZX2CF(A,f.mvar())*fc;
379      b=convertNTLZZX2CF(B,f.mvar())*gc;
380      CanonicalForm rr=convertZZ2CF(RR);
381      ASSERT (!rr.isZero(), "NTL:XGCD failed");
382      r*=rr;
383      if ( r.sign() < 0 ) { r= -r; a= -a; b= -b; }
384      return r;
385    }
386  }
387#endif
388  // may contain bug in the co-factors, see track 107
389  CanonicalForm contf = content( f );
390  CanonicalForm contg = content( g );
391
392  CanonicalForm p0 = f / contf, p1 = g / contg;
393  CanonicalForm f0 = 1, f1 = 0, g0 = 0, g1 = 1, q, r;
394
395  while ( ! p1.isZero() )
396  {
397      divrem( p0, p1, q, r );
398      p0 = p1; p1 = r;
399      r = g0 - g1 * q;
400      g0 = g1; g1 = r;
401      r = f0 - f1 * q;
402      f0 = f1; f1 = r;
403  }
404  CanonicalForm contp0 = content( p0 );
405  a = f0 / ( contf * contp0 );
406  b = g0 / ( contg * contp0 );
407  p0 /= contp0;
408  if ( p0.sign() < 0 )
409  {
410      p0 = -p0;
411      a = -a;
412      b = -b;
413  }
414  return p0;
415}
416//}}}
417
418//{{{ static CanonicalForm balance ( const CanonicalForm & f, const CanonicalForm & q )
419//{{{ docu
420//
421// balance() - map f from positive to symmetric representation
422//   mod q.
423//
424// This makes sense for univariate polynomials over Z only.
425// q should be an integer.
426//
427// Used by gcd_poly_univar0().
428//
429//}}}
430static CanonicalForm
431balance ( const CanonicalForm & f, const CanonicalForm & q )
432{
433    Variable x = f.mvar();
434    CanonicalForm result = 0, qh = q / 2;
435    CanonicalForm c;
436    CFIterator i;
437    for ( i = f; i.hasTerms(); i++ ) {
438        c = mod( i.coeff(), q );
439        if ( c > qh )
440            result += power( x, i.exp() ) * (c - q);
441        else
442            result += power( x, i.exp() ) * c;
443    }
444    return result;
445}
446//}}}
447
448static CanonicalForm gcd_poly_univar0( const CanonicalForm & F, const CanonicalForm & G, bool primitive )
449{
450  CanonicalForm f, g, c, cg, cl, BB, B, M, q, Dp, newD, D, newq;
451  int p, i;
452
453  if ( primitive )
454  {
455    f = F;
456    g = G;
457    c = 1;
458  }
459  else
460  {
461    CanonicalForm cF = content( F ), cG = content( G );
462    f = F / cF;
463    g = G / cG;
464    c = bgcd( cF, cG );
465  }
466  cg = gcd( f.lc(), g.lc() );
467  cl = ( f.lc() / cg ) * g.lc();
468//     B = 2 * cg * tmin(
469//         maxnorm(f)*power(CanonicalForm(2),f.degree())*isqrt(f.degree()+1),
470//         maxnorm(g)*power(CanonicalForm(2),g.degree())*isqrt(g.degree()+1)
471//         )+1;
472  M = tmin( maxNorm(f), maxNorm(g) );
473  BB = power(CanonicalForm(2),tmin(f.degree(),g.degree()))*M;
474  q = 0;
475  i = cf_getNumSmallPrimes() - 1;
476  while ( true )
477  {
478    B = BB;
479    while ( i >= 0 && q < B )
480    {
481      p = cf_getSmallPrime( i );
482      i--;
483      while ( i >= 0 && mod( cl, p ) == 0 )
484      {
485        p = cf_getSmallPrime( i );
486        i--;
487      }
488      setCharacteristic( p );
489      Dp = gcd( mapinto( f ), mapinto( g ) );
490      Dp = ( Dp / Dp.lc() ) * mapinto( cg );
491      setCharacteristic( 0 );
492      if ( Dp.degree() == 0 )
493        return c;
494      if ( q.isZero() )
495      {
496        D = mapinto( Dp );
497        q = p;
498        B = power(CanonicalForm(2),D.degree())*M+1;
499      }
500      else
501      {
502        if ( Dp.degree() == D.degree() )
503        {
504          chineseRemainder( D, q, mapinto( Dp ), p, newD, newq );
505          q = newq;
506          D = newD;
507        }
508        else if ( Dp.degree() < D.degree() )
509        {
510          // all previous p's are bad primes
511          q = p;
512          D = mapinto( Dp );
513          B = power(CanonicalForm(2),D.degree())*M+1;
514        }
515        // else p is a bad prime
516      }
517    }
518    if ( i >= 0 )
519    {
520      // now balance D mod q
521      D = pp( balance( D, q ) );
522      if ( fdivides( D, f ) && fdivides( D, g ) )
523        return D * c;
524      else
525        q = 0;
526    }
527    else
528      return gcd_poly( F, G );
529    DEBOUTLN( cerr, "another try ..." );
530  }
531}
532
533static CanonicalForm
534gcd_poly_p( const CanonicalForm & f, const CanonicalForm & g )
535{
536    CanonicalForm pi, pi1;
537    CanonicalForm C, Ci, Ci1, Hi, bi, pi2;
538    bool bpure;
539    int delta = degree( f ) - degree( g );
540
541    if ( delta >= 0 )
542    {
543        pi = f; pi1 = g;
544    }
545    else
546    {
547        pi = g; pi1 = f; delta = -delta;
548    }
549    Ci = content( pi ); Ci1 = content( pi1 );
550    pi1 = pi1 / Ci1; pi = pi / Ci;
551    C = gcd( Ci, Ci1 );
552    if ( !( pi.isUnivariate() && pi1.isUnivariate() ) )
553    {
554        if ( gcd_test_one( pi1, pi, true ) )
555        {
556          C=abs(C);
557          //out_cf("GCD:",C,"\n");
558          return C;
559        }
560        bpure = false;
561    }
562    else
563    {
564        bpure = isPurePoly(pi) && isPurePoly(pi1);
565#ifdef HAVE_FLINT
566        if (bpure && (CFFactory::gettype() != GaloisFieldDomain))
567          return gcd_univar_flintp(pi,pi1)*C;
568#else
569#ifdef HAVE_NTL
570        if ( isOn(SW_USE_NTL_GCD_P) && bpure && (CFFactory::gettype() != GaloisFieldDomain))
571            return gcd_univar_ntlp(pi, pi1 ) * C;
572#endif
573#endif
574    }
575    Variable v = f.mvar();
576    Hi = power( LC( pi1, v ), delta );
577    if ( (delta+1) % 2 )
578        bi = 1;
579    else
580        bi = -1;
581    int maxNumVars= tmax (getNumVars (pi), getNumVars (pi1));
582    CanonicalForm oldPi= pi, oldPi1= pi1;
583    while ( degree( pi1, v ) > 0 )
584    {
585        if (!(pi.isUnivariate() && pi1.isUnivariate()))
586        {
587          if (size (pi)/maxNumVars > 500 || size (pi1)/maxNumVars > 500)
588          {
589            On (SW_USE_FF_MOD_GCD);
590            C *= gcd (oldPi, oldPi1);
591            Off (SW_USE_FF_MOD_GCD);
592            return C;
593          }
594        }
595        pi2 = psr( pi, pi1, v );
596        pi2 = pi2 / bi;
597        pi = pi1; pi1 = pi2;
598        maxNumVars= tmax (getNumVars (pi), getNumVars (pi1));
599        if ( degree( pi1, v ) > 0 )
600        {
601            delta = degree( pi, v ) - degree( pi1, v );
602            if ( (delta+1) % 2 )
603                bi = LC( pi, v ) * power( Hi, delta );
604            else
605                bi = -LC( pi, v ) * power( Hi, delta );
606            Hi = power( LC( pi1, v ), delta ) / power( Hi, delta-1 );
607        }
608    }
609    if ( degree( pi1, v ) == 0 )
610    {
611      C=abs(C);
612      //out_cf("GCD:",C,"\n");
613      return C;
614    }
615    pi /= content( pi );
616    if ( bpure )
617        pi /= pi.lc();
618    C=abs(C*pi);
619    //out_cf("GCD:",C,"\n");
620    return C;
621}
622
623static CanonicalForm
624gcd_poly_0( const CanonicalForm & f, const CanonicalForm & g )
625{
626    CanonicalForm pi, pi1;
627    CanonicalForm C, Ci, Ci1, Hi, bi, pi2;
628    int delta = degree( f ) - degree( g );
629
630    if ( delta >= 0 )
631    {
632        pi = f; pi1 = g;
633    }
634    else
635    {
636        pi = g; pi1 = f; delta = -delta;
637    }
638    Ci = content( pi ); Ci1 = content( pi1 );
639    pi1 = pi1 / Ci1; pi = pi / Ci;
640    C = gcd( Ci, Ci1 );
641    if ( pi.isUnivariate() && pi1.isUnivariate() )
642    {
643/*#ifdef HAVE_FLINT
644        if (isPurePoly(pi) && isPurePoly(pi1) )
645            return gcd_univar_flint0(pi, pi1 ) * C;
646#else*/
647#ifdef HAVE_NTL
648        if ( isOn(SW_USE_NTL_GCD_0) && isPurePoly(pi) && isPurePoly(pi1) )
649            return gcd_univar_ntl0(pi, pi1 ) * C;
650#endif
651//#endif
652        return gcd_poly_univar0( pi, pi1, true ) * C;
653    }
654    else if ( gcd_test_one( pi1, pi, true ) )
655      return C;
656    Variable v = f.mvar();
657    Hi = power( LC( pi1, v ), delta );
658    if ( (delta+1) % 2 )
659        bi = 1;
660    else
661        bi = -1;
662    while ( degree( pi1, v ) > 0 )
663    {
664        pi2 = psr( pi, pi1, v );
665        pi2 = pi2 / bi;
666        pi = pi1; pi1 = pi2;
667        if ( degree( pi1, v ) > 0 )
668        {
669            delta = degree( pi, v ) - degree( pi1, v );
670            if ( (delta+1) % 2 )
671                bi = LC( pi, v ) * power( Hi, delta );
672            else
673                bi = -LC( pi, v ) * power( Hi, delta );
674            Hi = power( LC( pi1, v ), delta ) / power( Hi, delta-1 );
675        }
676    }
677    if ( degree( pi1, v ) == 0 )
678        return C;
679    else
680        return C * pp( pi );
681}
682
683//{{{ CanonicalForm gcd_poly ( const CanonicalForm & f, const CanonicalForm & g )
684//{{{ docu
685//
686// gcd_poly() - calculate polynomial gcd.
687//
688// This is the dispatcher for polynomial gcd calculation.  We call either
689// ezgcd(), sparsemod() or gcd_poly1() in dependecy on the current
690// characteristic and settings of SW_USE_EZGCD and SW_USE_SPARSEMOD, resp.
691//
692// Used by gcd() and gcd_poly_univar0().
693//
694//}}}
695#if 0
696int si_factor_reminder=1;
697#endif
698CanonicalForm gcd_poly ( const CanonicalForm & f, const CanonicalForm & g )
699{
700  CanonicalForm fc, gc, d1;
701  int mp, cc, p1, pe;
702  mp = f.level()+1;
703  bool fc_isUnivariate=f.isUnivariate();
704  bool gc_isUnivariate=g.isUnivariate();
705  bool fc_and_gc_Univariate=fc_isUnivariate && gc_isUnivariate;
706  cf_prepgcd( f, g, cc, p1, pe);
707  if ( cc != 0 )
708  {
709    if ( cc > 0 )
710    {
711      fc = replacevar( f, Variable(cc), Variable(mp) );
712      gc = g;
713    }
714    else
715    {
716      fc = replacevar( g, Variable(-cc), Variable(mp) );
717      gc = f;
718    }
719    return cf_content( fc, gc );
720  }
721// now each appearing variable is in f and g
722  fc = f;
723  gc = g;
724  if ( getCharacteristic() != 0 )
725  {
726    if ((!fc_and_gc_Univariate)
727    && isOn(SW_USE_fieldGCD)
728    && (getCharacteristic() >100))
729    {
730      return fieldGCD(f,g);
731    }
732    #ifdef HAVE_NTL
733    else if ((!fc_and_gc_Univariate) && (isOn( SW_USE_EZGCD_P )))
734    {
735      fc= EZGCD_P (fc, gc);
736    }
737    else if (isOn(SW_USE_FF_MOD_GCD) && !fc_and_gc_Univariate)
738    {
739      Variable a;
740      if (hasFirstAlgVar (fc, a) || hasFirstAlgVar (gc, a))
741        fc=GCD_Fp_extension (fc, gc, a);
742      else if (CFFactory::gettype() == GaloisFieldDomain)
743        fc=GCD_GF (fc, gc);
744      else
745        fc=GCD_small_p (fc, gc);
746    }
747    #endif
748    else if ( p1 == fc.level() )
749      fc = gcd_poly_p( fc, gc );
750    else
751    {
752      fc = replacevar( fc, Variable(p1), Variable(mp) );
753      gc = replacevar( gc, Variable(p1), Variable(mp) );
754      fc = replacevar( gcd_poly_p( fc, gc ), Variable(mp), Variable(p1) );
755    }
756  }
757  else if (!fc_and_gc_Univariate)
758  {
759    if ( isOn( SW_USE_EZGCD ) )
760    {
761      fc= ezgcd (fc, gc);
762      /*if ( pe == 1 )
763        fc = ezgcd( fc, gc );
764      else if ( pe > 0 )// no variable at position 1
765      {
766        fc = replacevar( fc, Variable(pe), Variable(1) );
767        gc = replacevar( gc, Variable(pe), Variable(1) );
768        fc = replacevar( ezgcd( fc, gc ), Variable(1), Variable(pe) );
769      }
770      else
771      {
772        pe = -pe;
773        fc = swapvar( fc, Variable(pe), Variable(1) );
774        gc = swapvar( gc, Variable(pe), Variable(1) );
775        fc = swapvar( ezgcd( fc, gc ), Variable(1), Variable(pe) );
776      }*/
777    }
778    else if (
779    isOn(SW_USE_CHINREM_GCD)
780    && (isPurePoly_m(fc)) && (isPurePoly_m(gc))
781    )
782    {
783    #if 0
784      if ( p1 == fc.level() )
785        fc = chinrem_gcd( fc, gc );
786      else
787      {
788        fc = replacevar( fc, Variable(p1), Variable(mp) );
789        gc = replacevar( gc, Variable(p1), Variable(mp) );
790        fc = replacevar( chinrem_gcd( fc, gc ), Variable(mp), Variable(p1) );
791      }
792    #else
793      fc = chinrem_gcd( fc, gc);
794    #endif
795    }
796    else
797    {
798       fc = gcd_poly_0( fc, gc );
799    }
800  }
801  else
802  {
803    fc = gcd_poly_0( fc, gc );
804  }
805  if ( d1.degree() > 0 )
806    fc *= d1;
807  return fc;
808}
809//}}}
810
811//{{{ static CanonicalForm cf_content ( const CanonicalForm & f, const CanonicalForm & g )
812//{{{ docu
813//
814// cf_content() - return gcd(g, content(f)).
815//
816// content(f) is calculated with respect to f's main variable.
817//
818// Used by gcd(), content(), content( CF, Variable ).
819//
820//}}}
821static CanonicalForm
822cf_content ( const CanonicalForm & f, const CanonicalForm & g )
823{
824    if ( f.inPolyDomain() || ( f.inExtension() && ! getReduce( f.mvar() ) ) )
825    {
826        CFIterator i = f;
827        CanonicalForm result = g;
828        while ( i.hasTerms() && ! result.isOne() )
829        {
830            result = gcd( i.coeff(), result );
831            i++;
832        }
833        return result;
834    }
835    else
836        return abs( f );
837}
838//}}}
839
840//{{{ CanonicalForm content ( const CanonicalForm & f )
841//{{{ docu
842//
843// content() - return content(f) with respect to main variable.
844//
845// Normalizes result.
846//
847//}}}
848CanonicalForm
849content ( const CanonicalForm & f )
850{
851    if ( f.inPolyDomain() || ( f.inExtension() && ! getReduce( f.mvar() ) ) )
852    {
853        CFIterator i = f;
854        CanonicalForm result = abs( i.coeff() );
855        i++;
856        while ( i.hasTerms() && ! result.isOne() )
857        {
858            result = gcd( i.coeff(), result );
859            i++;
860        }
861        return result;
862    }
863    else
864        return abs( f );
865}
866//}}}
867
868//{{{ CanonicalForm content ( const CanonicalForm & f, const Variable & x )
869//{{{ docu
870//
871// content() - return content(f) with respect to x.
872//
873// x should be a polynomial variable.
874//
875//}}}
876CanonicalForm
877content ( const CanonicalForm & f, const Variable & x )
878{
879    ASSERT( x.level() > 0, "cannot calculate content with respect to algebraic variable" );
880    Variable y = f.mvar();
881
882    if ( y == x )
883        return cf_content( f, 0 );
884    else  if ( y < x )
885        return f;
886    else
887        return swapvar( content( swapvar( f, y, x ), y ), y, x );
888}
889//}}}
890
891//{{{ CanonicalForm vcontent ( const CanonicalForm & f, const Variable & x )
892//{{{ docu
893//
894// vcontent() - return content of f with repect to variables >= x.
895//
896// The content is recursively calculated over all coefficients in
897// f having level less than x.  x should be a polynomial
898// variable.
899//
900//}}}
901CanonicalForm
902vcontent ( const CanonicalForm & f, const Variable & x )
903{
904    ASSERT( x.level() > 0, "cannot calculate vcontent with respect to algebraic variable" );
905
906    if ( f.mvar() <= x )
907        return content( f, x );
908    else {
909        CFIterator i;
910        CanonicalForm d = 0;
911        for ( i = f; i.hasTerms() && ! d.isOne(); i++ )
912            d = gcd( d, vcontent( i.coeff(), x ) );
913        return d;
914    }
915}
916//}}}
917
918//{{{ CanonicalForm pp ( const CanonicalForm & f )
919//{{{ docu
920//
921// pp() - return primitive part of f.
922//
923// Returns zero if f equals zero, otherwise f / content(f).
924//
925//}}}
926CanonicalForm
927pp ( const CanonicalForm & f )
928{
929    if ( f.isZero() )
930        return f;
931    else
932        return f / content( f );
933}
934//}}}
935
936CanonicalForm
937gcd ( const CanonicalForm & f, const CanonicalForm & g )
938{
939    bool b = f.isZero();
940    if ( b || g.isZero() )
941    {
942        if ( b )
943            return abs( g );
944        else
945            return abs( f );
946    }
947    if ( f.inPolyDomain() || g.inPolyDomain() )
948    {
949        if ( f.mvar() != g.mvar() )
950        {
951            if ( f.mvar() > g.mvar() )
952                return cf_content( f, g );
953            else
954                return cf_content( g, f );
955        }
956        if (isOn(SW_USE_QGCD))
957        {
958          Variable m;
959          if (
960          (getCharacteristic() == 0) &&
961          (hasFirstAlgVar(f,m) || hasFirstAlgVar(g,m))
962          )
963          {
964            bool on_rational = isOn(SW_RATIONAL);
965            CanonicalForm r=QGCD(f,g);
966            On(SW_RATIONAL);
967            CanonicalForm cdF = bCommonDen( r );
968            if (!on_rational) Off(SW_RATIONAL);
969            return cdF*r;
970          }
971        }
972
973        if ( f.inExtension() && getReduce( f.mvar() ) )
974            return CanonicalForm(1);
975        else
976        {
977            if ( fdivides( f, g ) )
978                return abs( f );
979            else  if ( fdivides( g, f ) )
980                return abs( g );
981            if ( !( getCharacteristic() == 0 && isOn( SW_RATIONAL ) ) )
982            {
983                CanonicalForm d;
984                d = gcd_poly( f, g );
985                return abs( d );
986            }
987            else
988            {
989                //printf ("here\n");
990                CanonicalForm cdF = bCommonDen( f );
991                CanonicalForm cdG = bCommonDen( g );
992                Off( SW_RATIONAL );
993                CanonicalForm l = lcm( cdF, cdG );
994                On( SW_RATIONAL );
995                CanonicalForm F = f * l, G = g * l;
996                Off( SW_RATIONAL );
997                l = gcd_poly( F, G );
998                On( SW_RATIONAL );
999                return abs( l );
1000            }
1001        }
1002    }
1003    if ( f.inBaseDomain() && g.inBaseDomain() )
1004        return bgcd( f, g );
1005    else
1006        return 1;
1007}
1008
1009//{{{ CanonicalForm lcm ( const CanonicalForm & f, const CanonicalForm & g )
1010//{{{ docu
1011//
1012// lcm() - return least common multiple of f and g.
1013//
1014// The lcm is calculated using the formula lcm(f, g) = f * g / gcd(f, g).
1015//
1016// Returns zero if one of f or g equals zero.
1017//
1018//}}}
1019CanonicalForm
1020lcm ( const CanonicalForm & f, const CanonicalForm & g )
1021{
1022    if ( f.isZero() || g.isZero() )
1023        return 0;
1024    else
1025        return ( f / gcd( f, g ) ) * g;
1026}
1027//}}}
1028
1029#ifdef HAVE_NTL
1030
1031static CanonicalForm
1032gcd_univar_ntl0( const CanonicalForm & F, const CanonicalForm & G )
1033{
1034    ZZX F1=convertFacCF2NTLZZX(F);
1035    ZZX G1=convertFacCF2NTLZZX(G);
1036    ZZX R=GCD(F1,G1);
1037    return convertNTLZZX2CF(R,F.mvar());
1038}
1039
1040static CanonicalForm
1041gcd_univar_ntlp( const CanonicalForm & F, const CanonicalForm & G )
1042{
1043  if (fac_NTL_char!=getCharacteristic())
1044  {
1045    fac_NTL_char=getCharacteristic();
1046    #ifdef NTL_ZZ
1047    ZZ r;
1048    r=getCharacteristic();
1049    ZZ_pContext ccc(r);
1050    #else
1051    zz_pContext ccc(getCharacteristic());
1052    #endif
1053    ccc.restore();
1054    #ifdef NTL_ZZ
1055    ZZ_p::init(r);
1056    #else
1057    zz_p::init(getCharacteristic());
1058    #endif
1059  }
1060  #ifdef NTL_ZZ
1061  ZZ_pX F1=convertFacCF2NTLZZpX(F);
1062  ZZ_pX G1=convertFacCF2NTLZZpX(G);
1063  ZZ_pX R=GCD(F1,G1);
1064  return  convertNTLZZpX2CF(R,F.mvar());
1065  #else
1066  zz_pX F1=convertFacCF2NTLzzpX(F);
1067  zz_pX G1=convertFacCF2NTLzzpX(G);
1068  zz_pX R=GCD(F1,G1);
1069  return  convertNTLzzpX2CF(R,F.mvar());
1070  #endif
1071}
1072
1073#endif
1074
1075#ifdef HAVE_FLINT
1076static CanonicalForm
1077gcd_univar_flintp (const CanonicalForm& F, const CanonicalForm& G)
1078{
1079  nmod_poly_t F1, G1;
1080  convertFacCF2nmod_poly_t (F1, F);
1081  convertFacCF2nmod_poly_t (G1, G);
1082  nmod_poly_gcd (F1, F1, G1);
1083  CanonicalForm result= convertnmod_poly_t2FacCF (F1, F.mvar());
1084  nmod_poly_clear (F1);
1085  nmod_poly_clear (G1);
1086  return result;
1087}
1088
1089static CanonicalForm
1090gcd_univar_flint0( const CanonicalForm & F, const CanonicalForm & G )
1091{
1092  fmpz_poly_t F1, G1;
1093  convertFacCF2Fmpz_poly_t(F1, F);
1094  convertFacCF2Fmpz_poly_t(G1, G);
1095  fmpz_poly_gcd (F1, F1, G1);
1096  CanonicalForm result= convertFmpz_poly_t2FacCF (F1, F.mvar());
1097  fmpz_poly_clear (F1);
1098  fmpz_poly_clear (G1);
1099  return result;
1100}
1101#endif
1102
1103
1104/*
1105*  compute positions p1 and pe of optimal variables:
1106*    pe is used in "ezgcd" and
1107*    p1 in "gcd_poly1"
1108*/
1109static
1110void optvalues ( const int * df, const int * dg, const int n, int & p1, int &pe )
1111{
1112    int i, o1, oe;
1113    if ( df[n] > dg[n] )
1114    {
1115        o1 = df[n]; oe = dg[n];
1116    }
1117    else
1118    {
1119        o1 = dg[n]; oe = df[n];
1120    }
1121    i = n-1;
1122    while ( i > 0 )
1123    {
1124        if ( df[i] != 0 )
1125        {
1126            if ( df[i] > dg[i] )
1127            {
1128                if ( o1 > df[i]) { o1 = df[i]; p1 = i; }
1129                if ( oe <= dg[i]) { oe = dg[i]; pe = i; }
1130            }
1131            else
1132            {
1133                if ( o1 > dg[i]) { o1 = dg[i]; p1 = i; }
1134                if ( oe <= df[i]) { oe = df[i]; pe = i; }
1135            }
1136        }
1137        i--;
1138    }
1139}
1140
1141/*
1142*  make some changes of variables, see optvalues
1143*/
1144static void
1145cf_prepgcd( const CanonicalForm & f, const CanonicalForm & g, int & cc, int & p1, int &pe )
1146{
1147    int i, k, n;
1148    n = f.level();
1149    cc = 0;
1150    p1 = pe = n;
1151    if ( n == 1 )
1152        return;
1153    int * degsf = new int[n+1];
1154    int * degsg = new int[n+1];
1155    for ( i = n; i > 0; i-- )
1156    {
1157        degsf[i] = degsg[i] = 0;
1158    }
1159    degsf = degrees( f, degsf );
1160    degsg = degrees( g, degsg );
1161
1162    k = 0;
1163    for ( i = n-1; i > 0; i-- )
1164    {
1165        if ( degsf[i] == 0 )
1166        {
1167            if ( degsg[i] != 0 )
1168            {
1169                cc = -i;
1170                break;
1171            }
1172        }
1173        else
1174        {
1175            if ( degsg[i] == 0 )
1176            {
1177                cc = i;
1178                break;
1179            }
1180            else k++;
1181        }
1182    }
1183
1184    if ( ( cc == 0 ) && ( k != 0 ) )
1185        optvalues( degsf, degsg, n, p1, pe );
1186    if ( ( pe != 1 ) && ( degsf[1] != 0 ) )
1187        pe = -pe;
1188
1189    delete [] degsf;
1190    delete [] degsg;
1191}
1192
1193
1194static CanonicalForm
1195balance_p ( const CanonicalForm & f, const CanonicalForm & q )
1196{
1197    Variable x = f.mvar();
1198    CanonicalForm result = 0, qh = q / 2;
1199    CanonicalForm c;
1200    CFIterator i;
1201    for ( i = f; i.hasTerms(); i++ )
1202    {
1203        c = i.coeff();
1204        if ( c.inCoeffDomain())
1205        {
1206          if ( c > qh )
1207            result += power( x, i.exp() ) * (c - q);
1208          else
1209            result += power( x, i.exp() ) * c;
1210        }
1211        else
1212          result += power( x, i.exp() ) * balance_p(c,q);
1213    }
1214    return result;
1215}
1216
1217CanonicalForm chinrem_gcd ( const CanonicalForm & FF, const CanonicalForm & GG )
1218{
1219  CanonicalForm f, g, cl, q(0), Dp, newD, D, newq;
1220  int p, i, dp_deg, d_deg=-1;
1221
1222  CanonicalForm cd ( bCommonDen( FF ));
1223  f=cd*FF;
1224  Variable x= Variable (1);
1225  CanonicalForm cf, cg;
1226  cf= icontent (f);
1227  f /= cf;
1228  //cd = bCommonDen( f ); f *=cd;
1229  //f /=vcontent(f,Variable(1));
1230
1231  cd = bCommonDen( GG );
1232  g=cd*GG;
1233  cg= icontent (g);
1234  g /= cg;
1235  //cd = bCommonDen( g ); g *=cd;
1236  //g /=vcontent(g,Variable(1));
1237
1238  CanonicalForm Dn, test= 0;
1239  cl =  gcd (f.lc(),g.lc());
1240  CanonicalForm b= 1;
1241  int minCommonDeg= 0;
1242  CanonicalForm gcdcfcg= gcd (cf, cg);
1243  for (i= tmax (f.level(), g.level()); i > 0; i--)
1244  {
1245    if (degree (f, i) <= 0 || degree (g, i) <= 0)
1246      continue;
1247    else
1248    {
1249      minCommonDeg= tmin (degree (g, i), degree (f, i));
1250      break;
1251    }
1252  }
1253  if (i == 0)
1254    return gcdcfcg;
1255  for (; i > 0; i--)
1256  {
1257    if (degree (f, i) <= 0 || degree (g, i) <= 0)
1258      continue;
1259    else
1260      minCommonDeg= tmin (minCommonDeg, tmin (degree (g, i), degree (f, i)));
1261  }
1262  b= 2*tmin (maxNorm (f), maxNorm (g))*abs (cl)*power (CanonicalForm (2), minCommonDeg);
1263  bool equal= false;
1264  i = cf_getNumBigPrimes() - 1;
1265
1266  //Off (SW_RATIONAL);
1267  while ( true )
1268  {
1269    p = cf_getBigPrime( i );
1270    i--;
1271    while ( i >= 0 && mod( cl, p ) == 0 )
1272    {
1273      p = cf_getBigPrime( i );
1274      i--;
1275    }
1276    //printf("try p=%d\n",p);
1277    setCharacteristic( p );
1278    Dp = gcd_poly( mapinto( f ), mapinto( g ) );
1279    Dp /=Dp.lc();
1280    setCharacteristic( 0 );
1281    dp_deg=totaldegree(Dp);
1282    if ( dp_deg == 0 )
1283    {
1284      //printf(" -> 1\n");
1285      return CanonicalForm(gcdcfcg);
1286    }
1287    if ( q.isZero() )
1288    {
1289      D = mapinto( Dp );
1290      d_deg=dp_deg;
1291      q = p;
1292    }
1293    else
1294    {
1295      if ( dp_deg == d_deg )
1296      {
1297        chineseRemainder( D, q, mapinto( Dp ), p, newD, newq );
1298        q = newq;
1299        D = newD;
1300      }
1301      else if ( dp_deg < d_deg )
1302      {
1303        //printf(" were all bad, try more\n");
1304        // all previous p's are bad primes
1305        q = p;
1306        D = mapinto( Dp );
1307        d_deg=dp_deg;
1308        test= 0;
1309        equal= false;
1310      }
1311      else
1312      {
1313        //printf(" was bad, try more\n");
1314      }
1315      //else dp_deg > d_deg: bad prime
1316    }
1317    if ( i >= 0 )
1318    {
1319      Dn= Farey(D,q);
1320      int is_rat= isOn (SW_RATIONAL);
1321      On (SW_RATIONAL);
1322      cd = bCommonDen( Dn ); // we need On(SW_RATIONAL)
1323      if (!is_rat)
1324        Off (SW_RATIONAL);
1325      Dn *=cd;
1326      if (test != Dn)
1327        test= Dn;
1328      else
1329        equal= true;
1330      //Dn /=vcontent(Dn,Variable(1));
1331      if ((equal || q > b) && fdivides( Dn, f ) && fdivides( Dn, g ) )
1332      {
1333        //printf(" -> success\n");
1334        return Dn*gcdcfcg;
1335      }
1336      equal= false;
1337      //else: try more primes
1338    }
1339    else
1340    { // try other method
1341      //printf("try other gcd\n");
1342      Off(SW_USE_CHINREM_GCD);
1343      D=gcd_poly( f, g );
1344      On(SW_USE_CHINREM_GCD);
1345      return D*gcdcfcg;
1346    }
1347  }
1348}
1349
Note: See TracBrowser for help on using the repository browser.