source: git/factory/cf_gcd.cc @ 0c4a34b

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