source: git/factory/cf_gcd.cc @ fc9f44

spielwiese
Last change on this file since fc9f44 was fc9f44, checked in by Hans Schoenemann <hannes@…>, 14 years ago
minor fixes to gcd stuff git-svn-id: file:///usr/local/Singular/svn/trunk@12760 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 29.8 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 "ftmpl_functions.h"
20#include "ffreval.h"
21#include "algext.h"
22#include "fieldGCD.h"
23#include "cf_gcd_smallp.h"
24
25
26#ifdef HAVE_NTL
27#include <NTL/ZZX.h>
28#include "NTLconvert.h"
29bool isPurePoly(const CanonicalForm & );
30static CanonicalForm gcd_univar_ntl0( const CanonicalForm &, const CanonicalForm & );
31static CanonicalForm gcd_univar_ntlp( const CanonicalForm &, const CanonicalForm & );
32#endif
33
34static CanonicalForm cf_content ( const CanonicalForm &, const CanonicalForm & );
35static bool gcd_avoid_mtaildegree ( CanonicalForm &, CanonicalForm &, CanonicalForm & );
36static void cf_prepgcd( const CanonicalForm &, const CanonicalForm &, int &, int &, int & );
37
38void out_cf(const char *s1,const CanonicalForm &f,const char *s2);
39
40CanonicalForm chinrem_gcd(const CanonicalForm & FF,const CanonicalForm & GG);
41CanonicalForm newGCD(CanonicalForm A, CanonicalForm B);
42
43bool
44gcd_test_one ( const CanonicalForm & f, const CanonicalForm & g, bool swap )
45{
46    int count = 0;
47    // assume polys have same level;
48    CFRandom * sample = CFRandomFactory::generate();
49    REvaluation e( 2, tmax( f.level(), g.level() ), *sample );
50    delete sample;
51    CanonicalForm lcf, lcg;
52    if ( swap )
53    {
54        lcf = swapvar( LC( f ), Variable(1), f.mvar() );
55        lcg = swapvar( LC( g ), Variable(1), f.mvar() );
56    }
57    else
58    {
59        lcf = LC( f, Variable(1) );
60        lcg = LC( g, Variable(1) );
61    }
62    #define TEST_ONE_MAX 50
63    while ( ( e( lcf ).isZero() || e( lcg ).isZero() ) && count < TEST_ONE_MAX )
64    {
65        e.nextpoint();
66        count++;
67    }
68    if ( count == TEST_ONE_MAX )
69        return false;
70    CanonicalForm F, G;
71    if ( swap )
72    {
73        F=swapvar( f, Variable(1), f.mvar() );
74        G=swapvar( g, Variable(1), g.mvar() );
75    }
76    else
77    {
78        F = f;
79        G = g;
80    }
81    if ( e(F).taildegree() > 0 && e(G).taildegree() > 0 )
82        return false;
83    return gcd( e( F ), e( G ) ).degree() < 1;
84}
85
86//{{{ static CanonicalForm icontent ( const CanonicalForm & f, const CanonicalForm & c )
87//{{{ docu
88//
89// icontent() - return gcd of c and all coefficients of f which
90//   are in a coefficient domain.
91//
92// Used by icontent().
93//
94//}}}
95static CanonicalForm
96icontent ( const CanonicalForm & f, const CanonicalForm & c )
97{
98    if ( f.inBaseDomain() )
99    {
100      if (c.isZero()) return abs(f);
101      return bgcd( f, c );
102    }
103    //else if ( f.inCoeffDomain() )
104    //   return gcd(f,c);
105    else
106    {
107        CanonicalForm g = c;
108        for ( CFIterator i = f; i.hasTerms() && ! g.isOne(); i++ )
109            g = icontent( i.coeff(), g );
110        return g;
111    }
112}
113//}}}
114
115//{{{ CanonicalForm icontent ( const CanonicalForm & f )
116//{{{ docu
117//
118// icontent() - return gcd over all coefficients of f which are
119//   in a coefficient domain.
120//
121//}}}
122CanonicalForm
123icontent ( const CanonicalForm & f )
124{
125    return icontent( f, 0 );
126}
127//}}}
128
129//{{{ CanonicalForm extgcd ( const CanonicalForm & f, const CanonicalForm & g, CanonicalForm & a, CanonicalForm & b )
130//{{{ docu
131//
132// extgcd() - returns polynomial extended gcd of f and g.
133//
134// Returns gcd(f, g) and a and b sucht that f*a+g*b=gcd(f, g).
135// The gcd is calculated using an extended euclidean polynomial
136// remainder sequence, so f and g should be polynomials over an
137// euclidean domain.  Normalizes result.
138//
139// Note: be sure that f and g have the same level!
140//
141//}}}
142CanonicalForm
143extgcd ( const CanonicalForm & f, const CanonicalForm & g, CanonicalForm & a, CanonicalForm & b )
144{
145#ifdef HAVE_NTL
146  if (isOn(SW_USE_NTL_GCD_P) && ( getCharacteristic() > 0 )
147  &&  (f.level()==g.level()) && isPurePoly(f) && isPurePoly(g))
148  {
149    if (fac_NTL_char!=getCharacteristic())
150    {
151      fac_NTL_char=getCharacteristic();
152      #ifdef NTL_ZZ
153      ZZ r;
154      r=getCharacteristic();
155      ZZ_pContext ccc(r);
156      #else
157      zz_pContext ccc(getCharacteristic());
158      #endif
159      ccc.restore();
160      #ifdef NTL_ZZ
161      ZZ_p::init(r);
162      #else
163      zz_p::init(getCharacteristic());
164      #endif
165    }
166    #ifdef NTL_ZZ
167    ZZ_pX F1=convertFacCF2NTLZZpX(f);
168    ZZ_pX G1=convertFacCF2NTLZZpX(g);
169    ZZ_pX R;
170    ZZ_pX A,B;
171    XGCD(R,A,B,F1,G1);
172    a=convertNTLZZpX2CF(A,f.mvar());
173    b=convertNTLZZpX2CF(B,f.mvar());
174    return convertNTLZZpX2CF(R,f.mvar());
175    #else
176    zz_pX F1=convertFacCF2NTLzzpX(f);
177    zz_pX G1=convertFacCF2NTLzzpX(g);
178    zz_pX R;
179    zz_pX A,B;
180    XGCD(R,A,B,F1,G1);
181    a=convertNTLzzpX2CF(A,f.mvar());
182    b=convertNTLzzpX2CF(B,f.mvar());
183    return convertNTLzzpX2CF(R,f.mvar());
184    #endif
185  }
186  if (isOn(SW_USE_NTL_GCD_0) && ( getCharacteristic() ==0)
187  && (f.level()==g.level()) && isPurePoly(f) && isPurePoly(g))
188  {
189    CanonicalForm fc=bCommonDen(f);
190    CanonicalForm gc=bCommonDen(g);
191    ZZX F1=convertFacCF2NTLZZX(f*fc);
192    ZZX G1=convertFacCF2NTLZZX(g*gc);
193    ZZX R=GCD(F1,G1);
194    CanonicalForm r=convertNTLZZX2CF(R,f.mvar());
195    ZZ RR;
196    ZZX A,B;
197    if (r.inCoeffDomain())
198    {
199      XGCD(RR,A,B,F1,G1,1);
200      CanonicalForm rr=convertZZ2CF(RR);
201      ASSERT (!rr.isZero(), "NTL:XGCD failed");
202      a=convertNTLZZX2CF(A,f.mvar())*fc/rr;
203      b=convertNTLZZX2CF(B,f.mvar())*gc/rr;
204      return CanonicalForm(1);
205    }
206    else
207    {
208      fc=bCommonDen(f);
209      gc=bCommonDen(g);
210      F1=convertFacCF2NTLZZX(f*fc/r);
211      G1=convertFacCF2NTLZZX(g*gc/r);
212      XGCD(RR,A,B,F1,G1,1);
213      a=convertNTLZZX2CF(A,f.mvar())*fc;
214      b=convertNTLZZX2CF(B,f.mvar())*gc;
215      CanonicalForm rr=convertZZ2CF(RR);
216      ASSERT (!rr.isZero(), "NTL:XGCD failed");
217      r*=rr;
218      if ( r.sign() < 0 ) { r= -r; a= -a; b= -b; }
219      return r;
220    }
221  }
222#endif
223  // may contain bug in the co-factors, see track 107
224  CanonicalForm contf = content( f );
225  CanonicalForm contg = content( g );
226
227  CanonicalForm p0 = f / contf, p1 = g / contg;
228  CanonicalForm f0 = 1, f1 = 0, g0 = 0, g1 = 1, q, r;
229
230  while ( ! p1.isZero() )
231  {
232      divrem( p0, p1, q, r );
233      p0 = p1; p1 = r;
234      r = g0 - g1 * q;
235      g0 = g1; g1 = r;
236      r = f0 - f1 * q;
237      f0 = f1; f1 = r;
238  }
239  CanonicalForm contp0 = content( p0 );
240  a = f0 / ( contf * contp0 );
241  b = g0 / ( contg * contp0 );
242  p0 /= contp0;
243  if ( p0.sign() < 0 )
244  {
245      p0 = -p0;
246      a = -a;
247      b = -b;
248  }
249  return p0;
250}
251//}}}
252
253//{{{ static CanonicalForm balance ( const CanonicalForm & f, const CanonicalForm & q )
254//{{{ docu
255//
256// balance() - map f from positive to symmetric representation
257//   mod q.
258//
259// This makes sense for univariate polynomials over Z only.
260// q should be an integer.
261//
262// Used by gcd_poly_univar0().
263//
264//}}}
265static CanonicalForm
266balance ( const CanonicalForm & f, const CanonicalForm & q )
267{
268    Variable x = f.mvar();
269    CanonicalForm result = 0, qh = q / 2;
270    CanonicalForm c;
271    CFIterator i;
272    for ( i = f; i.hasTerms(); i++ ) {
273        c = mod( i.coeff(), q );
274        if ( c > qh )
275            result += power( x, i.exp() ) * (c - q);
276        else
277            result += power( x, i.exp() ) * c;
278    }
279    return result;
280}
281//}}}
282
283static CanonicalForm gcd_poly_univar0( const CanonicalForm & F, const CanonicalForm & G, bool primitive )
284{
285  CanonicalForm f, g, c, cg, cl, BB, B, M, q, Dp, newD, D, newq;
286  int p, i;
287
288  if ( primitive )
289  {
290    f = F;
291    g = G;
292    c = 1;
293  }
294  else
295  {
296    CanonicalForm cF = content( F ), cG = content( G );
297    f = F / cF;
298    g = G / cG;
299    c = bgcd( cF, cG );
300  }
301  cg = gcd( f.lc(), g.lc() );
302  cl = ( f.lc() / cg ) * g.lc();
303//     B = 2 * cg * tmin(
304//         maxnorm(f)*power(CanonicalForm(2),f.degree())*isqrt(f.degree()+1),
305//         maxnorm(g)*power(CanonicalForm(2),g.degree())*isqrt(g.degree()+1)
306//         )+1;
307  M = tmin( maxNorm(f), maxNorm(g) );
308  BB = power(CanonicalForm(2),tmin(f.degree(),g.degree()))*M;
309  q = 0;
310  i = cf_getNumSmallPrimes() - 1;
311  while ( true )
312  {
313    B = BB;
314    while ( i >= 0 && q < B )
315    {
316      p = cf_getSmallPrime( i );
317      i--;
318      while ( i >= 0 && mod( cl, p ) == 0 )
319      {
320        p = cf_getSmallPrime( i );
321        i--;
322      }
323      setCharacteristic( p );
324      Dp = gcd( mapinto( f ), mapinto( g ) );
325      Dp = ( Dp / Dp.lc() ) * mapinto( cg );
326      setCharacteristic( 0 );
327      if ( Dp.degree() == 0 )
328        return c;
329      if ( q.isZero() )
330      {
331        D = mapinto( Dp );
332        q = p;
333        B = power(CanonicalForm(2),D.degree())*M+1;
334      }
335      else
336      {
337        if ( Dp.degree() == D.degree() )
338        {
339          chineseRemainder( D, q, mapinto( Dp ), p, newD, newq );
340          q = newq;
341          D = newD;
342        }
343        else if ( Dp.degree() < D.degree() )
344        {
345          // all previous p's are bad primes
346          q = p;
347          D = mapinto( Dp );
348          B = power(CanonicalForm(2),D.degree())*M+1;
349        }
350        // else p is a bad prime
351      }
352    }
353    if ( i >= 0 )
354    {
355      // now balance D mod q
356      D = pp( balance( D, q ) );
357      if ( fdivides( D, f ) && fdivides( D, g ) )
358        return D * c;
359      else
360        q = 0;
361    }
362    else
363      return gcd_poly( F, G );
364    DEBOUTLN( cerr, "another try ..." );
365  }
366}
367
368static CanonicalForm
369gcd_poly_p( const CanonicalForm & f, const CanonicalForm & g )
370{
371    CanonicalForm pi, pi1;
372    CanonicalForm C, Ci, Ci1, Hi, bi, pi2;
373    bool bpure;
374    int delta = degree( f ) - degree( g );
375
376    if ( delta >= 0 )
377    {
378        pi = f; pi1 = g;
379    }
380    else
381    {
382        pi = g; pi1 = f; delta = -delta;
383    }
384    Ci = content( pi ); Ci1 = content( pi1 );
385    pi1 = pi1 / Ci1; pi = pi / Ci;
386    C = gcd( Ci, Ci1 );
387    if ( !( pi.isUnivariate() && pi1.isUnivariate() ) )
388    {
389        //out_cf("F:",f,"\n");
390        //out_cf("G:",g,"\n");
391        //out_cf("newGCD:",newGCD(f,g),"\n");
392        if (isOn(SW_USE_GCD_P) && (getCharacteristic()>0))
393        {
394          return newGCD(f,g);
395        }
396        if ( gcd_test_one( pi1, pi, true ) )
397        {
398          C=abs(C);
399          //out_cf("GCD:",C,"\n");
400          return C;
401        }
402        bpure = false;
403    }
404    else
405    {
406        bpure = isPurePoly(pi) && isPurePoly(pi1);
407#ifdef HAVE_NTL
408        if ( isOn(SW_USE_NTL_GCD_P) && bpure )
409            return gcd_univar_ntlp(pi, pi1 ) * C;
410#endif
411    }
412    Variable v = f.mvar();
413    Hi = power( LC( pi1, v ), delta );
414    if ( (delta+1) % 2 )
415        bi = 1;
416    else
417        bi = -1;
418    while ( degree( pi1, v ) > 0 )
419    {
420        pi2 = psr( pi, pi1, v );
421        pi2 = pi2 / bi;
422        pi = pi1; pi1 = pi2;
423        if ( degree( pi1, v ) > 0 )
424        {
425            delta = degree( pi, v ) - degree( pi1, v );
426            if ( (delta+1) % 2 )
427                bi = LC( pi, v ) * power( Hi, delta );
428            else
429                bi = -LC( pi, v ) * power( Hi, delta );
430            Hi = power( LC( pi1, v ), delta ) / power( Hi, delta-1 );
431        }
432    }
433    if ( degree( pi1, v ) == 0 )
434    {
435      C=abs(C);
436      //out_cf("GCD:",C,"\n");
437      return C;
438    }
439    pi /= content( pi );
440    if ( bpure )
441        pi /= pi.lc();
442    C=abs(C*pi);
443    //out_cf("GCD:",C,"\n");
444    return C;
445}
446
447static CanonicalForm
448gcd_poly_0( const CanonicalForm & f, const CanonicalForm & g )
449{
450    CanonicalForm pi, pi1;
451    CanonicalForm C, Ci, Ci1, Hi, bi, pi2;
452    int delta = degree( f ) - degree( g );
453
454    if ( delta >= 0 )
455    {
456        pi = f; pi1 = g;
457    }
458    else
459    {
460        pi = g; pi1 = f; delta = -delta;
461    }
462    Ci = content( pi ); Ci1 = content( pi1 );
463    pi1 = pi1 / Ci1; pi = pi / Ci;
464    C = gcd( Ci, Ci1 );
465    if ( pi.isUnivariate() && pi1.isUnivariate() )
466    {
467#ifdef HAVE_NTL
468        if ( isOn(SW_USE_NTL_GCD_0) && isPurePoly(pi) && isPurePoly(pi1) )
469            return gcd_univar_ntl0(pi, pi1 ) * C;
470#endif
471        return gcd_poly_univar0( pi, pi1, true ) * C;
472    }
473    else if ( gcd_test_one( pi1, pi, true ) )
474      return C;
475    Variable v = f.mvar();
476    Hi = power( LC( pi1, v ), delta );
477    if ( (delta+1) % 2 )
478        bi = 1;
479    else
480        bi = -1;
481    while ( degree( pi1, v ) > 0 )
482    {
483        pi2 = psr( pi, pi1, v );
484        pi2 = pi2 / bi;
485        pi = pi1; pi1 = pi2;
486        if ( degree( pi1, v ) > 0 )
487        {
488            delta = degree( pi, v ) - degree( pi1, v );
489            if ( (delta+1) % 2 )
490                bi = LC( pi, v ) * power( Hi, delta );
491            else
492                bi = -LC( pi, v ) * power( Hi, delta );
493            Hi = power( LC( pi1, v ), delta ) / power( Hi, delta-1 );
494        }
495    }
496    if ( degree( pi1, v ) == 0 )
497        return C;
498    else
499        return C * pp( pi );
500}
501
502//{{{ CanonicalForm gcd_poly ( const CanonicalForm & f, const CanonicalForm & g )
503//{{{ docu
504//
505// gcd_poly() - calculate polynomial gcd.
506//
507// This is the dispatcher for polynomial gcd calculation.  We call either
508// ezgcd(), sparsemod() or gcd_poly1() in dependecy on the current
509// characteristic and settings of SW_USE_EZGCD and SW_USE_SPARSEMOD, resp.
510//
511// Used by gcd() and gcd_poly_univar0().
512//
513//}}}
514#if 0
515int si_factor_reminder=1;
516#endif
517CanonicalForm gcd_poly ( const CanonicalForm & f, const CanonicalForm & g )
518{
519  CanonicalForm fc, gc, d1;
520  int mp, cc, p1, pe;
521  mp = f.level()+1;
522  bool fc_isUnivariate=f.isUnivariate();
523  bool gc_isUnivariate=g.isUnivariate();
524  bool fc_and_gc_Univariate=fc_isUnivariate && gc_isUnivariate;
525#if 1
526  if (( getCharacteristic() == 0 )
527  && (f.level() >4)
528  && (g.level() >4)
529  && isOn( SW_USE_CHINREM_GCD)
530  && (!fc_and_gc_Univariate)
531  && (isPurePoly_m(f))
532  && (isPurePoly_m(g))
533  )
534  {
535      return chinrem_gcd( f, g );
536  }
537#endif
538  cf_prepgcd( f, g, cc, p1, pe);
539  if ( cc != 0 )
540  {
541    if ( cc > 0 )
542    {
543      fc = replacevar( f, Variable(cc), Variable(mp) );
544      gc = g;
545    }
546    else
547    {
548      fc = replacevar( g, Variable(-cc), Variable(mp) );
549      gc = f;
550    }
551    return cf_content( fc, gc );
552  }
553// now each appearing variable is in f and g
554  fc = f;
555  gc = g;
556  if( gcd_avoid_mtaildegree ( fc, gc, d1 ) )
557      return d1;
558  if ( getCharacteristic() != 0 )
559  {
560    if (isOn(SW_USE_fieldGCD)
561    && (!fc_and_gc_Univariate)
562    && (getCharacteristic() >100))
563    {
564      return fieldGCD(f,g);
565    }
566    else if (isOn( SW_USE_EZGCD_P ) && (!fc_and_gc_Univariate))
567    {
568      if ( pe == 1 )
569        fc = fin_ezgcd( fc, gc );
570      else if ( pe > 0 )// no variable at position 1
571      {
572        fc = replacevar( fc, Variable(pe), Variable(1) );
573        gc = replacevar( gc, Variable(pe), Variable(1) );
574        fc = replacevar( fin_ezgcd( fc, gc ), Variable(1), Variable(pe) );
575      }
576      else
577      {
578        pe = -pe;
579        fc = swapvar( fc, Variable(pe), Variable(1) );
580        gc = swapvar( gc, Variable(pe), Variable(1) );
581        fc = swapvar( fin_ezgcd( fc, gc ), Variable(1), Variable(pe) );
582      }
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    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.