source: git/factory/cf_gcd.cc @ b5c084

jengelh-datetimespielwiese
Last change on this file since b5c084 was b5c084, checked in by Martin Lee <martinlee84@…>, 12 years ago
minor fixes in gcd git-svn-id: file:///usr/local/Singular/svn/trunk@14369 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 33.5 KB
Line 
1/* emacs edit mode for this file is -*- C++ -*- */
2/* $Id$ */
3
4#include <config.h>
5
6#include "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 "ffreval.h"
19#include "algext.h"
20#include "fieldGCD.h"
21#include "cf_gcd_smallp.h"
22#include "cf_map_ext.h"
23#include "cf_util.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 bool gcd_avoid_mtaildegree ( CanonicalForm &, CanonicalForm &, CanonicalForm & );
35static void cf_prepgcd( const CanonicalForm &, const CanonicalForm &, int &, int &, int & );
36
37void out_cf(const char *s1,const CanonicalForm &f,const char *s2);
38
39CanonicalForm chinrem_gcd(const CanonicalForm & FF,const CanonicalForm & GG);
40CanonicalForm newGCD(CanonicalForm A, CanonicalForm B);
41
42bool
43gcd_test_one ( const CanonicalForm & f, const CanonicalForm & g, bool swap )
44{
45    int count = 0;
46    // assume polys have same level;
47
48    Variable v= Variable (1);
49    bool algExtension= (hasFirstAlgVar (f, v) || hasFirstAlgVar (g, v));
50    CanonicalForm lcf, lcg;
51    if ( swap )
52    {
53        lcf = swapvar( LC( f ), Variable(1), f.mvar() );
54        lcg = swapvar( LC( g ), Variable(1), f.mvar() );
55    }
56    else
57    {
58        lcf = LC( f, Variable(1) );
59        lcg = LC( g, Variable(1) );
60    }
61
62    CanonicalForm F, G;
63    if ( swap )
64    {
65        F=swapvar( f, Variable(1), f.mvar() );
66        G=swapvar( g, Variable(1), g.mvar() );
67    }
68    else
69    {
70        F = f;
71        G = g;
72    }
73
74    #define TEST_ONE_MAX 50
75    int p= getCharacteristic();
76    bool passToGF= false;
77    int k= 1;
78    if (p > 0 && p < TEST_ONE_MAX && CFFactory::gettype() != GaloisFieldDomain && !algExtension)
79    {
80      if (p == 2)
81        setCharacteristic (2, 6, 'Z');
82      else if (p == 3)
83        setCharacteristic (3, 4, 'Z');
84      else if (p == 5 || p == 7)
85        setCharacteristic (p, 3, 'Z');
86      else
87        setCharacteristic (p, 2, 'Z');
88      passToGF= true;
89    }
90    else if (p > 0 && CFFactory::gettype() == GaloisFieldDomain && ipower (p , getGFDegree()) < TEST_ONE_MAX)
91    {
92      k= getGFDegree();
93      if (ipower (p, 2*k) > TEST_ONE_MAX)
94        setCharacteristic (p, 2*k, gf_name);
95      else
96        setCharacteristic (p, 3*k, gf_name);
97      F= GFMapUp (F, k);
98      G= GFMapUp (G, k);
99      lcf= GFMapUp (lcf, k);
100      lcg= GFMapUp (lcg, k);
101    }
102    else if (p > 0 && p < TEST_ONE_MAX && algExtension)
103    {
104      bool extOfExt= false;
105      int d= degree (getMipo (v));
106      CFList source, dest;
107      Variable v2;
108      CanonicalForm primElem, imPrimElem;
109      if (p == 2 && d < 6)
110      {
111        zz_p::init (p);
112        bool primFail= false;
113        Variable vBuf;
114        primElem= primitiveElement (v, vBuf, primFail);
115        ASSERT (!primFail, "failure in integer factorizer");
116        if (d < 3)
117        {
118          zz_pX NTLIrredpoly;
119          BuildIrred (NTLIrredpoly, d*3);
120          CanonicalForm newMipo= convertNTLzzpX2CF (NTLIrredpoly, Variable (1));
121          v2= rootOf (newMipo);
122        }
123        else
124        {
125          zz_pX NTLIrredpoly;
126          BuildIrred (NTLIrredpoly, d*2);
127          CanonicalForm newMipo= convertNTLzzpX2CF (NTLIrredpoly, Variable (1));
128          v2= rootOf (newMipo);
129        }
130        imPrimElem= mapPrimElem (primElem, v, v2);
131        extOfExt= true;
132      }
133      else if ((p == 3 && d < 4) || ((p == 5 || p == 7) && d < 3))
134      {
135        zz_p::init (p);
136        bool primFail= false;
137        Variable vBuf;
138        primElem= primitiveElement (v, vBuf, primFail);
139        ASSERT (!primFail, "failure in integer factorizer");
140        zz_pX NTLIrredpoly;
141        BuildIrred (NTLIrredpoly, d*2);
142        CanonicalForm newMipo= convertNTLzzpX2CF (NTLIrredpoly, Variable (1));
143        v2= rootOf (newMipo);
144        imPrimElem= mapPrimElem (primElem, v, v2);
145        extOfExt= true;
146      }
147      if (extOfExt)
148      {
149        F= mapUp (F, v, v2, primElem, imPrimElem, source, dest);
150        G= mapUp (G, v, v2, primElem, imPrimElem, source, dest);
151        lcf= mapUp (lcf, v, v2, primElem, imPrimElem, source, dest);
152        lcg= mapUp (lcg, v, v2, primElem, imPrimElem, source, dest);
153        v= v2;
154      }
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        //out_cf("F:",f,"\n");
537        //out_cf("G:",g,"\n");
538        //out_cf("newGCD:",newGCD(f,g),"\n");
539        if (isOn(SW_USE_GCD_P) && (getCharacteristic()>0))
540        {
541          return newGCD(f,g);
542        }
543        if ( gcd_test_one( pi1, pi, true ) )
544        {
545          C=abs(C);
546          //out_cf("GCD:",C,"\n");
547          return C;
548        }
549        bpure = false;
550    }
551    else
552    {
553        bpure = isPurePoly(pi) && isPurePoly(pi1);
554#ifdef HAVE_NTL
555        if ( isOn(SW_USE_NTL_GCD_P) && bpure && (CFFactory::gettype() != GaloisFieldDomain))
556            return gcd_univar_ntlp(pi, pi1 ) * C;
557#endif
558    }
559    Variable v = f.mvar();
560    Hi = power( LC( pi1, v ), delta );
561    if ( (delta+1) % 2 )
562        bi = 1;
563    else
564        bi = -1;
565    int maxNumVars= tmax (getNumVars (pi), getNumVars (pi1));
566    CanonicalForm oldPi= pi, oldPi1= pi1;
567    while ( degree( pi1, v ) > 0 )
568    {
569        if (!(pi.isUnivariate() && pi1.isUnivariate()))
570        {
571          if (size (pi)/maxNumVars > 500 || size (pi1)/maxNumVars > 500)
572          {
573            On (SW_USE_FF_MOD_GCD);
574            C *= gcd (oldPi, oldPi1);
575            Off (SW_USE_FF_MOD_GCD);
576            return C;
577          }
578        }
579        pi2 = psr( pi, pi1, v );
580        pi2 = pi2 / bi;
581        pi = pi1; pi1 = pi2;
582        maxNumVars= tmax (getNumVars (pi), getNumVars (pi1));
583        if ( degree( pi1, v ) > 0 )
584        {
585            delta = degree( pi, v ) - degree( pi1, v );
586            if ( (delta+1) % 2 )
587                bi = LC( pi, v ) * power( Hi, delta );
588            else
589                bi = -LC( pi, v ) * power( Hi, delta );
590            Hi = power( LC( pi1, v ), delta ) / power( Hi, delta-1 );
591        }
592    }
593    if ( degree( pi1, v ) == 0 )
594    {
595      C=abs(C);
596      //out_cf("GCD:",C,"\n");
597      return C;
598    }
599    pi /= content( pi );
600    if ( bpure )
601        pi /= pi.lc();
602    C=abs(C*pi);
603    //out_cf("GCD:",C,"\n");
604    return C;
605}
606
607static CanonicalForm
608gcd_poly_0( const CanonicalForm & f, const CanonicalForm & g )
609{
610    CanonicalForm pi, pi1;
611    CanonicalForm C, Ci, Ci1, Hi, bi, pi2;
612    int delta = degree( f ) - degree( g );
613
614    if ( delta >= 0 )
615    {
616        pi = f; pi1 = g;
617    }
618    else
619    {
620        pi = g; pi1 = f; delta = -delta;
621    }
622    Ci = content( pi ); Ci1 = content( pi1 );
623    pi1 = pi1 / Ci1; pi = pi / Ci;
624    C = gcd( Ci, Ci1 );
625    if ( pi.isUnivariate() && pi1.isUnivariate() )
626    {
627#ifdef HAVE_NTL
628        if ( isOn(SW_USE_NTL_GCD_0) && isPurePoly(pi) && isPurePoly(pi1) )
629            return gcd_univar_ntl0(pi, pi1 ) * C;
630#endif
631        return gcd_poly_univar0( pi, pi1, true ) * C;
632    }
633    else if ( gcd_test_one( pi1, pi, true ) )
634      return C;
635    Variable v = f.mvar();
636    Hi = power( LC( pi1, v ), delta );
637    if ( (delta+1) % 2 )
638        bi = 1;
639    else
640        bi = -1;
641    while ( degree( pi1, v ) > 0 )
642    {
643        pi2 = psr( pi, pi1, v );
644        pi2 = pi2 / bi;
645        pi = pi1; pi1 = pi2;
646        if ( degree( pi1, v ) > 0 )
647        {
648            delta = degree( pi, v ) - degree( pi1, v );
649            if ( (delta+1) % 2 )
650                bi = LC( pi, v ) * power( Hi, delta );
651            else
652                bi = -LC( pi, v ) * power( Hi, delta );
653            Hi = power( LC( pi1, v ), delta ) / power( Hi, delta-1 );
654        }
655    }
656    if ( degree( pi1, v ) == 0 )
657        return C;
658    else
659        return C * pp( pi );
660}
661
662//{{{ CanonicalForm gcd_poly ( const CanonicalForm & f, const CanonicalForm & g )
663//{{{ docu
664//
665// gcd_poly() - calculate polynomial gcd.
666//
667// This is the dispatcher for polynomial gcd calculation.  We call either
668// ezgcd(), sparsemod() or gcd_poly1() in dependecy on the current
669// characteristic and settings of SW_USE_EZGCD and SW_USE_SPARSEMOD, resp.
670//
671// Used by gcd() and gcd_poly_univar0().
672//
673//}}}
674#if 0
675int si_factor_reminder=1;
676#endif
677CanonicalForm gcd_poly ( const CanonicalForm & f, const CanonicalForm & g )
678{
679  CanonicalForm fc, gc, d1;
680  int mp, cc, p1, pe;
681  mp = f.level()+1;
682  bool fc_isUnivariate=f.isUnivariate();
683  bool gc_isUnivariate=g.isUnivariate();
684  bool fc_and_gc_Univariate=fc_isUnivariate && gc_isUnivariate;
685#if 1
686  if (( getCharacteristic() == 0 )
687  && (f.level() >4)
688  && (g.level() >4)
689  && isOn( SW_USE_CHINREM_GCD)
690  && (!fc_and_gc_Univariate)
691  && (isPurePoly_m(f))
692  && (isPurePoly_m(g))
693  )
694  {
695      return chinrem_gcd( f, g );
696  }
697#endif
698  cf_prepgcd( f, g, cc, p1, pe);
699  if ( cc != 0 )
700  {
701    if ( cc > 0 )
702    {
703      fc = replacevar( f, Variable(cc), Variable(mp) );
704      gc = g;
705    }
706    else
707    {
708      fc = replacevar( g, Variable(-cc), Variable(mp) );
709      gc = f;
710    }
711    return cf_content( fc, gc );
712  }
713// now each appearing variable is in f and g
714  fc = f;
715  gc = g;
716  if( gcd_avoid_mtaildegree ( fc, gc, d1 ) )
717      return d1;
718  if ( getCharacteristic() != 0 )
719  {
720    if ((!fc_and_gc_Univariate)
721    && isOn(SW_USE_fieldGCD)
722    && (getCharacteristic() >100))
723    {
724      return fieldGCD(f,g);
725    }
726    #ifdef HAVE_NTL
727    else if ((!fc_and_gc_Univariate) && (isOn( SW_USE_EZGCD_P )))
728    {
729      /*if ( pe == 1 )
730        fc = fin_ezgcd( fc, gc );
731      else if ( pe > 0 )// no variable at position 1
732      {
733        fc = replacevar( fc, Variable(pe), Variable(1) );
734        gc = replacevar( gc, Variable(pe), Variable(1) );
735        fc = replacevar( fin_ezgcd( fc, gc ), Variable(1), Variable(pe) );
736      }
737      else
738      {
739        pe = -pe;
740        fc = swapvar( fc, Variable(pe), Variable(1) );
741        gc = swapvar( gc, Variable(pe), Variable(1) );
742        fc = swapvar( fin_ezgcd( fc, gc ), Variable(1), Variable(pe) );
743      }*/
744      fc= EZGCD_P (fc, gc);
745    }
746    #endif
747    else if (isOn(SW_USE_GCD_P))
748    {
749      fc=newGCD(fc,gc);
750    }
751    #ifdef HAVE_NTL
752    else if (isOn(SW_USE_FF_MOD_GCD) && !fc_and_gc_Univariate)
753    {
754      Variable a;
755      if (hasFirstAlgVar (fc, a) || hasFirstAlgVar (gc, a))
756        fc=GCD_Fp_extension (fc, gc, a);
757      else if (CFFactory::gettype() == GaloisFieldDomain)
758        fc=GCD_GF (fc, gc);
759      else
760        fc=GCD_small_p (fc, gc);
761    }
762    #endif
763    else if ( p1 == fc.level() )
764      fc = gcd_poly_p( fc, gc );
765    else
766    {
767      fc = replacevar( fc, Variable(p1), Variable(mp) );
768      gc = replacevar( gc, Variable(p1), Variable(mp) );
769      fc = replacevar( gcd_poly_p( fc, gc ), Variable(mp), Variable(p1) );
770    }
771  }
772  else if (!fc_and_gc_Univariate)
773  {
774    if (
775    isOn(SW_USE_CHINREM_GCD)
776    && (gc.level() >5)
777    && (fc.level() >5)
778    && (isPurePoly_m(fc)) && (isPurePoly_m(gc))
779    )
780    {
781    #if 0
782      if ( p1 == fc.level() )
783        fc = chinrem_gcd( fc, gc );
784      else
785      {
786        fc = replacevar( fc, Variable(p1), Variable(mp) );
787        gc = replacevar( gc, Variable(p1), Variable(mp) );
788        fc = replacevar( chinrem_gcd( fc, gc ), Variable(mp), Variable(p1) );
789      }
790    #else
791      fc = chinrem_gcd( fc, gc);
792    #endif
793    }
794    else if ( isOn( SW_USE_EZGCD ) )
795    {
796      if ( pe == 1 )
797        fc = ezgcd( fc, gc );
798      else if ( pe > 0 )// no variable at position 1
799      {
800        fc = replacevar( fc, Variable(pe), Variable(1) );
801        gc = replacevar( gc, Variable(pe), Variable(1) );
802        fc = replacevar( ezgcd( fc, gc ), Variable(1), Variable(pe) );
803      }
804      else
805      {
806        pe = -pe;
807        fc = swapvar( fc, Variable(pe), Variable(1) );
808        gc = swapvar( gc, Variable(pe), Variable(1) );
809        fc = swapvar( ezgcd( fc, gc ), Variable(1), Variable(pe) );
810      }
811    }
812    else if (
813    isOn(SW_USE_CHINREM_GCD)
814    && (isPurePoly_m(fc)) && (isPurePoly_m(gc))
815    )
816    {
817    #if 0
818      if ( p1 == fc.level() )
819        fc = chinrem_gcd( fc, gc );
820      else
821      {
822        fc = replacevar( fc, Variable(p1), Variable(mp) );
823        gc = replacevar( gc, Variable(p1), Variable(mp) );
824        fc = replacevar( chinrem_gcd( fc, gc ), Variable(mp), Variable(p1) );
825      }
826    #else
827      fc = chinrem_gcd( fc, gc);
828    #endif
829    }
830    else
831    {
832       fc = gcd_poly_0( fc, gc );
833    }
834  }
835  else
836  {
837    fc = gcd_poly_0( fc, gc );
838  }
839  if ( d1.degree() > 0 )
840    fc *= d1;
841  return fc;
842}
843//}}}
844
845//{{{ static CanonicalForm cf_content ( const CanonicalForm & f, const CanonicalForm & g )
846//{{{ docu
847//
848// cf_content() - return gcd(g, content(f)).
849//
850// content(f) is calculated with respect to f's main variable.
851//
852// Used by gcd(), content(), content( CF, Variable ).
853//
854//}}}
855static CanonicalForm
856cf_content ( const CanonicalForm & f, const CanonicalForm & g )
857{
858    if ( f.inPolyDomain() || ( f.inExtension() && ! getReduce( f.mvar() ) ) )
859    {
860        CFIterator i = f;
861        CanonicalForm result = g;
862        while ( i.hasTerms() && ! result.isOne() )
863        {
864            result = gcd( i.coeff(), result );
865            i++;
866        }
867        return result;
868    }
869    else
870        return abs( f );
871}
872//}}}
873
874//{{{ CanonicalForm content ( const CanonicalForm & f )
875//{{{ docu
876//
877// content() - return content(f) with respect to main variable.
878//
879// Normalizes result.
880//
881//}}}
882CanonicalForm
883content ( const CanonicalForm & f )
884{
885    if ( f.inPolyDomain() || ( f.inExtension() && ! getReduce( f.mvar() ) ) )
886    {
887        CFIterator i = f;
888        CanonicalForm result = abs( i.coeff() );
889        i++;
890        while ( i.hasTerms() && ! result.isOne() )
891        {
892            result = gcd( i.coeff(), result );
893            i++;
894        }
895        return result;
896    }
897    else
898        return abs( f );
899}
900//}}}
901
902//{{{ CanonicalForm content ( const CanonicalForm & f, const Variable & x )
903//{{{ docu
904//
905// content() - return content(f) with respect to x.
906//
907// x should be a polynomial variable.
908//
909//}}}
910CanonicalForm
911content ( const CanonicalForm & f, const Variable & x )
912{
913    ASSERT( x.level() > 0, "cannot calculate content with respect to algebraic variable" );
914    Variable y = f.mvar();
915
916    if ( y == x )
917        return cf_content( f, 0 );
918    else  if ( y < x )
919        return f;
920    else
921        return swapvar( content( swapvar( f, y, x ), y ), y, x );
922}
923//}}}
924
925//{{{ CanonicalForm vcontent ( const CanonicalForm & f, const Variable & x )
926//{{{ docu
927//
928// vcontent() - return content of f with repect to variables >= x.
929//
930// The content is recursively calculated over all coefficients in
931// f having level less than x.  x should be a polynomial
932// variable.
933//
934//}}}
935CanonicalForm
936vcontent ( const CanonicalForm & f, const Variable & x )
937{
938    ASSERT( x.level() > 0, "cannot calculate vcontent with respect to algebraic variable" );
939
940    if ( f.mvar() <= x )
941        return content( f, x );
942    else {
943        CFIterator i;
944        CanonicalForm d = 0;
945        for ( i = f; i.hasTerms() && ! d.isOne(); i++ )
946            d = gcd( d, vcontent( i.coeff(), x ) );
947        return d;
948    }
949}
950//}}}
951
952//{{{ CanonicalForm pp ( const CanonicalForm & f )
953//{{{ docu
954//
955// pp() - return primitive part of f.
956//
957// Returns zero if f equals zero, otherwise f / content(f).
958//
959//}}}
960CanonicalForm
961pp ( const CanonicalForm & f )
962{
963    if ( f.isZero() )
964        return f;
965    else
966        return f / content( f );
967}
968//}}}
969
970CanonicalForm
971gcd ( const CanonicalForm & f, const CanonicalForm & g )
972{
973    bool b = f.isZero();
974    if ( b || g.isZero() )
975    {
976        if ( b )
977            return abs( g );
978        else
979            return abs( f );
980    }
981    if ( f.inPolyDomain() || g.inPolyDomain() )
982    {
983        if ( f.mvar() != g.mvar() )
984        {
985            if ( f.mvar() > g.mvar() )
986                return cf_content( f, g );
987            else
988                return cf_content( g, f );
989        }
990        if (isOn(SW_USE_QGCD))
991        {
992          Variable m;
993          if (
994          (getCharacteristic() == 0) &&
995          (hasFirstAlgVar(f,m) || hasFirstAlgVar(g,m))
996          )
997          {
998            bool on_rational = isOn(SW_RATIONAL);
999            CanonicalForm r=QGCD(f,g);
1000            On(SW_RATIONAL);
1001            CanonicalForm cdF = bCommonDen( r );
1002            if (!on_rational) Off(SW_RATIONAL);
1003            return cdF*r;
1004          }
1005        }
1006
1007        if ( f.inExtension() && getReduce( f.mvar() ) )
1008            return CanonicalForm(1);
1009        else
1010        {
1011            if ( fdivides( f, g ) )
1012                return abs( f );
1013            else  if ( fdivides( g, f ) )
1014                return abs( g );
1015            if ( !( getCharacteristic() == 0 && isOn( SW_RATIONAL ) ) )
1016            {
1017                CanonicalForm d;
1018                d = gcd_poly( f, g );
1019                return abs( d );
1020            }
1021            else
1022            {
1023                CanonicalForm cdF = bCommonDen( f );
1024                CanonicalForm cdG = bCommonDen( g );
1025                Off( SW_RATIONAL );
1026                CanonicalForm l = lcm( cdF, cdG );
1027                On( SW_RATIONAL );
1028                CanonicalForm F = f * l, G = g * l;
1029                Off( SW_RATIONAL );
1030                l = gcd_poly( F, G );
1031                On( SW_RATIONAL );
1032                return abs( l );
1033            }
1034        }
1035    }
1036    if ( f.inBaseDomain() && g.inBaseDomain() )
1037        return bgcd( f, g );
1038    else
1039        return 1;
1040}
1041
1042//{{{ CanonicalForm lcm ( const CanonicalForm & f, const CanonicalForm & g )
1043//{{{ docu
1044//
1045// lcm() - return least common multiple of f and g.
1046//
1047// The lcm is calculated using the formula lcm(f, g) = f * g / gcd(f, g).
1048//
1049// Returns zero if one of f or g equals zero.
1050//
1051//}}}
1052CanonicalForm
1053lcm ( const CanonicalForm & f, const CanonicalForm & g )
1054{
1055    if ( f.isZero() || g.isZero() )
1056        return 0;
1057    else
1058        return ( f / gcd( f, g ) ) * g;
1059}
1060//}}}
1061
1062#ifdef HAVE_NTL
1063
1064static CanonicalForm
1065gcd_univar_ntl0( const CanonicalForm & F, const CanonicalForm & G )
1066{
1067    ZZX F1=convertFacCF2NTLZZX(F);
1068    ZZX G1=convertFacCF2NTLZZX(G);
1069    ZZX R=GCD(F1,G1);
1070    return convertNTLZZX2CF(R,F.mvar());
1071}
1072
1073static CanonicalForm
1074gcd_univar_ntlp( const CanonicalForm & F, const CanonicalForm & G )
1075{
1076  if (fac_NTL_char!=getCharacteristic())
1077  {
1078    fac_NTL_char=getCharacteristic();
1079    #ifdef NTL_ZZ
1080    ZZ r;
1081    r=getCharacteristic();
1082    ZZ_pContext ccc(r);
1083    #else
1084    zz_pContext ccc(getCharacteristic());
1085    #endif
1086    ccc.restore();
1087    #ifdef NTL_ZZ
1088    ZZ_p::init(r);
1089    #else
1090    zz_p::init(getCharacteristic());
1091    #endif
1092  }
1093  #ifdef NTL_ZZ
1094  ZZ_pX F1=convertFacCF2NTLZZpX(F);
1095  ZZ_pX G1=convertFacCF2NTLZZpX(G);
1096  ZZ_pX R=GCD(F1,G1);
1097  return  convertNTLZZpX2CF(R,F.mvar());
1098  #else
1099  zz_pX F1=convertFacCF2NTLzzpX(F);
1100  zz_pX G1=convertFacCF2NTLzzpX(G);
1101  zz_pX R=GCD(F1,G1);
1102  return  convertNTLzzpX2CF(R,F.mvar());
1103  #endif
1104}
1105
1106#endif
1107
1108static bool
1109gcd_avoid_mtaildegree ( CanonicalForm & f1, CanonicalForm & g1, CanonicalForm & d1 )
1110{
1111    bool rdy = true;
1112    int df = f1.taildegree();
1113    int dg = g1.taildegree();
1114
1115    d1 = d1.genOne();
1116    if ( dg == 0 )
1117    {
1118        if ( df == 0 )
1119            return false;
1120        else
1121        {
1122            if ( f1.degree() == df )
1123                d1 = cf_content( g1, LC( f1 ) );
1124            else
1125            {
1126                f1 /= power( f1.mvar(), df );
1127                rdy = false;
1128            }
1129        }
1130    }
1131    else
1132    {
1133        if ( df == 0)
1134        {
1135            if ( g1.degree() == dg )
1136                d1 = cf_content( f1, LC( g1 ) );
1137            else
1138            {
1139                g1 /= power( g1.mvar(), dg );
1140                rdy = false;
1141            }
1142        }
1143        else
1144        {
1145            if ( df > dg )
1146                d1 = power( f1.mvar(), dg );
1147            else
1148                d1 = power( f1.mvar(), df );
1149            if ( f1.degree() == df )
1150            {
1151                if (g1.degree() == dg)
1152                    d1 *= gcd( LC( f1 ), LC( g1 ) );
1153                else
1154                {
1155                    g1 /= power( g1.mvar(), dg);
1156                    d1 *= cf_content( g1, LC( f1 ) );
1157                }
1158            }
1159            else
1160            {
1161                f1 /= power( f1.mvar(), df );
1162                if ( g1.degree() == dg )
1163                    d1 *= cf_content( f1, LC( g1 ) );
1164                else
1165                {
1166                    g1 /= power( g1.mvar(), dg );
1167                    rdy = false;
1168                }
1169            }
1170        }
1171    }
1172    return rdy;
1173}
1174
1175/*
1176*  compute positions p1 and pe of optimal variables:
1177*    pe is used in "ezgcd" and
1178*    p1 in "gcd_poly1"
1179*/
1180static
1181void optvalues ( const int * df, const int * dg, const int n, int & p1, int &pe )
1182{
1183    int i, o1, oe;
1184    if ( df[n] > dg[n] )
1185    {
1186        o1 = df[n]; oe = dg[n];
1187    }
1188    else
1189    {
1190        o1 = dg[n]; oe = df[n];
1191    }
1192    i = n-1;
1193    while ( i > 0 )
1194    {
1195        if ( df[i] != 0 )
1196        {
1197            if ( df[i] > dg[i] )
1198            {
1199                if ( o1 > df[i]) { o1 = df[i]; p1 = i; }
1200                if ( oe <= dg[i]) { oe = dg[i]; pe = i; }
1201            }
1202            else
1203            {
1204                if ( o1 > dg[i]) { o1 = dg[i]; p1 = i; }
1205                if ( oe <= df[i]) { oe = df[i]; pe = i; }
1206            }
1207        }
1208        i--;
1209    }
1210}
1211
1212/*
1213*  make some changes of variables, see optvalues
1214*/
1215static void
1216cf_prepgcd( const CanonicalForm & f, const CanonicalForm & g, int & cc, int & p1, int &pe )
1217{
1218    int i, k, n;
1219    n = f.level();
1220    cc = 0;
1221    p1 = pe = n;
1222    if ( n == 1 )
1223        return;
1224    int * degsf = new int[n+1];
1225    int * degsg = new int[n+1];
1226    for ( i = n; i > 0; i-- )
1227    {
1228        degsf[i] = degsg[i] = 0;
1229    }
1230    degsf = degrees( f, degsf );
1231    degsg = degrees( g, degsg );
1232
1233    k = 0;
1234    for ( i = n-1; i > 0; i-- )
1235    {
1236        if ( degsf[i] == 0 )
1237        {
1238            if ( degsg[i] != 0 )
1239            {
1240                cc = -i;
1241                break;
1242            }
1243        }
1244        else
1245        {
1246            if ( degsg[i] == 0 )
1247            {
1248                cc = i;
1249                break;
1250            }
1251            else k++;
1252        }
1253    }
1254
1255    if ( ( cc == 0 ) && ( k != 0 ) )
1256        optvalues( degsf, degsg, n, p1, pe );
1257    if ( ( pe != 1 ) && ( degsf[1] != 0 ) )
1258        pe = -pe;
1259
1260    delete [] degsf;
1261    delete [] degsg;
1262}
1263
1264
1265static CanonicalForm
1266balance_p ( const CanonicalForm & f, const CanonicalForm & q )
1267{
1268    Variable x = f.mvar();
1269    CanonicalForm result = 0, qh = q / 2;
1270    CanonicalForm c;
1271    CFIterator i;
1272    for ( i = f; i.hasTerms(); i++ )
1273    {
1274        c = i.coeff();
1275        if ( c.inCoeffDomain())
1276        {
1277          if ( c > qh )
1278            result += power( x, i.exp() ) * (c - q);
1279          else
1280            result += power( x, i.exp() ) * c;
1281        }
1282        else
1283          result += power( x, i.exp() ) * balance_p(c,q);
1284    }
1285    return result;
1286}
1287
1288CanonicalForm chinrem_gcd ( const CanonicalForm & FF, const CanonicalForm & GG )
1289{
1290  CanonicalForm f, g, cg, cl, q(0), Dp, newD, D, newq;
1291  int p, i, dp_deg, d_deg=-1;
1292
1293  CanonicalForm cd ( bCommonDen( FF ));
1294  f=cd*FF;
1295  f /=vcontent(f,Variable(1));
1296  //cd = bCommonDen( f ); f *=cd;
1297  //f /=vcontent(f,Variable(1));
1298
1299  cd = bCommonDen( GG );
1300  g=cd*GG;
1301  g /=vcontent(g,Variable(1));
1302  //cd = bCommonDen( g ); g *=cd;
1303  //g /=vcontent(g,Variable(1));
1304
1305  i = cf_getNumBigPrimes() - 1;
1306  cl =  f.lc()* g.lc();
1307
1308  while ( true )
1309  {
1310    p = cf_getBigPrime( i );
1311    i--;
1312    while ( i >= 0 && mod( cl, p ) == 0 )
1313    {
1314      p = cf_getBigPrime( i );
1315      i--;
1316    }
1317    //printf("try p=%d\n",p);
1318    setCharacteristic( p );
1319    Dp = gcd_poly( mapinto( f ), mapinto( g ) );
1320    Dp /=Dp.lc();
1321    setCharacteristic( 0 );
1322    dp_deg=totaldegree(Dp);
1323    if ( dp_deg == 0 )
1324    {
1325      //printf(" -> 1\n");
1326      return CanonicalForm(1);
1327    }
1328    if ( q.isZero() )
1329    {
1330      D = mapinto( Dp );
1331      d_deg=dp_deg;
1332      q = p;
1333    }
1334    else
1335    {
1336      if ( dp_deg == d_deg )
1337      {
1338        chineseRemainder( D, q, mapinto( Dp ), p, newD, newq );
1339        q = newq;
1340        D = newD;
1341      }
1342      else if ( dp_deg < d_deg )
1343      {
1344        //printf(" were all bad, try more\n");
1345        // all previous p's are bad primes
1346        q = p;
1347        D = mapinto( Dp );
1348        d_deg=dp_deg;
1349      }
1350      else
1351      {
1352        //printf(" was bad, try more\n");
1353      }
1354      //else dp_deg > d_deg: bad prime
1355    }
1356    if ( i >= 0 )
1357    {
1358      CanonicalForm Dn= Farey(D,q);
1359      int is_rat=isOn(SW_RATIONAL);
1360      On(SW_RATIONAL);
1361      CanonicalForm cd = bCommonDen( Dn ); // we need On(SW_RATIONAL)
1362      if (!is_rat) Off(SW_RATIONAL);
1363      Dn *=cd;
1364      //Dn /=vcontent(Dn,Variable(1));
1365      if ( fdivides( Dn, f ) && fdivides( Dn, g ) )
1366      {
1367        //printf(" -> success\n");
1368        return Dn;
1369      }
1370      //else: try more primes
1371    }
1372    else
1373    { // try other method
1374      //printf("try other gcd\n");
1375      Off(SW_USE_CHINREM_GCD);
1376      D=gcd_poly( f, g );
1377      On(SW_USE_CHINREM_GCD);
1378      return D;
1379    }
1380  }
1381}
1382
Note: See TracBrowser for help on using the repository browser.