source: git/factory/cf_gcd.cc @ c6eecb

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