source: git/factory/cf_gcd.cc @ 2072126

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