source: git/factory/cf_gcd.cc @ 094eed

spielwiese
Last change on this file since 094eed was 094eed, checked in by Hans Schönemann <hannes@…>, 18 years ago
*hannes: ASSERT for cf_content git-svn-id: file:///usr/local/Singular/svn/trunk@8841 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 17.0 KB
Line 
1/* emacs edit mode for this file is -*- C++ -*- */
2/* $Id: cf_gcd.cc,v 1.35 2005-12-09 15:00:59 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 "cf_map.h"
16#include "sm_sparsemod.h"
17#include "fac_util.h"
18#include "ftmpl_functions.h"
19
20#ifdef HAVE_NTL
21#include <NTL/ZZX.h>
22#include "NTLconvert.h"
23bool isPurePoly(const CanonicalForm & f);
24#endif
25
26static CanonicalForm gcd_poly( const CanonicalForm & f, const CanonicalForm& g, bool modularflag );
27
28bool
29gcd_test_one ( const CanonicalForm & f, const CanonicalForm & g, bool swap )
30{
31    int count = 0;
32    // assume polys have same level;
33    CFRandom * sample = CFRandomFactory::generate();
34    REvaluation e( 2, tmax( f.level(), g.level() ), *sample );
35    delete sample;
36    CanonicalForm lcf, lcg;
37    if ( swap ) {
38        lcf = swapvar( LC( f ), Variable(1), f.mvar() );
39        lcg = swapvar( LC( g ), Variable(1), f.mvar() );
40    }
41    else {
42        lcf = LC( f, Variable(1) );
43        lcg = LC( g, Variable(1) );
44    }
45    while ( ( e( lcf ).isZero() || e( lcg ).isZero() ) && count < 100 ) {
46        e.nextpoint();
47        count++;
48    }
49    if ( count == 100 )
50        return false;
51    CanonicalForm F, G;
52    if ( swap ) {
53        F=swapvar( f, Variable(1), f.mvar() );
54        G=swapvar( g, Variable(1), g.mvar() );
55    }
56    else {
57        F = f;
58        G = g;
59    }
60    if ( e(F).taildegree() > 0 && e(G).taildegree() > 0 )
61        return false;
62    return gcd( e( F ), e( G ) ).degree() < 1;
63}
64
65//{{{ static CanonicalForm balance ( const CanonicalForm & f, const CanonicalForm & q )
66//{{{ docu
67//
68// balance() - map f from positive to symmetric representation
69//   mod q.
70//
71// This makes sense for univariate polynomials over Z only.
72// q should be an integer.
73//
74// Used by gcd_poly_univar0().
75//
76//}}}
77static CanonicalForm
78balance ( const CanonicalForm & f, const CanonicalForm & q )
79{
80    Variable x = f.mvar();
81    CanonicalForm result = 0, qh = q / 2;
82    CanonicalForm c;
83    CFIterator i;
84    for ( i = f; i.hasTerms(); i++ ) {
85        c = mod( i.coeff(), q );
86        if ( c > qh )
87            result += power( x, i.exp() ) * (c - q);
88        else
89            result += power( x, i.exp() ) * c;
90    }
91    return result;
92}
93//}}}
94
95//{{{ static CanonicalForm icontent ( const CanonicalForm & f, const CanonicalForm & c )
96//{{{ docu
97//
98// icontent() - return gcd of c and all coefficients of f which
99//   are in a coefficient domain.
100//
101// Used by icontent().
102//
103//}}}
104static CanonicalForm
105icontent ( const CanonicalForm & f, const CanonicalForm & c )
106{
107    if ( f.inCoeffDomain() )
108        return gcd( f, c );
109    else {
110        CanonicalForm g = c;
111        for ( CFIterator i = f; i.hasTerms() && ! g.isOne(); i++ )
112            g = icontent( i.coeff(), g );
113        return g;
114    }
115}
116//}}}
117
118//{{{ static CanonicalForm bcontent ( const CanonicalForm & f, const CanonicalForm & b )
119//{{{ docu
120//
121// bcontent() - return gcd of b and all coefficients of f which
122//   are in a basic domain.
123//
124// Used by gcd().
125//
126//}}}
127static CanonicalForm
128bcontent ( const CanonicalForm & f, const CanonicalForm & c )
129{
130    if ( f.inBaseDomain() )
131        return bgcd( f, c );
132    else if ( f.inCoeffDomain() )
133        return f.genOne();
134    else {
135        CanonicalForm g = c;
136        for ( CFIterator i = f; i.hasTerms() && ! g.isOne(); i++ )
137            g = bcontent( i.coeff(), g );
138        if( g.lc().sign() < 0 )
139            return -g;
140        else
141            return g;
142    }
143}
144//}}}
145
146
147//{{{ CanonicalForm icontent ( const CanonicalForm & f )
148//{{{ docu
149//
150// icontent() - return gcd over all coefficients of f which are
151//   in a coefficient domain.
152//
153//}}}
154CanonicalForm
155icontent ( const CanonicalForm & f )
156{
157    return icontent( f, 0 );
158}
159//}}}
160
161//{{{ CanonicalForm extgcd ( const CanonicalForm & f, const CanonicalForm & g, CanonicalForm & a, CanonicalForm & b )
162//{{{ docu
163//
164// extgcd() - returns polynomial extended gcd of f and g.
165//
166// Returns gcd(f, g) and a and b sucht that f*a+g*b=gcd(f, g).
167// The gcd is calculated using an extended euclidean polynomial
168// remainder sequence, so f and g should be polynomials over an
169// euclidean domain.  Normalizes result.
170//
171// Note: be sure that f and g have the same level!
172//
173//}}}
174CanonicalForm
175extgcd ( const CanonicalForm & f, const CanonicalForm & g, CanonicalForm & a, CanonicalForm & b )
176{
177#ifdef HAVE_NTL
178  if (isOn(SW_USE_NTL_GCD_P) && ( getCharacteristic() > 0 )
179  && isPurePoly(f) && isPurePoly(g))
180  {
181    zz_pContext ccc(getCharacteristic());
182    ccc.restore();
183    zz_p::init(getCharacteristic());
184    zz_pX F1=convertFacCF2NTLzzpX(f);
185    zz_pX G1=convertFacCF2NTLzzpX(g);
186    zz_pX R;
187    zz_pX A,B;
188    XGCD(R,A,B,F1,G1);
189    a=convertNTLzzpX2CF(A,f.mvar());
190    b=convertNTLzzpX2CF(B,f.mvar());
191    return convertNTLzzpX2CF(R,f.mvar());
192  }
193#endif
194  CanonicalForm contf = content( f );
195  CanonicalForm contg = content( g );
196
197  CanonicalForm p0 = f / contf, p1 = g / contg;
198  CanonicalForm f0 = 1, f1 = 0, g0 = 0, g1 = 1, q, r;
199
200  while ( ! p1.isZero() ) {
201      divrem( p0, p1, q, r );
202      p0 = p1; p1 = r;
203      r = g0 - g1 * q;
204      g0 = g1; g1 = r;
205      r = f0 - f1 * q;
206      f0 = f1; f1 = r;
207  }
208  CanonicalForm contp0 = content( p0 );
209  a = f0 / ( contf * contp0 );
210  b = g0 / ( contg * contp0 );
211  p0 /= contp0;
212  if ( p0.sign() < 0 ) {
213      p0 = -p0;
214      a = -a;
215      b = -b;
216  }
217  return p0;
218}
219//}}}
220
221static CanonicalForm
222gcd_poly_univar0( const CanonicalForm & F, const CanonicalForm & G, bool primitive )
223{
224#ifdef HAVE_NTL
225  if (isOn(SW_USE_NTL_GCD_P) && isPurePoly(F) && isPurePoly(G))
226  {
227    if ( getCharacteristic() > 0 )
228    {
229      //CanonicalForm cf=F.lc();
230      //CanonicalForm f=F / cf;
231      //CanonicalForm cg=G.lc();
232      //CanonicalForm g= G / cg;
233      zz_pContext ccc(getCharacteristic());
234      ccc.restore();
235      zz_p::init(getCharacteristic());
236      zz_pX F1=convertFacCF2NTLzzpX(F);
237      zz_pX G1=convertFacCF2NTLzzpX(G);
238      zz_pX R=GCD(F1,G1);
239      return  convertNTLzzpX2CF(R,F.mvar());
240    }
241    else
242    {
243      CanonicalForm f=F ;
244      CanonicalForm g=G ;
245      bool rat=isOn( SW_RATIONAL );
246      if ( isOn( SW_RATIONAL ) )
247      {
248         DEBOUTLN( cerr, "NTL_gcd: ..." );
249         CanonicalForm cdF = bCommonDen( F );
250         CanonicalForm cdG = bCommonDen( G );
251         Off( SW_RATIONAL );
252         CanonicalForm l = lcm( cdF, cdG );
253         On( SW_RATIONAL );
254         f *= l, g *= l;
255      }
256      DEBOUTLN( cerr, "NTL_gcd: f=" << f );
257      DEBOUTLN( cerr, "NTL_gcd: g=" << g );
258      ZZX F1=convertFacCF2NTLZZX(f);
259      ZZX G1=convertFacCF2NTLZZX(g);
260      ZZX R=GCD(F1,G1);
261      CanonicalForm r=convertNTLZZX2CF(R,F.mvar());
262      DEBOUTLN( cerr, "NTL_gcd: -> " << r );
263      if (rat)
264      {
265        r /= r.lc();
266        DEBOUTLN( cerr, "NTL_gcd2: -> " << r );
267      }
268      return r;
269    }
270  }
271#endif
272  CanonicalForm f, g, c, cg, cl, BB, B, M, q, Dp, newD, D, newq;
273  int p, i, n;
274
275  if ( primitive )
276  {
277    f = F;
278    g = G;
279    c = 1;
280  }
281  else
282  {
283    CanonicalForm cF = content( F ), cG = content( G );
284    f = F / cF;
285    g = G / cG;
286    c = bgcd( cF, cG );
287  }
288  cg = gcd( f.lc(), g.lc() );
289  cl = ( f.lc() / cg ) * g.lc();
290//     B = 2 * cg * tmin(
291//         maxnorm(f)*power(CanonicalForm(2),f.degree())*isqrt(f.degree()+1),
292//         maxnorm(g)*power(CanonicalForm(2),g.degree())*isqrt(g.degree()+1)
293//         )+1;
294  M = tmin( maxNorm(f), maxNorm(g) );
295  BB = power(CanonicalForm(2),tmin(f.degree(),g.degree()))*M;
296  q = 0;
297  i = cf_getNumSmallPrimes() - 1;
298  while ( true )
299  {
300    B = BB;
301    while ( i >= 0 && q < B )
302    {
303      p = cf_getSmallPrime( i );
304      i--;
305      while ( i >= 0 && mod( cl, p ) == 0 )
306      {
307        p = cf_getSmallPrime( i );
308        i--;
309      }
310      setCharacteristic( p );
311      Dp = gcd( mapinto( f ), mapinto( g ) );
312      Dp = ( Dp / Dp.lc() ) * mapinto( cg );
313      setCharacteristic( 0 );
314      if ( Dp.degree() == 0 )
315        return c;
316      if ( q.isZero() )
317      {
318        D = mapinto( Dp );
319        q = p;
320        B = power(CanonicalForm(2),D.degree())*M+1;
321      }
322      else
323      {
324        if ( Dp.degree() == D.degree() )
325        {
326          chineseRemainder( D, q, mapinto( Dp ), p, newD, newq );
327          q = newq;
328          D = newD;
329        }
330        else if ( Dp.degree() < D.degree() )
331        {
332          // all previous p's are bad primes
333          q = p;
334          D = mapinto( Dp );
335          B = power(CanonicalForm(2),D.degree())*M+1;
336        }
337        // else p is a bad prime
338      }
339    }
340    if ( i >= 0 )
341    {
342      // now balance D mod q
343      D = pp( balance( D, q ) );
344      if ( divides( D, f ) && divides( D, g ) )
345        return D * c;
346      else
347        q = 0;
348    }
349    else
350      return gcd_poly( F, G, false );
351    DEBOUTLN( cerr, "another try ..." );
352  }
353}
354
355CanonicalForm
356gcd_poly1( const CanonicalForm & f, const CanonicalForm & g, bool modularflag )
357{
358    CanonicalForm C, Ci, Ci1, Hi, bi, pi, pi1, pi2;
359    int delta;
360    Variable v = f.mvar();
361
362    if ( f.degree( v ) >= g.degree( v ) ) {
363        pi = f; pi1 = g;
364    }
365    else {
366        pi = g; pi1 = f;
367    }
368    Ci = content( pi ); Ci1 = content( pi1 );
369    C = gcd( Ci, Ci1 );
370    pi1 = pi1 / Ci1; pi = pi / Ci;
371    if ( pi.isUnivariate() && pi1.isUnivariate() )
372    {
373#ifdef HAVE_NTL
374      if ((isOn(SW_USE_NTL_GCD_P)||isOn(SW_USE_NTL_GCD_0))
375       && isPurePoly(pi) && isPurePoly(pi1))
376         return gcd_poly_univar0(f, g, true);
377#endif
378      if ( modularflag)
379        return gcd_poly_univar0( pi, pi1, true ) * C;
380    }
381    else if ( gcd_test_one( pi1, pi, true ) )
382      return C;
383    delta = degree( pi, v ) - degree( pi1, v );
384    Hi = power( LC( pi1, v ), delta );
385    if ( (delta+1) % 2 )
386        bi = 1;
387    else
388        bi = -1;
389    while ( degree( pi1, v ) > 0 ) {
390        pi2 = psr( pi, pi1, v );
391        pi2 = pi2 / bi;
392        pi = pi1; pi1 = pi2;
393        if ( degree( pi1, v ) > 0 ) {
394            delta = degree( pi, v ) - degree( pi1, v );
395            if ( (delta+1) % 2 )
396                bi = LC( pi, v ) * power( Hi, delta );
397            else
398                bi = -LC( pi, v ) * power( Hi, delta );
399            Hi = power( LC( pi1, v ), delta ) / power( Hi, delta-1 );
400        }
401    }
402    if ( degree( pi1, v ) == 0 )
403        return C;
404    else {
405        return C * pp( pi );
406    }
407}
408
409//{{{ static CanonicalForm gcd_poly ( const CanonicalForm & f, const CanonicalForm & g, bool modularflag )
410//{{{ docu
411//
412// gcd_poly() - calculate polynomial gcd.
413//
414// This is the dispatcher for polynomial gcd calculation.  We call either
415// ezgcd(), sparsemod() or gcd_poly1() in dependecy on the current
416// characteristic and settings of SW_USE_EZGCD and SW_USE_SPARSEMOD, resp.
417//
418// modularflag is reached down to gcd_poly1() without change in case of
419// zero characteristic.
420//
421// Used by gcd() and gcd_poly_univar0().
422//
423//}}}
424int si_factor_reminder=1;
425static CanonicalForm
426gcd_poly ( const CanonicalForm & f, const CanonicalForm & g, bool modularflag )
427{
428    if ( getCharacteristic() != 0 ) {
429        return gcd_poly1( f, g, false );
430    }
431    else if ( isOn( SW_USE_EZGCD ) && ! ( f.isUnivariate() && g.isUnivariate() ) ) {
432        CFMap M, N;
433        compress( f, g, M, N );
434#if 0
435        CanonicalForm r=N( ezgcd( M(f), M(g) ) );
436        if ((f%r!=0) || (g % r !=0))
437        {
438           if (si_factor_reminder)
439           printf("ezgcd failed, trying gcd_poly1\n");
440           return gcd_poly1( f, g, modularflag);
441        }
442        else
443           return r;
444#else
445         return N( ezgcd( M(f), M(g) ) );
446#endif
447    }
448    else if ( isOn( SW_USE_SPARSEMOD )
449    && ! ( f.isUnivariate() && g.isUnivariate() ) )
450    {
451#if 0
452        CanonicalForm r=sparsemod( f, g );
453        if ((f%r!=0) || (g % r !=0))
454        {
455           if (si_factor_reminder)
456           printf("sparsemod failed, trying gcd_poly1\n");
457           return r;
458           //return gcd_poly1( f, g, modularflag);
459        }
460        else
461          return r;
462#else
463        return sparsemod( f, g );
464#endif
465    }
466    else
467    {
468        return gcd_poly1( f, g, modularflag );
469    }
470}
471//}}}
472
473//{{{ static CanonicalForm cf_content ( const CanonicalForm & f, const CanonicalForm & g )
474//{{{ docu
475//
476// cf_content() - return gcd(g, content(f)).
477//
478// content(f) is calculated with respect to f's main variable.
479//
480// Used by gcd(), content(), content( CF, Variable ).
481//
482//}}}
483static CanonicalForm
484cf_content ( const CanonicalForm & f, const CanonicalForm & g )
485{
486    if ( f.inPolyDomain() || ( f.inExtension() && ! getReduce( f.mvar() ) ) ) {
487        CFIterator i = f;
488        CanonicalForm result = g;
489        while ( i.hasTerms() && ! result.isOne() ) {
490            result = gcd( result, i.coeff() );
491            i++;
492        }
493        return result;
494    }
495    else
496    {
497        ASSERT(g==0,"invalid call of cf_gcd");
498        if ( f.sign() < 0 )
499            return -f;
500        else
501            return f;
502    }
503}
504//}}}
505
506//{{{ CanonicalForm content ( const CanonicalForm & f )
507//{{{ docu
508//
509// content() - return content(f) with respect to main variable.
510//
511// Normalizes result.
512//
513//}}}
514CanonicalForm
515content ( const CanonicalForm & f )
516{
517    return cf_content( f, 0 );
518}
519//}}}
520
521//{{{ CanonicalForm content ( const CanonicalForm & f, const Variable & x )
522//{{{ docu
523//
524// content() - return content(f) with respect to x.
525//
526// x should be a polynomial variable.
527//
528//}}}
529CanonicalForm
530content ( const CanonicalForm & f, const Variable & x )
531{
532    ASSERT( x.level() > 0, "cannot calculate content with respect to algebraic variable" );
533    Variable y = f.mvar();
534
535    if ( y == x )
536        return cf_content( f, 0 );
537    else  if ( y < x )
538        return f;
539    else
540        return swapvar( content( swapvar( f, y, x ), y ), y, x );
541}
542//}}}
543
544//{{{ CanonicalForm vcontent ( const CanonicalForm & f, const Variable & x )
545//{{{ docu
546//
547// vcontent() - return content of f with repect to variables >= x.
548//
549// The content is recursively calculated over all coefficients in
550// f having level less than x.  x should be a polynomial
551// variable.
552//
553//}}}
554CanonicalForm
555vcontent ( const CanonicalForm & f, const Variable & x )
556{
557    ASSERT( x.level() > 0, "cannot calculate vcontent with respect to algebraic variable" );
558
559    if ( f.mvar() <= x )
560        return content( f, x );
561    else {
562        CFIterator i;
563        CanonicalForm d = 0;
564        for ( i = f; i.hasTerms() && ! d.isOne(); i++ )
565            d = gcd( d, vcontent( i.coeff(), x ) );
566        return d;
567    }
568}
569//}}}
570
571//{{{ CanonicalForm pp ( const CanonicalForm & f )
572//{{{ docu
573//
574// pp() - return primitive part of f.
575//
576// Returns zero if f equals zero, otherwise f / content(f).
577//
578//}}}
579CanonicalForm
580pp ( const CanonicalForm & f )
581{
582    if ( f.isZero() )
583        return f;
584    else
585        return f / content( f );
586}
587//}}}
588
589CanonicalForm
590gcd ( const CanonicalForm & f, const CanonicalForm & g )
591{
592    if ( f.isZero() )
593    {
594      if ( g.lc().sign() < 0 )
595        return -g;
596      else
597       return g;
598    }
599    else  if ( g.isZero() )
600    {
601      if ( f.lc().sign() < 0 )
602        return -f;
603      else
604        return f;
605    }
606    else  if ( f.inBaseDomain() )
607        return bcontent( g, f );
608    else  if ( g.inBaseDomain() )
609        return bcontent( f, g );
610    else  if ( f.mvar() == g.mvar() )
611        if ( f.inExtension() && getReduce( f.mvar() ) )
612            return 1;
613        else {
614            if ( divides( f, g ) )
615                if ( f.lc().sign() < 0 )
616                    return -f;
617                else
618                    return f;
619            else  if ( divides( g, f ) )
620                if ( g.lc().sign() < 0 )
621                    return -g;
622                else
623                    return g;
624            if ( getCharacteristic() == 0 && isOn( SW_RATIONAL ) ) {
625                CanonicalForm cdF = bCommonDen( f );
626                CanonicalForm cdG = bCommonDen( g );
627                Off( SW_RATIONAL );
628                CanonicalForm l = lcm( cdF, cdG );
629                On( SW_RATIONAL );
630                CanonicalForm F = f * l, G = g * l;
631                Off( SW_RATIONAL );
632                do { l = gcd_poly( F, G, true ); }
633                while ((!divides(l,F)) || (!divides(l,G)));
634                On( SW_RATIONAL );
635                if ( l.lc().sign() < 0 )
636                    return -l;
637                else
638                    return l;
639            }
640            else {
641                CanonicalForm d;
642                do{ d = gcd_poly( f, g, getCharacteristic()==0 ); }
643                while ((!divides(d,f)) || (!divides(d,g)));
644                if ( d.lc().sign() < 0 )
645                    return -d;
646                else
647                    return d;
648            }
649        }
650    else  if ( f.mvar() > g.mvar() )
651        return cf_content( f, g );
652    else
653        return cf_content( g, f );
654}
655
656//{{{ CanonicalForm lcm ( const CanonicalForm & f, const CanonicalForm & g )
657//{{{ docu
658//
659// lcm() - return least common multiple of f and g.
660//
661// The lcm is calculated using the formula lcm(f, g) = f * g / gcd(f, g).
662//
663// Returns zero if one of f or g equals zero.
664//
665//}}}
666CanonicalForm
667lcm ( const CanonicalForm & f, const CanonicalForm & g )
668{
669    if ( f.isZero() || g.isZero() )
670        return f;
671    else
672        return ( f / gcd( f, g ) ) * g;
673}
674//}}}
Note: See TracBrowser for help on using the repository browser.