source: git/factory/cf_gcd.cc @ efdcc5

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