source: git/factory/cf_gcd.cc @ f06059

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