source: git/factory/cf_gcd.cc @ 5944b4

spielwiese
Last change on this file since 5944b4 was 5944b4, checked in by Hans Schönemann <hannes@…>, 17 years ago
*hannes: debug possibility git-svn-id: file:///usr/local/Singular/svn/trunk@9425 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 21.3 KB
Line 
1/* emacs edit mode for this file is -*- C++ -*- */
2/* $Id: cf_gcd.cc,v 1.48 2006-09-25 12:33:11 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 ( fdivides( D, f ) && fdivides( 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    bool bpure;
289    int delta = degree( f ) - degree( g );
290
291    if ( delta >= 0 )
292    {
293        pi = f; pi1 = g;
294    }
295    else
296    {
297        pi = g; pi1 = f; delta = -delta;
298    }
299    Ci = content( pi ); Ci1 = content( pi1 );
300    pi1 = pi1 / Ci1; pi = pi / Ci;
301    C = gcd( Ci, Ci1 );
302    if ( !( pi.isUnivariate() && pi1.isUnivariate() ) )
303    {
304        if ( gcd_test_one( pi1, pi, true ) )
305            return C;
306        bpure = false;
307     
308    }
309    else
310    {
311        bpure = isPurePoly(pi) && isPurePoly(pi1);
312#ifdef HAVE_NTL
313        if ( isOn(SW_USE_NTL_GCD_P) && bpure )
314            return gcd_univar_ntlp(pi, pi1 ) * C;
315#endif
316    }
317    Variable v = f.mvar();
318    Hi = power( LC( pi1, v ), delta );
319    if ( (delta+1) % 2 )
320        bi = 1;
321    else
322        bi = -1;
323    while ( degree( pi1, v ) > 0 ) {
324        pi2 = psr( pi, pi1, v );
325        pi2 = pi2 / bi;
326        pi = pi1; pi1 = pi2;
327        if ( degree( pi1, v ) > 0 ) {
328            delta = degree( pi, v ) - degree( pi1, v );
329            if ( (delta+1) % 2 )
330                bi = LC( pi, v ) * power( Hi, delta );
331            else
332                bi = -LC( pi, v ) * power( Hi, delta );
333            Hi = power( LC( pi1, v ), delta ) / power( Hi, delta-1 );
334        }
335    }
336    if ( degree( pi1, v ) == 0 )
337        return C;
338    pi /= content( pi );
339    if ( bpure )
340        pi /= pi.lc();
341    return C * pi;
342}
343
344static CanonicalForm
345gcd_poly_0( const CanonicalForm & f, const CanonicalForm & g )
346{
347    CanonicalForm pi, pi1;
348    CanonicalForm C, Ci, Ci1, Hi, bi, pi2;
349    int delta = degree( f ) - degree( g );
350
351    if ( delta >= 0 )
352    {
353        pi = f; pi1 = g;
354    }
355    else
356    {
357        pi = g; pi1 = f; delta = -delta;
358    }
359    Ci = content( pi ); Ci1 = content( pi1 );
360    pi1 = pi1 / Ci1; pi = pi / Ci;
361    C = gcd( Ci, Ci1 );
362    if ( pi.isUnivariate() && pi1.isUnivariate() )
363    {
364#ifdef HAVE_NTL
365        if ( isOn(SW_USE_NTL_GCD_0) && isPurePoly(pi) && isPurePoly(pi1) )
366            return gcd_univar_ntl0(pi, pi1 ) * C;
367#endif
368        return gcd_poly_univar0( pi, pi1, true ) * C;
369    }
370    else if ( gcd_test_one( pi1, pi, true ) )
371      return C;
372    Variable v = f.mvar();
373    Hi = power( LC( pi1, v ), delta );
374    if ( (delta+1) % 2 )
375        bi = 1;
376    else
377        bi = -1;
378    while ( degree( pi1, v ) > 0 ) {
379        pi2 = psr( pi, pi1, v );
380        pi2 = pi2 / bi;
381        pi = pi1; pi1 = pi2;
382        if ( degree( pi1, v ) > 0 ) {
383            delta = degree( pi, v ) - degree( pi1, v );
384            if ( (delta+1) % 2 )
385                bi = LC( pi, v ) * power( Hi, delta );
386            else
387                bi = -LC( pi, v ) * power( Hi, delta );
388            Hi = power( LC( pi1, v ), delta ) / power( Hi, delta-1 );
389        }
390    }
391    if ( degree( pi1, v ) == 0 )
392        return C;
393    else
394        return C * pp( pi );
395}
396
397//{{{ static CanonicalForm gcd_poly ( const CanonicalForm & f, const CanonicalForm & g )
398//{{{ docu
399//
400// gcd_poly() - calculate polynomial gcd.
401//
402// This is the dispatcher for polynomial gcd calculation.  We call either
403// ezgcd(), sparsemod() or gcd_poly1() in dependecy on the current
404// characteristic and settings of SW_USE_EZGCD and SW_USE_SPARSEMOD, resp.
405//
406// Used by gcd() and gcd_poly_univar0().
407//
408//}}}
409#if 0
410int si_factor_reminder=1;
411#endif
412static CanonicalForm
413gcd_poly ( const CanonicalForm & f, const CanonicalForm & g )
414{
415    CanonicalForm fc, gc, d1;
416    int mp, cc, p1, pe;
417    mp = f.level()+1;
418    cf_prepgcd( f, g, cc, p1, pe);
419    if ( cc != 0 )
420    {
421        if ( cc > 0 )
422        {
423            fc = replacevar( f, Variable(cc), Variable(mp) );
424            gc = g;
425        }
426        else
427        {
428            fc = replacevar( g, Variable(-cc), Variable(mp) );
429            gc = f;
430        }
431        return cf_content( fc, gc );
432    }
433// now each appearing variable is in f and g
434    fc = f;
435    gc = g;
436    if( gcd_avoid_mtaildegree ( fc, gc, d1 ) )
437        return d1;
438    if ( getCharacteristic() != 0 )
439    {
440        if ( p1 == fc.level() )
441            fc = gcd_poly_p( fc, gc );
442        else
443        {
444            fc = replacevar( fc, Variable(p1), Variable(mp) );
445            gc = replacevar( gc, Variable(p1), Variable(mp) );
446            fc = replacevar( gcd_poly_p( fc, gc ), Variable(mp), Variable(p1) );
447        }
448    }
449    else if ( isOn( SW_USE_EZGCD ) && !f.isUnivariate() )
450    {
451        if ( pe == 1 )
452            fc = ezgcd( fc, gc );
453        else if ( pe > 0 )// no variable at position 1
454        {
455            fc = replacevar( fc, Variable(pe), Variable(1) );
456            gc = replacevar( gc, Variable(pe), Variable(1) );
457            fc = replacevar( ezgcd( fc, gc ), Variable(1), Variable(pe) );
458        }
459        else
460        {
461            pe = -pe;
462            fc = swapvar( fc, Variable(pe), Variable(1) );
463            gc = swapvar( gc, Variable(pe), Variable(1) );
464            fc = swapvar( ezgcd( fc, gc ), Variable(1), Variable(pe) );
465        }
466    }
467    else
468    {
469        fc = gcd_poly_0( fc, gc );
470    }
471    if ( d1.degree() > 0 )
472        fc *= d1;
473    return fc;
474}
475//}}}
476
477//{{{ static CanonicalForm cf_content ( const CanonicalForm & f, const CanonicalForm & g )
478//{{{ docu
479//
480// cf_content() - return gcd(g, content(f)).
481//
482// content(f) is calculated with respect to f's main variable.
483//
484// Used by gcd(), content(), content( CF, Variable ).
485//
486//}}}
487static CanonicalForm
488cf_content ( const CanonicalForm & f, const CanonicalForm & g )
489{
490    if ( f.inPolyDomain() || ( f.inExtension() && ! getReduce( f.mvar() ) ) ) {
491        CFIterator i = f;
492        CanonicalForm result = g;
493        while ( i.hasTerms() && ! result.isOne() ) {
494            result = gcd( i.coeff(), result );
495            i++;
496        }
497        return result;
498    }
499    else
500        return abs( f );
501}
502//}}}
503
504//{{{ CanonicalForm content ( const CanonicalForm & f )
505//{{{ docu
506//
507// content() - return content(f) with respect to main variable.
508//
509// Normalizes result.
510//
511//}}}
512CanonicalForm
513content ( const CanonicalForm & f )
514{
515    if ( f.inPolyDomain() || ( f.inExtension() && ! getReduce( f.mvar() ) ) ) {
516        CFIterator i = f;
517        CanonicalForm result = abs( i.coeff() );
518        i++;
519        while ( i.hasTerms() && ! result.isOne() ) {
520            result = gcd( i.coeff(), result );
521            i++;
522        }
523        return result;
524    }
525    else
526        return abs( f );
527}
528//}}}
529
530//{{{ CanonicalForm content ( const CanonicalForm & f, const Variable & x )
531//{{{ docu
532//
533// content() - return content(f) with respect to x.
534//
535// x should be a polynomial variable.
536//
537//}}}
538CanonicalForm
539content ( const CanonicalForm & f, const Variable & x )
540{
541    ASSERT( x.level() > 0, "cannot calculate content with respect to algebraic variable" );
542    Variable y = f.mvar();
543
544    if ( y == x )
545        return cf_content( f, 0 );
546    else  if ( y < x )
547        return f;
548    else
549        return swapvar( content( swapvar( f, y, x ), y ), y, x );
550}
551//}}}
552
553//{{{ CanonicalForm vcontent ( const CanonicalForm & f, const Variable & x )
554//{{{ docu
555//
556// vcontent() - return content of f with repect to variables >= x.
557//
558// The content is recursively calculated over all coefficients in
559// f having level less than x.  x should be a polynomial
560// variable.
561//
562//}}}
563CanonicalForm
564vcontent ( const CanonicalForm & f, const Variable & x )
565{
566    ASSERT( x.level() > 0, "cannot calculate vcontent with respect to algebraic variable" );
567
568    if ( f.mvar() <= x )
569        return content( f, x );
570    else {
571        CFIterator i;
572        CanonicalForm d = 0;
573        for ( i = f; i.hasTerms() && ! d.isOne(); i++ )
574            d = gcd( d, vcontent( i.coeff(), x ) );
575        return d;
576    }
577}
578//}}}
579
580//{{{ CanonicalForm pp ( const CanonicalForm & f )
581//{{{ docu
582//
583// pp() - return primitive part of f.
584//
585// Returns zero if f equals zero, otherwise f / content(f).
586//
587//}}}
588CanonicalForm
589pp ( const CanonicalForm & f )
590{
591    if ( f.isZero() )
592        return f;
593    else
594        return f / content( f );
595}
596//}}}
597
598CanonicalForm
599gcd ( const CanonicalForm & f, const CanonicalForm & g )
600{
601    bool b = f.isZero();
602    if ( b || g.isZero() )
603    {
604        if ( b )
605            return abs( g );
606        else
607            return abs( f );
608    }
609    if ( f.inPolyDomain() || g.inPolyDomain() )
610    {
611        if ( f.mvar() != g.mvar() )
612        {
613            if ( f.mvar() > g.mvar() )
614                return cf_content( f, g );
615            else
616                return cf_content( g, f );
617        }
618        if ( f.inExtension() && getReduce( f.mvar() ) )
619            return 1;
620        else
621        {
622            if ( fdivides( f, g ) )
623                return abs( f );
624            else  if ( fdivides( g, f ) )
625                return abs( g );
626            if ( !( getCharacteristic() == 0 && isOn( SW_RATIONAL ) ) )
627            {
628                CanonicalForm d;
629#if 1
630                do{ d = gcd_poly( f, g ); }
631                while ((!fdivides(d,f)) || (!fdivides(d,g)));
632#else
633                while(1)
634                { 
635                  d = gcd_poly( f, g ); 
636                  if ((fdivides(d,f)) && (fdivides(d,g))) break;
637                  printf("g"); fflush(stdout);
638                }
639#endif
640                return abs( d );
641            }
642            else
643            {
644                CanonicalForm cdF = bCommonDen( f );
645                CanonicalForm cdG = bCommonDen( g );
646                Off( SW_RATIONAL );
647                CanonicalForm l = lcm( cdF, cdG );
648                On( SW_RATIONAL );
649                CanonicalForm F = f * l, G = g * l;
650                Off( SW_RATIONAL );
651                do { l = gcd_poly( F, G ); }
652                while ((!fdivides(l,F)) || (!fdivides(l,G)));
653                On( SW_RATIONAL );
654                return abs( l );
655            }
656        }
657    }
658    if ( f.inBaseDomain() && g.inBaseDomain() )
659        return bgcd( f, g );
660    else
661        return 1;
662}
663
664//{{{ CanonicalForm lcm ( const CanonicalForm & f, const CanonicalForm & g )
665//{{{ docu
666//
667// lcm() - return least common multiple of f and g.
668//
669// The lcm is calculated using the formula lcm(f, g) = f * g / gcd(f, g).
670//
671// Returns zero if one of f or g equals zero.
672//
673//}}}
674CanonicalForm
675lcm ( const CanonicalForm & f, const CanonicalForm & g )
676{
677    if ( f.isZero() || g.isZero() )
678        return 0;
679    else
680        return ( f / gcd( f, g ) ) * g;
681}
682//}}}
683
684#ifdef HAVE_NTL
685
686static CanonicalForm
687gcd_univar_ntl0( const CanonicalForm & F, const CanonicalForm & G )
688{
689    ZZX F1=convertFacCF2NTLZZX(F);
690    ZZX G1=convertFacCF2NTLZZX(G);
691    ZZX R=GCD(F1,G1);
692    return convertNTLZZX2CF(R,F.mvar());
693}
694
695static CanonicalForm
696gcd_univar_ntlp( const CanonicalForm & F, const CanonicalForm & G )
697{
698    zz_pContext ccc(getCharacteristic());
699    ccc.restore();
700    zz_p::init(getCharacteristic());
701    zz_pX F1=convertFacCF2NTLzzpX(F);
702    zz_pX G1=convertFacCF2NTLzzpX(G);
703    zz_pX R=GCD(F1,G1);
704    return  convertNTLzzpX2CF(R,F.mvar());
705}
706
707#endif
708
709static bool
710gcd_avoid_mtaildegree ( CanonicalForm & f1, CanonicalForm & g1, CanonicalForm & d1 )
711{
712    bool rdy = true;
713    int df = f1.taildegree();
714    int dg = g1.taildegree();
715
716    d1 = d1.genOne();
717    if ( dg == 0 )
718    {
719        if ( df == 0 )
720            return false;
721        else
722        {
723            if ( f1.degree() == df )
724                d1 = cf_content( g1, LC( f1 ) );
725            else
726            {
727                f1 /= power( f1.mvar(), df );
728                rdy = false;
729            }
730        }
731    }
732    else
733    {
734        if ( df == 0)
735        {
736            if ( g1.degree() == dg )
737                d1 = cf_content( f1, LC( g1 ) );
738            else
739            {
740                g1 /= power( g1.mvar(), dg );
741                rdy = false;
742            }
743        }
744        else
745        {
746            if ( df > dg )
747                d1 = power( f1.mvar(), dg );
748            else
749                d1 = power( f1.mvar(), df );
750            if ( f1.degree() == df )
751            {
752                if (g1.degree() == dg)
753                    d1 *= gcd( LC( f1 ), LC( g1 ) );
754                else
755                {
756                    g1 /= power( g1.mvar(), dg);
757                    d1 *= cf_content( g1, LC( f1 ) );
758                }
759            }
760            else
761            {
762                f1 /= power( f1.mvar(), df );
763                if ( g1.degree() == dg )
764                    d1 *= cf_content( f1, LC( g1 ) );
765                else
766                {
767                    g1 /= power( g1.mvar(), dg );
768                    rdy = false;
769                }
770            }
771        }
772    }
773    return rdy;
774}
775
776/*
777*  compute positions p1 and pe of optimal variables:
778*    pe is used in "ezgcd" and
779*    p1 in "gcd_poly1"
780*/
781static
782void optvalues ( const int * df, const int * dg, const int n, int & p1, int &pe )
783{
784    int i, o1, oe;
785    if ( df[n] > dg[n] )
786    {
787        o1 = df[n]; oe = dg[n];
788    }
789    else
790    {
791        o1 = dg[n]; oe = df[n];
792    }
793    i = n-1;
794    while ( i > 0 )
795    {
796        if ( df[i] != 0 )
797        {
798            if ( df[i] > dg[i] )
799            {
800                if ( o1 > df[i]) { o1 = df[i]; p1 = i; }
801                if ( oe <= dg[i]) { oe = dg[i]; pe = i; }
802            }
803            else
804            {
805                if ( o1 > dg[i]) { o1 = dg[i]; p1 = i; }
806                if ( oe <= df[i]) { oe = df[i]; pe = i; }
807            }
808        }
809        i--;
810    }
811}
812
813/*
814*  make some changes of variables, see optvalues
815*/
816static void
817cf_prepgcd( const CanonicalForm & f, const CanonicalForm & g, int & cc, int & p1, int &pe )
818{
819    int i, k, n;
820    n = f.level();
821    cc = 0;
822    p1 = pe = n;
823    if ( n == 1 )
824        return;
825    int * degsf = new int[n+1];
826    int * degsg = new int[n+1];
827    for ( i = n; i > 0; i-- )
828    {
829        degsf[i] = degsg[i] = 0;
830    }
831    degsf = degrees( f, degsf );
832    degsg = degrees( g, degsg );
833
834    k = 0;
835    for ( i = n-1; i > 0; i-- )
836    {
837        if ( degsf[i] == 0 )
838        {
839            if ( degsg[i] != 0 )
840            {
841                cc = -i;
842                break;
843            }
844        }
845        else
846        {
847            if ( degsg[i] == 0 )
848            {
849                cc = i;
850                break;
851            }
852            else k++;
853        }
854    }
855
856    if ( ( cc == 0 ) && ( k != 0 ) )
857        optvalues( degsf, degsg, n, p1, pe );
858    if ( ( pe != 1 ) && ( degsf[1] != 0 ) )
859        pe = -pe;
860   
861    delete [] degsf;
862    delete [] degsg;
863}
Note: See TracBrowser for help on using the repository browser.