source: git/factory/cf_gcd.cc @ c97005

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