source: git/factory/cf_gcd.cc @ 09f10e

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