source: git/factory/cf_gcd.cc @ 98938c

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