source: git/factory/cf_gcd.cc @ 7cb5590

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