source: git/factory/cf_gcd.cc @ 276c3f

spielwiese
Last change on this file since 276c3f was c1b9927, checked in by Hans Schoenemann <hannes@…>, 13 years ago
- removed some unsed variables - never put static inline routine without a body in a .h file git-svn-id: file:///usr/local/Singular/svn/trunk@14265 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 33.9 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      {
757        fc=GCD_Fp_extension (fc, gc, a);
758      }
759      if (CFFactory::gettype() == GaloisFieldDomain)
760      {
761        fc=GCD_GF (fc, gc);
762      }
763      fc=GCD_small_p (fc, gc);
764    }
765    #endif
766    else if ( p1 == fc.level() )
767      fc = gcd_poly_p( fc, gc );
768    else
769    {
770      fc = replacevar( fc, Variable(p1), Variable(mp) );
771      gc = replacevar( gc, Variable(p1), Variable(mp) );
772      fc = replacevar( gcd_poly_p( fc, gc ), Variable(mp), Variable(p1) );
773    }
774  }
775  else if (!fc_and_gc_Univariate)
776  {
777    if (
778    isOn(SW_USE_CHINREM_GCD)
779    && (gc.level() >5)
780    && (fc.level() >5)
781    && (isPurePoly_m(fc)) && (isPurePoly_m(gc))
782    )
783    {
784    #if 0
785      if ( p1 == fc.level() )
786        fc = chinrem_gcd( fc, gc );
787      else
788      {
789        fc = replacevar( fc, Variable(p1), Variable(mp) );
790        gc = replacevar( gc, Variable(p1), Variable(mp) );
791        fc = replacevar( chinrem_gcd( fc, gc ), Variable(mp), Variable(p1) );
792      }
793    #else
794      fc = chinrem_gcd( fc, gc);
795    #endif
796    }
797    else if ( isOn( SW_USE_EZGCD ) )
798    {
799      if ( pe == 1 )
800        fc = ezgcd( fc, gc );
801      else if ( pe > 0 )// no variable at position 1
802      {
803        fc = replacevar( fc, Variable(pe), Variable(1) );
804        gc = replacevar( gc, Variable(pe), Variable(1) );
805        fc = replacevar( ezgcd( fc, gc ), Variable(1), Variable(pe) );
806      }
807      else
808      {
809        pe = -pe;
810        fc = swapvar( fc, Variable(pe), Variable(1) );
811        gc = swapvar( gc, Variable(pe), Variable(1) );
812        fc = swapvar( ezgcd( fc, gc ), Variable(1), Variable(pe) );
813      }
814    }
815    else if (
816    isOn(SW_USE_CHINREM_GCD)
817    && (isPurePoly_m(fc)) && (isPurePoly_m(gc))
818    )
819    {
820    #if 0
821      if ( p1 == fc.level() )
822        fc = chinrem_gcd( fc, gc );
823      else
824      {
825        fc = replacevar( fc, Variable(p1), Variable(mp) );
826        gc = replacevar( gc, Variable(p1), Variable(mp) );
827        fc = replacevar( chinrem_gcd( fc, gc ), Variable(mp), Variable(p1) );
828      }
829    #else
830      fc = chinrem_gcd( fc, gc);
831    #endif
832    }
833    else
834    {
835       fc = gcd_poly_0( fc, gc );
836    }
837  }
838  else
839  {
840    fc = gcd_poly_0( fc, gc );
841  }
842  if ( d1.degree() > 0 )
843    fc *= d1;
844  return fc;
845}
846//}}}
847
848//{{{ static CanonicalForm cf_content ( const CanonicalForm & f, const CanonicalForm & g )
849//{{{ docu
850//
851// cf_content() - return gcd(g, content(f)).
852//
853// content(f) is calculated with respect to f's main variable.
854//
855// Used by gcd(), content(), content( CF, Variable ).
856//
857//}}}
858static CanonicalForm
859cf_content ( const CanonicalForm & f, const CanonicalForm & g )
860{
861    if ( f.inPolyDomain() || ( f.inExtension() && ! getReduce( f.mvar() ) ) )
862    {
863        CFIterator i = f;
864        CanonicalForm result = g;
865        while ( i.hasTerms() && ! result.isOne() )
866        {
867            result = gcd( i.coeff(), result );
868            i++;
869        }
870        return result;
871    }
872    else
873        return abs( f );
874}
875//}}}
876
877//{{{ CanonicalForm content ( const CanonicalForm & f )
878//{{{ docu
879//
880// content() - return content(f) with respect to main variable.
881//
882// Normalizes result.
883//
884//}}}
885CanonicalForm
886content ( const CanonicalForm & f )
887{
888    if ( f.inPolyDomain() || ( f.inExtension() && ! getReduce( f.mvar() ) ) )
889    {
890        CFIterator i = f;
891        CanonicalForm result = abs( i.coeff() );
892        i++;
893        while ( i.hasTerms() && ! result.isOne() )
894        {
895            result = gcd( i.coeff(), result );
896            i++;
897        }
898        return result;
899    }
900    else
901        return abs( f );
902}
903//}}}
904
905//{{{ CanonicalForm content ( const CanonicalForm & f, const Variable & x )
906//{{{ docu
907//
908// content() - return content(f) with respect to x.
909//
910// x should be a polynomial variable.
911//
912//}}}
913CanonicalForm
914content ( const CanonicalForm & f, const Variable & x )
915{
916    ASSERT( x.level() > 0, "cannot calculate content with respect to algebraic variable" );
917    Variable y = f.mvar();
918
919    if ( y == x )
920        return cf_content( f, 0 );
921    else  if ( y < x )
922        return f;
923    else
924        return swapvar( content( swapvar( f, y, x ), y ), y, x );
925}
926//}}}
927
928//{{{ CanonicalForm vcontent ( const CanonicalForm & f, const Variable & x )
929//{{{ docu
930//
931// vcontent() - return content of f with repect to variables >= x.
932//
933// The content is recursively calculated over all coefficients in
934// f having level less than x.  x should be a polynomial
935// variable.
936//
937//}}}
938CanonicalForm
939vcontent ( const CanonicalForm & f, const Variable & x )
940{
941    ASSERT( x.level() > 0, "cannot calculate vcontent with respect to algebraic variable" );
942
943    if ( f.mvar() <= x )
944        return content( f, x );
945    else {
946        CFIterator i;
947        CanonicalForm d = 0;
948        for ( i = f; i.hasTerms() && ! d.isOne(); i++ )
949            d = gcd( d, vcontent( i.coeff(), x ) );
950        return d;
951    }
952}
953//}}}
954
955//{{{ CanonicalForm pp ( const CanonicalForm & f )
956//{{{ docu
957//
958// pp() - return primitive part of f.
959//
960// Returns zero if f equals zero, otherwise f / content(f).
961//
962//}}}
963CanonicalForm
964pp ( const CanonicalForm & f )
965{
966    if ( f.isZero() )
967        return f;
968    else
969        return f / content( f );
970}
971//}}}
972
973CanonicalForm
974gcd ( const CanonicalForm & f, const CanonicalForm & g )
975{
976    bool b = f.isZero();
977    if ( b || g.isZero() )
978    {
979        if ( b )
980            return abs( g );
981        else
982            return abs( f );
983    }
984    if ( f.inPolyDomain() || g.inPolyDomain() )
985    {
986        if ( f.mvar() != g.mvar() )
987        {
988            if ( f.mvar() > g.mvar() )
989                return cf_content( f, g );
990            else
991                return cf_content( g, f );
992        }
993        if (isOn(SW_USE_QGCD))
994        {
995          Variable m;
996          if (
997          (getCharacteristic() == 0) &&
998          (hasFirstAlgVar(f,m) || hasFirstAlgVar(g,m))
999          )
1000          {
1001            bool on_rational = isOn(SW_RATIONAL);
1002            CanonicalForm r=QGCD(f,g);
1003            On(SW_RATIONAL);
1004            CanonicalForm cdF = bCommonDen( r );
1005            if (!on_rational) Off(SW_RATIONAL);
1006            return cdF*r;
1007          }
1008        }
1009
1010        if ( f.inExtension() && getReduce( f.mvar() ) )
1011            return CanonicalForm(1);
1012        else
1013        {
1014            if ( fdivides( f, g ) )
1015                return abs( f );
1016            else  if ( fdivides( g, f ) )
1017                return abs( g );
1018            if ( !( getCharacteristic() == 0 && isOn( SW_RATIONAL ) ) )
1019            {
1020                CanonicalForm d;
1021#if 1
1022                do{ d = gcd_poly( f, g ); }
1023                while ((!fdivides(d,f)) || (!fdivides(d,g)));
1024#else
1025                while(1)
1026                {
1027                  d = gcd_poly( f, g );
1028                  if ((fdivides(d,f)) && (fdivides(d,g))) break;
1029                  printf("g"); fflush(stdout);
1030                }
1031#endif
1032                return abs( d );
1033            }
1034            else
1035            {
1036                CanonicalForm cdF = bCommonDen( f );
1037                CanonicalForm cdG = bCommonDen( g );
1038                Off( SW_RATIONAL );
1039                CanonicalForm l = lcm( cdF, cdG );
1040                On( SW_RATIONAL );
1041                CanonicalForm F = f * l, G = g * l;
1042                Off( SW_RATIONAL );
1043                do { l = gcd_poly( F, G ); }
1044                while ((!fdivides(l,F)) || (!fdivides(l,G)));
1045                On( SW_RATIONAL );
1046                return abs( l );
1047            }
1048        }
1049    }
1050    if ( f.inBaseDomain() && g.inBaseDomain() )
1051        return bgcd( f, g );
1052    else
1053        return 1;
1054}
1055
1056//{{{ CanonicalForm lcm ( const CanonicalForm & f, const CanonicalForm & g )
1057//{{{ docu
1058//
1059// lcm() - return least common multiple of f and g.
1060//
1061// The lcm is calculated using the formula lcm(f, g) = f * g / gcd(f, g).
1062//
1063// Returns zero if one of f or g equals zero.
1064//
1065//}}}
1066CanonicalForm
1067lcm ( const CanonicalForm & f, const CanonicalForm & g )
1068{
1069    if ( f.isZero() || g.isZero() )
1070        return 0;
1071    else
1072        return ( f / gcd( f, g ) ) * g;
1073}
1074//}}}
1075
1076#ifdef HAVE_NTL
1077
1078static CanonicalForm
1079gcd_univar_ntl0( const CanonicalForm & F, const CanonicalForm & G )
1080{
1081    ZZX F1=convertFacCF2NTLZZX(F);
1082    ZZX G1=convertFacCF2NTLZZX(G);
1083    ZZX R=GCD(F1,G1);
1084    return convertNTLZZX2CF(R,F.mvar());
1085}
1086
1087static CanonicalForm
1088gcd_univar_ntlp( const CanonicalForm & F, const CanonicalForm & G )
1089{
1090  if (fac_NTL_char!=getCharacteristic())
1091  {
1092    fac_NTL_char=getCharacteristic();
1093    #ifdef NTL_ZZ
1094    ZZ r;
1095    r=getCharacteristic();
1096    ZZ_pContext ccc(r);
1097    #else
1098    zz_pContext ccc(getCharacteristic());
1099    #endif
1100    ccc.restore();
1101    #ifdef NTL_ZZ
1102    ZZ_p::init(r);
1103    #else
1104    zz_p::init(getCharacteristic());
1105    #endif
1106  }
1107  #ifdef NTL_ZZ
1108  ZZ_pX F1=convertFacCF2NTLZZpX(F);
1109  ZZ_pX G1=convertFacCF2NTLZZpX(G);
1110  ZZ_pX R=GCD(F1,G1);
1111  return  convertNTLZZpX2CF(R,F.mvar());
1112  #else
1113  zz_pX F1=convertFacCF2NTLzzpX(F);
1114  zz_pX G1=convertFacCF2NTLzzpX(G);
1115  zz_pX R=GCD(F1,G1);
1116  return  convertNTLzzpX2CF(R,F.mvar());
1117  #endif
1118}
1119
1120#endif
1121
1122static bool
1123gcd_avoid_mtaildegree ( CanonicalForm & f1, CanonicalForm & g1, CanonicalForm & d1 )
1124{
1125    bool rdy = true;
1126    int df = f1.taildegree();
1127    int dg = g1.taildegree();
1128
1129    d1 = d1.genOne();
1130    if ( dg == 0 )
1131    {
1132        if ( df == 0 )
1133            return false;
1134        else
1135        {
1136            if ( f1.degree() == df )
1137                d1 = cf_content( g1, LC( f1 ) );
1138            else
1139            {
1140                f1 /= power( f1.mvar(), df );
1141                rdy = false;
1142            }
1143        }
1144    }
1145    else
1146    {
1147        if ( df == 0)
1148        {
1149            if ( g1.degree() == dg )
1150                d1 = cf_content( f1, LC( g1 ) );
1151            else
1152            {
1153                g1 /= power( g1.mvar(), dg );
1154                rdy = false;
1155            }
1156        }
1157        else
1158        {
1159            if ( df > dg )
1160                d1 = power( f1.mvar(), dg );
1161            else
1162                d1 = power( f1.mvar(), df );
1163            if ( f1.degree() == df )
1164            {
1165                if (g1.degree() == dg)
1166                    d1 *= gcd( LC( f1 ), LC( g1 ) );
1167                else
1168                {
1169                    g1 /= power( g1.mvar(), dg);
1170                    d1 *= cf_content( g1, LC( f1 ) );
1171                }
1172            }
1173            else
1174            {
1175                f1 /= power( f1.mvar(), df );
1176                if ( g1.degree() == dg )
1177                    d1 *= cf_content( f1, LC( g1 ) );
1178                else
1179                {
1180                    g1 /= power( g1.mvar(), dg );
1181                    rdy = false;
1182                }
1183            }
1184        }
1185    }
1186    return rdy;
1187}
1188
1189/*
1190*  compute positions p1 and pe of optimal variables:
1191*    pe is used in "ezgcd" and
1192*    p1 in "gcd_poly1"
1193*/
1194static
1195void optvalues ( const int * df, const int * dg, const int n, int & p1, int &pe )
1196{
1197    int i, o1, oe;
1198    if ( df[n] > dg[n] )
1199    {
1200        o1 = df[n]; oe = dg[n];
1201    }
1202    else
1203    {
1204        o1 = dg[n]; oe = df[n];
1205    }
1206    i = n-1;
1207    while ( i > 0 )
1208    {
1209        if ( df[i] != 0 )
1210        {
1211            if ( df[i] > dg[i] )
1212            {
1213                if ( o1 > df[i]) { o1 = df[i]; p1 = i; }
1214                if ( oe <= dg[i]) { oe = dg[i]; pe = i; }
1215            }
1216            else
1217            {
1218                if ( o1 > dg[i]) { o1 = dg[i]; p1 = i; }
1219                if ( oe <= df[i]) { oe = df[i]; pe = i; }
1220            }
1221        }
1222        i--;
1223    }
1224}
1225
1226/*
1227*  make some changes of variables, see optvalues
1228*/
1229static void
1230cf_prepgcd( const CanonicalForm & f, const CanonicalForm & g, int & cc, int & p1, int &pe )
1231{
1232    int i, k, n;
1233    n = f.level();
1234    cc = 0;
1235    p1 = pe = n;
1236    if ( n == 1 )
1237        return;
1238    int * degsf = new int[n+1];
1239    int * degsg = new int[n+1];
1240    for ( i = n; i > 0; i-- )
1241    {
1242        degsf[i] = degsg[i] = 0;
1243    }
1244    degsf = degrees( f, degsf );
1245    degsg = degrees( g, degsg );
1246
1247    k = 0;
1248    for ( i = n-1; i > 0; i-- )
1249    {
1250        if ( degsf[i] == 0 )
1251        {
1252            if ( degsg[i] != 0 )
1253            {
1254                cc = -i;
1255                break;
1256            }
1257        }
1258        else
1259        {
1260            if ( degsg[i] == 0 )
1261            {
1262                cc = i;
1263                break;
1264            }
1265            else k++;
1266        }
1267    }
1268
1269    if ( ( cc == 0 ) && ( k != 0 ) )
1270        optvalues( degsf, degsg, n, p1, pe );
1271    if ( ( pe != 1 ) && ( degsf[1] != 0 ) )
1272        pe = -pe;
1273
1274    delete [] degsf;
1275    delete [] degsg;
1276}
1277
1278
1279static CanonicalForm
1280balance_p ( const CanonicalForm & f, const CanonicalForm & q )
1281{
1282    Variable x = f.mvar();
1283    CanonicalForm result = 0, qh = q / 2;
1284    CanonicalForm c;
1285    CFIterator i;
1286    for ( i = f; i.hasTerms(); i++ )
1287    {
1288        c = i.coeff();
1289        if ( c.inCoeffDomain())
1290        {
1291          if ( c > qh )
1292            result += power( x, i.exp() ) * (c - q);
1293          else
1294            result += power( x, i.exp() ) * c;
1295        }
1296        else
1297          result += power( x, i.exp() ) * balance_p(c,q);
1298    }
1299    return result;
1300}
1301
1302CanonicalForm chinrem_gcd ( const CanonicalForm & FF, const CanonicalForm & GG )
1303{
1304  CanonicalForm f, g, cg, cl, q(0), Dp, newD, D, newq;
1305  int p, i, dp_deg, d_deg=-1;
1306
1307  CanonicalForm cd ( bCommonDen( FF ));
1308  f=cd*FF;
1309  f /=vcontent(f,Variable(1));
1310  //cd = bCommonDen( f ); f *=cd;
1311  //f /=vcontent(f,Variable(1));
1312
1313  cd = bCommonDen( GG );
1314  g=cd*GG;
1315  g /=vcontent(g,Variable(1));
1316  //cd = bCommonDen( g ); g *=cd;
1317  //g /=vcontent(g,Variable(1));
1318
1319  i = cf_getNumBigPrimes() - 1;
1320  cl =  f.lc()* g.lc();
1321
1322  while ( true )
1323  {
1324    p = cf_getBigPrime( i );
1325    i--;
1326    while ( i >= 0 && mod( cl, p ) == 0 )
1327    {
1328      p = cf_getBigPrime( i );
1329      i--;
1330    }
1331    //printf("try p=%d\n",p);
1332    setCharacteristic( p );
1333    Dp = gcd_poly( mapinto( f ), mapinto( g ) );
1334    Dp /=Dp.lc();
1335    setCharacteristic( 0 );
1336    dp_deg=totaldegree(Dp);
1337    if ( dp_deg == 0 )
1338    {
1339      //printf(" -> 1\n");
1340      return CanonicalForm(1);
1341    }
1342    if ( q.isZero() )
1343    {
1344      D = mapinto( Dp );
1345      d_deg=dp_deg;
1346      q = p;
1347    }
1348    else
1349    {
1350      if ( dp_deg == d_deg )
1351      {
1352        chineseRemainder( D, q, mapinto( Dp ), p, newD, newq );
1353        q = newq;
1354        D = newD;
1355      }
1356      else if ( dp_deg < d_deg )
1357      {
1358        //printf(" were all bad, try more\n");
1359        // all previous p's are bad primes
1360        q = p;
1361        D = mapinto( Dp );
1362        d_deg=dp_deg;
1363      }
1364      else
1365      {
1366        //printf(" was bad, try more\n");
1367      }
1368      //else dp_deg > d_deg: bad prime
1369    }
1370    if ( i >= 0 )
1371    {
1372      CanonicalForm Dn= Farey(D,q);
1373      int is_rat=isOn(SW_RATIONAL);
1374      On(SW_RATIONAL);
1375      CanonicalForm cd = bCommonDen( Dn ); // we need On(SW_RATIONAL)
1376      if (!is_rat) Off(SW_RATIONAL);
1377      Dn *=cd;
1378      //Dn /=vcontent(Dn,Variable(1));
1379      if ( fdivides( Dn, f ) && fdivides( Dn, g ) )
1380      {
1381        //printf(" -> success\n");
1382        return Dn;
1383      }
1384      //else: try more primes
1385    }
1386    else
1387    { // try other method
1388      //printf("try other gcd\n");
1389      Off(SW_USE_CHINREM_GCD);
1390      D=gcd_poly( f, g );
1391      On(SW_USE_CHINREM_GCD);
1392      return D;
1393    }
1394  }
1395}
1396
Note: See TracBrowser for help on using the repository browser.