source: git/factory/cf_gcd.cc @ 297e92

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