source: git/factory/cf_gcd.cc @ ef20c7

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