source: git/factory/cf_gcd.cc @ f45bfb

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