source: git/factory/cf_gcd.cc @ 2711ab

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