source: git/factory/cf_gcd.cc @ 7061c1

spielwiese
Last change on this file since 7061c1 was 8c00057, checked in by Martin Lee <martinlee84@…>, 11 years ago
chg: normalize gcd in extgcd computation via NTL
  • 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 ( const CanonicalForm & f, const CanonicalForm & q )
466//{{{ docu
467//
468// balance() - 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 ( 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        c = mod( i.coeff(), q );
486        if ( c > qh )
487            result += power( x, i.exp() ) * (c - q);
488        else
489            result += power( x, i.exp() ) * c;
490    }
491    return result;
492}
493//}}}
494
495static CanonicalForm gcd_poly_univar0( const CanonicalForm & F, const CanonicalForm & G, bool primitive )
496{
497  CanonicalForm f, g, c, cg, cl, BB, B, M, q, Dp, newD, D, newq;
498  int p, i;
499
500  if ( primitive )
501  {
502    f = F;
503    g = G;
504    c = 1;
505  }
506  else
507  {
508    CanonicalForm cF = content( F ), cG = content( G );
509    f = F / cF;
510    g = G / cG;
511    c = bgcd( cF, cG );
512  }
513  cg = gcd( f.lc(), g.lc() );
514  cl = ( f.lc() / cg ) * g.lc();
515//     B = 2 * cg * tmin(
516//         maxnorm(f)*power(CanonicalForm(2),f.degree())*isqrt(f.degree()+1),
517//         maxnorm(g)*power(CanonicalForm(2),g.degree())*isqrt(g.degree()+1)
518//         )+1;
519  M = tmin( maxNorm(f), maxNorm(g) );
520  BB = power(CanonicalForm(2),tmin(f.degree(),g.degree()))*M;
521  q = 0;
522  i = cf_getNumSmallPrimes() - 1;
523  while ( true )
524  {
525    B = BB;
526    while ( i >= 0 && q < B )
527    {
528      p = cf_getSmallPrime( i );
529      i--;
530      while ( i >= 0 && mod( cl, p ) == 0 )
531      {
532        p = cf_getSmallPrime( i );
533        i--;
534      }
535      setCharacteristic( p );
536      Dp = gcd( mapinto( f ), mapinto( g ) );
537      Dp = ( Dp / Dp.lc() ) * mapinto( cg );
538      setCharacteristic( 0 );
539      if ( Dp.degree() == 0 )
540        return c;
541      if ( q.isZero() )
542      {
543        D = mapinto( Dp );
544        q = p;
545        B = power(CanonicalForm(2),D.degree())*M+1;
546      }
547      else
548      {
549        if ( Dp.degree() == D.degree() )
550        {
551          chineseRemainder( D, q, mapinto( Dp ), p, newD, newq );
552          q = newq;
553          D = newD;
554        }
555        else if ( Dp.degree() < D.degree() )
556        {
557          // all previous p's are bad primes
558          q = p;
559          D = mapinto( Dp );
560          B = power(CanonicalForm(2),D.degree())*M+1;
561        }
562        // else p is a bad prime
563      }
564    }
565    if ( i >= 0 )
566    {
567      // now balance D mod q
568      D = pp( balance( D, q ) );
569      if ( fdivides( D, f ) && fdivides( D, g ) )
570        return D * c;
571      else
572        q = 0;
573    }
574    else
575      return gcd_poly( F, G );
576    DEBOUTLN( cerr, "another try ..." );
577  }
578}
579
580static CanonicalForm
581gcd_poly_p( const CanonicalForm & f, const CanonicalForm & g )
582{
583    CanonicalForm pi, pi1;
584    CanonicalForm C, Ci, Ci1, Hi, bi, pi2;
585    bool bpure;
586    int delta = degree( f ) - degree( g );
587
588    if ( delta >= 0 )
589    {
590        pi = f; pi1 = g;
591    }
592    else
593    {
594        pi = g; pi1 = f; delta = -delta;
595    }
596    Ci = content( pi ); Ci1 = content( pi1 );
597    pi1 = pi1 / Ci1; pi = pi / Ci;
598    C = gcd( Ci, Ci1 );
599    if ( !( pi.isUnivariate() && pi1.isUnivariate() ) )
600    {
601        if ( gcd_test_one( pi1, pi, true ) )
602        {
603          C=abs(C);
604          //out_cf("GCD:",C,"\n");
605          return C;
606        }
607        bpure = false;
608    }
609    else
610    {
611        bpure = isPurePoly(pi) && isPurePoly(pi1);
612#ifdef HAVE_FLINT
613        if (bpure && (CFFactory::gettype() != GaloisFieldDomain))
614          return gcd_univar_flintp(pi,pi1)*C;
615#else
616#ifdef HAVE_NTL
617        if ( isOn(SW_USE_NTL_GCD_P) && bpure && (CFFactory::gettype() != GaloisFieldDomain))
618            return gcd_univar_ntlp(pi, pi1 ) * C;
619#endif
620#endif
621    }
622    Variable v = f.mvar();
623    Hi = power( LC( pi1, v ), delta );
624    int maxNumVars= tmax (getNumVars (pi), getNumVars (pi1));
625
626    if (!(pi.isUnivariate() && pi1.isUnivariate()))
627    {
628      if (size (Hi)*size (pi)/(maxNumVars*3) > 500) //maybe this needs more tuning
629      {
630        On (SW_USE_FF_MOD_GCD);
631        C *= gcd (pi, pi1);
632        Off (SW_USE_FF_MOD_GCD);
633        return C;
634      }
635    }
636
637    if ( (delta+1) % 2 )
638        bi = 1;
639    else
640        bi = -1;
641    CanonicalForm oldPi= pi, oldPi1= pi1, powHi;
642    while ( degree( pi1, v ) > 0 )
643    {
644        if (!(pi.isUnivariate() && pi1.isUnivariate()))
645        {
646          if (size (pi)/maxNumVars > 500 || size (pi1)/maxNumVars > 500)
647          {
648            On (SW_USE_FF_MOD_GCD);
649            C *= gcd (oldPi, oldPi1);
650            Off (SW_USE_FF_MOD_GCD);
651            return C;
652          }
653        }
654        pi2 = psr( pi, pi1, v );
655        pi2 = pi2 / bi;
656        pi = pi1; pi1 = pi2;
657        maxNumVars= tmax (getNumVars (pi), getNumVars (pi1));
658        if (!(pi1.isUnivariate()) && (size (pi1)/maxNumVars > 500))
659        {
660            On (SW_USE_FF_MOD_GCD);
661            C *= gcd (oldPi, oldPi1);
662            Off (SW_USE_FF_MOD_GCD);
663            return C;
664        }
665        if ( degree( pi1, v ) > 0 )
666        {
667            delta = degree( pi, v ) - degree( pi1, v );
668            powHi= power (Hi, delta-1);
669            if ( (delta+1) % 2 )
670                bi = LC( pi, v ) * powHi*Hi;
671            else
672                bi = -LC( pi, v ) * powHi*Hi;
673            Hi = power( LC( pi1, v ), delta ) / powHi;
674            if (!(pi.isUnivariate() && pi1.isUnivariate()))
675            {
676              if (size (Hi)*size (pi)/(maxNumVars*3) > 500) //maybe this needs more tuning
677              {
678                On (SW_USE_FF_MOD_GCD);
679                C *= gcd (oldPi, oldPi1);
680                Off (SW_USE_FF_MOD_GCD);
681                return C;
682              }
683            }
684        }
685    }
686    if ( degree( pi1, v ) == 0 )
687    {
688      C=abs(C);
689      //out_cf("GCD:",C,"\n");
690      return C;
691    }
692    pi /= content( pi );
693    if ( bpure )
694        pi /= pi.lc();
695    C=abs(C*pi);
696    //out_cf("GCD:",C,"\n");
697    return C;
698}
699
700static CanonicalForm
701gcd_poly_0( const CanonicalForm & f, const CanonicalForm & g )
702{
703    CanonicalForm pi, pi1;
704    CanonicalForm C, Ci, Ci1, Hi, bi, pi2;
705    int delta = degree( f ) - degree( g );
706
707    if ( delta >= 0 )
708    {
709        pi = f; pi1 = g;
710    }
711    else
712    {
713        pi = g; pi1 = f; delta = -delta;
714    }
715    Ci = content( pi ); Ci1 = content( pi1 );
716    pi1 = pi1 / Ci1; pi = pi / Ci;
717    C = gcd( Ci, Ci1 );
718    if ( pi.isUnivariate() && pi1.isUnivariate() )
719    {
720#ifdef HAVE_FLINT
721        if (isPurePoly(pi) && isPurePoly(pi1) )
722            return gcd_univar_flint0(pi, pi1 ) * C;
723#else
724#ifdef HAVE_NTL
725        if ( isOn(SW_USE_NTL_GCD_0) && isPurePoly(pi) && isPurePoly(pi1) )
726            return gcd_univar_ntl0(pi, pi1 ) * C;
727#endif
728#endif
729        return gcd_poly_univar0( pi, pi1, true ) * C;
730    }
731    else if ( gcd_test_one( pi1, pi, true ) )
732      return C;
733    Variable v = f.mvar();
734    Hi = power( LC( pi1, v ), delta );
735    if ( (delta+1) % 2 )
736        bi = 1;
737    else
738        bi = -1;
739    while ( degree( pi1, v ) > 0 )
740    {
741        pi2 = psr( pi, pi1, v );
742        pi2 = pi2 / bi;
743        pi = pi1; pi1 = pi2;
744        if ( degree( pi1, v ) > 0 )
745        {
746            delta = degree( pi, v ) - degree( pi1, v );
747            if ( (delta+1) % 2 )
748                bi = LC( pi, v ) * power( Hi, delta );
749            else
750                bi = -LC( pi, v ) * power( Hi, delta );
751            Hi = power( LC( pi1, v ), delta ) / power( Hi, delta-1 );
752        }
753    }
754    if ( degree( pi1, v ) == 0 )
755        return C;
756    else
757        return C * pp( pi );
758}
759
760//{{{ CanonicalForm gcd_poly ( const CanonicalForm & f, const CanonicalForm & g )
761//{{{ docu
762//
763// gcd_poly() - calculate polynomial gcd.
764//
765// This is the dispatcher for polynomial gcd calculation.  We call either
766// ezgcd(), sparsemod() or gcd_poly1() in dependecy on the current
767// characteristic and settings of SW_USE_EZGCD.
768//
769// Used by gcd() and gcd_poly_univar0().
770//
771//}}}
772#if 0
773int si_factor_reminder=1;
774#endif
775CanonicalForm gcd_poly ( const CanonicalForm & f, const CanonicalForm & g )
776{
777  CanonicalForm fc, gc, d1;
778  bool fc_isUnivariate=f.isUnivariate();
779  bool gc_isUnivariate=g.isUnivariate();
780  bool fc_and_gc_Univariate=fc_isUnivariate && gc_isUnivariate;
781  fc = f;
782  gc = g;
783  if ( getCharacteristic() != 0 )
784  {
785    #ifdef HAVE_NTL
786    if ((!fc_and_gc_Univariate) && (isOn( SW_USE_EZGCD_P )))
787    {
788      fc= EZGCD_P (fc, gc);
789    }
790    else if (isOn(SW_USE_FF_MOD_GCD) && !fc_and_gc_Univariate)
791    {
792      Variable a;
793      if (hasFirstAlgVar (fc, a) || hasFirstAlgVar (gc, a))
794        fc=GCD_Fp_extension (fc, gc, a);
795      else if (CFFactory::gettype() == GaloisFieldDomain)
796        fc=GCD_GF (fc, gc);
797      else
798        fc=GCD_small_p (fc, gc);
799    }
800    else
801    #endif
802    fc = gcd_poly_p( fc, gc );
803  }
804  else if (!fc_and_gc_Univariate)
805  {
806    if ( isOn( SW_USE_EZGCD ) )
807      fc= ezgcd (fc, gc);
808    else if (isOn(SW_USE_CHINREM_GCD))
809      fc = chinrem_gcd( fc, gc);
810    else
811    {
812       fc = gcd_poly_0( fc, gc );
813    }
814  }
815  else
816  {
817    fc = gcd_poly_0( fc, gc );
818  }
819  if ( d1.degree() > 0 )
820    fc *= d1;
821  return fc;
822}
823//}}}
824
825//{{{ static CanonicalForm cf_content ( const CanonicalForm & f, const CanonicalForm & g )
826//{{{ docu
827//
828// cf_content() - return gcd(g, content(f)).
829//
830// content(f) is calculated with respect to f's main variable.
831//
832// Used by gcd(), content(), content( CF, Variable ).
833//
834//}}}
835static CanonicalForm
836cf_content ( const CanonicalForm & f, const CanonicalForm & g )
837{
838    if ( f.inPolyDomain() || ( f.inExtension() && ! getReduce( f.mvar() ) ) )
839    {
840        CFIterator i = f;
841        CanonicalForm result = g;
842        while ( i.hasTerms() && ! result.isOne() )
843        {
844            result = gcd( i.coeff(), result );
845            i++;
846        }
847        return result;
848    }
849    else
850        return abs( f );
851}
852//}}}
853
854//{{{ CanonicalForm content ( const CanonicalForm & f )
855//{{{ docu
856//
857// content() - return content(f) with respect to main variable.
858//
859// Normalizes result.
860//
861//}}}
862CanonicalForm
863content ( const CanonicalForm & f )
864{
865    if ( f.inPolyDomain() || ( f.inExtension() && ! getReduce( f.mvar() ) ) )
866    {
867        CFIterator i = f;
868        CanonicalForm result = abs( i.coeff() );
869        i++;
870        while ( i.hasTerms() && ! result.isOne() )
871        {
872            result = gcd( i.coeff(), result );
873            i++;
874        }
875        return result;
876    }
877    else
878        return abs( f );
879}
880//}}}
881
882//{{{ CanonicalForm content ( const CanonicalForm & f, const Variable & x )
883//{{{ docu
884//
885// content() - return content(f) with respect to x.
886//
887// x should be a polynomial variable.
888//
889//}}}
890CanonicalForm
891content ( const CanonicalForm & f, const Variable & x )
892{
893    ASSERT( x.level() > 0, "cannot calculate content with respect to algebraic variable" );
894    Variable y = f.mvar();
895
896    if ( y == x )
897        return cf_content( f, 0 );
898    else  if ( y < x )
899        return f;
900    else
901        return swapvar( content( swapvar( f, y, x ), y ), y, x );
902}
903//}}}
904
905//{{{ CanonicalForm vcontent ( const CanonicalForm & f, const Variable & x )
906//{{{ docu
907//
908// vcontent() - return content of f with repect to variables >= x.
909//
910// The content is recursively calculated over all coefficients in
911// f having level less than x.  x should be a polynomial
912// variable.
913//
914//}}}
915CanonicalForm
916vcontent ( const CanonicalForm & f, const Variable & x )
917{
918    ASSERT( x.level() > 0, "cannot calculate vcontent with respect to algebraic variable" );
919
920    if ( f.mvar() <= x )
921        return content( f, x );
922    else {
923        CFIterator i;
924        CanonicalForm d = 0;
925        for ( i = f; i.hasTerms() && ! d.isOne(); i++ )
926            d = gcd( d, vcontent( i.coeff(), x ) );
927        return d;
928    }
929}
930//}}}
931
932//{{{ CanonicalForm pp ( const CanonicalForm & f )
933//{{{ docu
934//
935// pp() - return primitive part of f.
936//
937// Returns zero if f equals zero, otherwise f / content(f).
938//
939//}}}
940CanonicalForm
941pp ( const CanonicalForm & f )
942{
943    if ( f.isZero() )
944        return f;
945    else
946        return f / content( f );
947}
948//}}}
949
950CanonicalForm
951gcd ( const CanonicalForm & f, const CanonicalForm & g )
952{
953    bool b = f.isZero();
954    if ( b || g.isZero() )
955    {
956        if ( b )
957            return abs( g );
958        else
959            return abs( f );
960    }
961    if ( f.inPolyDomain() || g.inPolyDomain() )
962    {
963        if ( f.mvar() != g.mvar() )
964        {
965            if ( f.mvar() > g.mvar() )
966                return cf_content( f, g );
967            else
968                return cf_content( g, f );
969        }
970        if (isOn(SW_USE_QGCD))
971        {
972          Variable m;
973          if (
974          (getCharacteristic() == 0) &&
975          (hasFirstAlgVar(f,m) || hasFirstAlgVar(g,m))
976          )
977          {
978            bool on_rational = isOn(SW_RATIONAL);
979            CanonicalForm r=QGCD(f,g);
980            On(SW_RATIONAL);
981            CanonicalForm cdF = bCommonDen( r );
982            if (!on_rational) Off(SW_RATIONAL);
983            return cdF*r;
984          }
985        }
986
987        if ( f.inExtension() && getReduce( f.mvar() ) )
988            return CanonicalForm(1);
989        else
990        {
991            if ( fdivides( f, g ) )
992                return abs( f );
993            else  if ( fdivides( g, f ) )
994                return abs( g );
995            if ( !( getCharacteristic() == 0 && isOn( SW_RATIONAL ) ) )
996            {
997                CanonicalForm d;
998                d = gcd_poly( f, g );
999                return abs( d );
1000            }
1001            else
1002            {
1003                //printf ("here\n");
1004                CanonicalForm cdF = bCommonDen( f );
1005                CanonicalForm cdG = bCommonDen( g );
1006                Off( SW_RATIONAL );
1007                CanonicalForm l = lcm( cdF, cdG );
1008                On( SW_RATIONAL );
1009                CanonicalForm F = f * l, G = g * l;
1010                Off( SW_RATIONAL );
1011                l = gcd_poly( F, G );
1012                On( SW_RATIONAL );
1013                return abs( l );
1014            }
1015        }
1016    }
1017    if ( f.inBaseDomain() && g.inBaseDomain() )
1018        return bgcd( f, g );
1019    else
1020        return 1;
1021}
1022
1023//{{{ CanonicalForm lcm ( const CanonicalForm & f, const CanonicalForm & g )
1024//{{{ docu
1025//
1026// lcm() - return least common multiple of f and g.
1027//
1028// The lcm is calculated using the formula lcm(f, g) = f * g / gcd(f, g).
1029//
1030// Returns zero if one of f or g equals zero.
1031//
1032//}}}
1033CanonicalForm
1034lcm ( const CanonicalForm & f, const CanonicalForm & g )
1035{
1036    if ( f.isZero() || g.isZero() )
1037        return 0;
1038    else
1039        return ( f / gcd( f, g ) ) * g;
1040}
1041//}}}
1042
1043#ifdef HAVE_NTL
1044#ifndef HAVE_FLINT
1045static CanonicalForm
1046gcd_univar_ntl0( const CanonicalForm & F, const CanonicalForm & G )
1047{
1048    ZZX F1=convertFacCF2NTLZZX(F);
1049    ZZX G1=convertFacCF2NTLZZX(G);
1050    ZZX R=GCD(F1,G1);
1051    return convertNTLZZX2CF(R,F.mvar());
1052}
1053
1054static CanonicalForm
1055gcd_univar_ntlp( const CanonicalForm & F, const CanonicalForm & G )
1056{
1057  if (fac_NTL_char!=getCharacteristic())
1058  {
1059    fac_NTL_char=getCharacteristic();
1060    #ifdef NTL_ZZ
1061    ZZ r;
1062    r=getCharacteristic();
1063    ZZ_pContext ccc(r);
1064    #else
1065    zz_pContext ccc(getCharacteristic());
1066    #endif
1067    ccc.restore();
1068    #ifdef NTL_ZZ
1069    ZZ_p::init(r);
1070    #else
1071    zz_p::init(getCharacteristic());
1072    #endif
1073  }
1074  #ifdef NTL_ZZ
1075  ZZ_pX F1=convertFacCF2NTLZZpX(F);
1076  ZZ_pX G1=convertFacCF2NTLZZpX(G);
1077  ZZ_pX R=GCD(F1,G1);
1078  return  convertNTLZZpX2CF(R,F.mvar());
1079  #else
1080  zz_pX F1=convertFacCF2NTLzzpX(F);
1081  zz_pX G1=convertFacCF2NTLzzpX(G);
1082  zz_pX R=GCD(F1,G1);
1083  return  convertNTLzzpX2CF(R,F.mvar());
1084  #endif
1085}
1086#endif
1087#endif
1088
1089#ifdef HAVE_FLINT
1090static CanonicalForm
1091gcd_univar_flintp (const CanonicalForm& F, const CanonicalForm& G)
1092{
1093  nmod_poly_t F1, G1;
1094  convertFacCF2nmod_poly_t (F1, F);
1095  convertFacCF2nmod_poly_t (G1, G);
1096  nmod_poly_gcd (F1, F1, G1);
1097  CanonicalForm result= convertnmod_poly_t2FacCF (F1, F.mvar());
1098  nmod_poly_clear (F1);
1099  nmod_poly_clear (G1);
1100  return result;
1101}
1102
1103static CanonicalForm
1104gcd_univar_flint0( const CanonicalForm & F, const CanonicalForm & G )
1105{
1106  fmpz_poly_t F1, G1;
1107  convertFacCF2Fmpz_poly_t(F1, F);
1108  convertFacCF2Fmpz_poly_t(G1, G);
1109  fmpz_poly_gcd (F1, F1, G1);
1110  CanonicalForm result= convertFmpz_poly_t2FacCF (F1, F.mvar());
1111  fmpz_poly_clear (F1);
1112  fmpz_poly_clear (G1);
1113  return result;
1114}
1115#endif
1116
1117
1118/*
1119*  compute positions p1 and pe of optimal variables:
1120*    pe is used in "ezgcd" and
1121*    p1 in "gcd_poly1"
1122*/
1123/*static
1124void optvalues ( const int * df, const int * dg, const int n, int & p1, int &pe )
1125{
1126    int i, o1, oe;
1127    if ( df[n] > dg[n] )
1128    {
1129        o1 = df[n]; oe = dg[n];
1130    }
1131    else
1132    {
1133        o1 = dg[n]; oe = df[n];
1134    }
1135    i = n-1;
1136    while ( i > 0 )
1137    {
1138        if ( df[i] != 0 )
1139        {
1140            if ( df[i] > dg[i] )
1141            {
1142                if ( o1 > df[i]) { o1 = df[i]; p1 = i; }
1143                if ( oe <= dg[i]) { oe = dg[i]; pe = i; }
1144            }
1145            else
1146            {
1147                if ( o1 > dg[i]) { o1 = dg[i]; p1 = i; }
1148                if ( oe <= df[i]) { oe = df[i]; pe = i; }
1149            }
1150        }
1151        i--;
1152    }
1153}*/
1154
1155/*
1156*  make some changes of variables, see optvalues
1157*/
1158/*static void
1159cf_prepgcd( const CanonicalForm & f, const CanonicalForm & g, int & cc, int & p1, int &pe )
1160{
1161    int i, k, n;
1162    n = f.level();
1163    cc = 0;
1164    p1 = pe = n;
1165    if ( n == 1 )
1166        return;
1167    int * degsf = new int[n+1];
1168    int * degsg = new int[n+1];
1169    for ( i = n; i > 0; i-- )
1170    {
1171        degsf[i] = degsg[i] = 0;
1172    }
1173    degsf = degrees( f, degsf );
1174    degsg = degrees( g, degsg );
1175
1176    k = 0;
1177    for ( i = n-1; i > 0; i-- )
1178    {
1179        if ( degsf[i] == 0 )
1180        {
1181            if ( degsg[i] != 0 )
1182            {
1183                cc = -i;
1184                break;
1185            }
1186        }
1187        else
1188        {
1189            if ( degsg[i] == 0 )
1190            {
1191                cc = i;
1192                break;
1193            }
1194            else k++;
1195        }
1196    }
1197
1198    if ( ( cc == 0 ) && ( k != 0 ) )
1199        optvalues( degsf, degsg, n, p1, pe );
1200    if ( ( pe != 1 ) && ( degsf[1] != 0 ) )
1201        pe = -pe;
1202
1203    delete [] degsf;
1204    delete [] degsg;
1205}*/
1206
1207
1208static CanonicalForm
1209balance_p ( const CanonicalForm & f, const CanonicalForm & q )
1210{
1211    Variable x = f.mvar();
1212    CanonicalForm result = 0, qh = q / 2;
1213    CanonicalForm c;
1214    CFIterator i;
1215    for ( i = f; i.hasTerms(); i++ )
1216    {
1217        c = i.coeff();
1218        if ( c.inCoeffDomain())
1219        {
1220          if ( c > qh )
1221            result += power( x, i.exp() ) * (c - q);
1222          else
1223            result += power( x, i.exp() ) * c;
1224        }
1225        else
1226          result += power( x, i.exp() ) * balance_p(c,q);
1227    }
1228    return result;
1229}
1230
1231TIMING_DEFINE_PRINT(chinrem_termination)
1232TIMING_DEFINE_PRINT(chinrem_recursion)
1233TIMING_DEFINE_PRINT(chinrem_reconstruction)
1234
1235CanonicalForm chinrem_gcd ( const CanonicalForm & FF, const CanonicalForm & GG )
1236{
1237  CanonicalForm f, g, cl, q(0), Dp, newD, D, newq;
1238  int p, i, dp_deg, d_deg=-1;
1239
1240  CanonicalForm cd ( bCommonDen( FF ));
1241  f=cd*FF;
1242  Variable x= Variable (1);
1243  CanonicalForm cf, cg;
1244  cf= icontent (f);
1245  f /= cf;
1246  //cd = bCommonDen( f ); f *=cd;
1247  //f /=vcontent(f,Variable(1));
1248
1249  cd = bCommonDen( GG );
1250  g=cd*GG;
1251  cg= icontent (g);
1252  g /= cg;
1253  //cd = bCommonDen( g ); g *=cd;
1254  //g /=vcontent(g,Variable(1));
1255
1256  CanonicalForm Dn, test= 0;
1257  cl =  gcd (f.lc(),g.lc());
1258  CanonicalForm gcdcfcg= gcd (cf, cg);
1259  CanonicalForm fp, gp;
1260  CanonicalForm b= 1;
1261  int minCommonDeg= 0;
1262  for (i= tmax (f.level(), g.level()); i > 0; i--)
1263  {
1264    if (degree (f, i) <= 0 || degree (g, i) <= 0)
1265      continue;
1266    else
1267    {
1268      minCommonDeg= tmin (degree (g, i), degree (f, i));
1269      break;
1270    }
1271  }
1272  if (i == 0)
1273    return gcdcfcg;
1274  for (; i > 0; i--)
1275  {
1276    if (degree (f, i) <= 0 || degree (g, i) <= 0)
1277      continue;
1278    else
1279      minCommonDeg= tmin (minCommonDeg, tmin (degree (g, i), degree (f, i)));
1280  }
1281  b= 2*tmin (maxNorm (f), maxNorm (g))*abs (cl)*
1282     power (CanonicalForm (2), minCommonDeg);
1283  bool equal= false;
1284  i = cf_getNumBigPrimes() - 1;
1285
1286  CanonicalForm cof, cog, cofp, cogp, newCof, newCog, cofn, cogn;
1287  int maxNumVars= tmax (getNumVars (f), getNumVars (g));
1288  //Off (SW_RATIONAL);
1289  while ( true )
1290  {
1291    p = cf_getBigPrime( i );
1292    i--;
1293    while ( i >= 0 && mod( cl*(lc(f)/cl)*(lc(g)/cl), p ) == 0 )
1294    {
1295      p = cf_getBigPrime( i );
1296      i--;
1297    }
1298    //printf("try p=%d\n",p);
1299    setCharacteristic( p );
1300    fp= mapinto (f);
1301    gp= mapinto (g);
1302    TIMING_START (chinrem_recursion)
1303#ifdef HAVE_NTL
1304    if (size (fp)/maxNumVars > 500 && size (gp)/maxNumVars > 500)
1305      Dp = GCD_small_p (fp, gp, cofp, cogp);
1306    else
1307    {
1308      Dp= gcd_poly (fp, gp);
1309      cofp= fp/Dp;
1310      cogp= gp/Dp;
1311    }
1312#else
1313    Dp= gcd_poly (fp, gp);
1314    cofp= fp/Dp;
1315    cogp= gp/Dp;
1316#endif
1317    TIMING_END_AND_PRINT (chinrem_recursion,
1318                          "time for gcd mod p in modular gcd: ");
1319    Dp /=Dp.lc();
1320    cofp /= lc (cofp);
1321    cogp /= lc (cogp);
1322    setCharacteristic( 0 );
1323    dp_deg=totaldegree(Dp);
1324    if ( dp_deg == 0 )
1325    {
1326      //printf(" -> 1\n");
1327      return CanonicalForm(gcdcfcg);
1328    }
1329    if ( q.isZero() )
1330    {
1331      D = mapinto( Dp );
1332      cof= mapinto (cofp);
1333      cog= mapinto (cogp);
1334      d_deg=dp_deg;
1335      q = p;
1336    }
1337    else
1338    {
1339      if ( dp_deg == d_deg )
1340      {
1341        chineseRemainder( D, q, mapinto( Dp ), p, newD, newq );
1342        chineseRemainder( cof, q, mapinto (cofp), p, newCof, newq);
1343        chineseRemainder( cog, q, mapinto (cogp), p, newCog, newq);
1344        cof= newCof;
1345        cog= newCog;
1346        q = newq;
1347        D = newD;
1348      }
1349      else if ( dp_deg < d_deg )
1350      {
1351        //printf(" were all bad, try more\n");
1352        // all previous p's are bad primes
1353        q = p;
1354        D = mapinto( Dp );
1355        cof= mapinto (cof);
1356        cog= mapinto (cog);
1357        d_deg=dp_deg;
1358        test= 0;
1359        equal= false;
1360      }
1361      else
1362      {
1363        //printf(" was bad, try more\n");
1364      }
1365      //else dp_deg > d_deg: bad prime
1366    }
1367    if ( i >= 0 )
1368    {
1369      TIMING_START (chinrem_reconstruction);
1370      Dn= Farey(D,q);
1371      cofn= Farey(cof,q);
1372      cogn= Farey(cog,q);
1373      TIMING_END_AND_PRINT (chinrem_reconstruction,
1374                           "time for rational reconstruction in modular gcd: ");
1375      int is_rat= isOn (SW_RATIONAL);
1376      On (SW_RATIONAL);
1377      cd = bCommonDen( Dn ); // we need On(SW_RATIONAL)
1378      cofn *= bCommonDen (cofn);
1379      cogn *= bCommonDen (cogn);
1380      if (!is_rat)
1381        Off (SW_RATIONAL);
1382      Dn *=cd;
1383      if (test != Dn)
1384        test= Dn;
1385      else
1386        equal= true;
1387      //Dn /=vcontent(Dn,Variable(1));
1388      TIMING_START (chinrem_termination);
1389      if ((terminationTest (f,g, cofn, cogn, Dn)) ||
1390          ((equal || q > b) && fdivides (Dn, f) && fdivides (Dn, g)))
1391      {
1392        TIMING_END_AND_PRINT (chinrem_termination,
1393                            "time for successful termination in modular gcd: ");
1394        //printf(" -> success\n");
1395        return Dn*gcdcfcg;
1396      }
1397      TIMING_END_AND_PRINT (chinrem_termination,
1398                          "time for unsuccessful termination in modular gcd: ");
1399      equal= false;
1400      //else: try more primes
1401    }
1402    else
1403    { // try other method
1404      //printf("try other gcd\n");
1405      Off(SW_USE_CHINREM_GCD);
1406      D=gcd_poly( f, g );
1407      On(SW_USE_CHINREM_GCD);
1408      return D*gcdcfcg;
1409    }
1410  }
1411}
1412
1413
Note: See TracBrowser for help on using the repository browser.