source: git/factory/cf_gcd.cc @ c4682e0

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