source: git/factory/cf_gcd.cc @ 6db552

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