source: git/factory/cf_gcd.cc @ 1baf4f

spielwiese
Last change on this file since 1baf4f was 1baf4f, checked in by Hans Schönemann <hannes@…>, 20 years ago
*hannes: retry for gcd git-svn-id: file:///usr/local/Singular/svn/trunk@7489 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 14.9 KB
Line 
1/* emacs edit mode for this file is -*- C++ -*- */
2/* $Id: cf_gcd.cc,v 1.24 2004-09-23 16:51:48 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//}}}
368static CanonicalForm
369gcd_poly ( const CanonicalForm & f, const CanonicalForm & g, bool modularflag )
370{
371    if ( getCharacteristic() != 0 ) {
372        return gcd_poly1( f, g, false );
373    }
374    else if ( isOn( SW_USE_EZGCD ) && ! ( f.isUnivariate() && g.isUnivariate() ) ) {
375        CFMap M, N;
376        compress( f, g, M, N );
377        return N( ezgcd( M(f), M(g) ) );
378    }
379    else if ( isOn( SW_USE_SPARSEMOD ) && ! ( f.isUnivariate() && g.isUnivariate() ) ) {
380        return sparsemod( f, g );
381    }
382    else {
383        return gcd_poly1( f, g, modularflag );
384    }
385}
386//}}}
387
388//{{{ static CanonicalForm cf_content ( const CanonicalForm & f, const CanonicalForm & g )
389//{{{ docu
390//
391// cf_content() - return gcd(g, content(f)).
392//
393// content(f) is calculated with respect to f's main variable.
394//
395// Used by gcd(), content(), content( CF, Variable ).
396//
397//}}}
398static CanonicalForm
399cf_content ( const CanonicalForm & f, const CanonicalForm & g )
400{
401    if ( f.inPolyDomain() || ( f.inExtension() && ! getReduce( f.mvar() ) ) ) {
402        CFIterator i = f;
403        CanonicalForm result = g;
404        while ( i.hasTerms() && ! result.isOne() ) {
405            result = gcd( result, i.coeff() );
406            i++;
407        }
408        return result;
409    }
410    else
411        if ( f.sign() < 0 )
412            return -f;
413        else
414            return f;
415}
416//}}}
417
418//{{{ CanonicalForm content ( const CanonicalForm & f )
419//{{{ docu
420//
421// content() - return content(f) with respect to main variable.
422//
423// Normalizes result.
424//
425//}}}
426CanonicalForm
427content ( const CanonicalForm & f )
428{
429    return cf_content( f, 0 );
430}
431//}}}
432
433//{{{ CanonicalForm content ( const CanonicalForm & f, const Variable & x )
434//{{{ docu
435//
436// content() - return content(f) with respect to x.
437//
438// x should be a polynomial variable.
439//
440//}}}
441CanonicalForm
442content ( const CanonicalForm & f, const Variable & x )
443{
444    ASSERT( x.level() > 0, "cannot calculate content with respect to algebraic variable" );
445    Variable y = f.mvar();
446
447    if ( y == x )
448        return cf_content( f, 0 );
449    else  if ( y < x )
450        return f;
451    else
452        return swapvar( content( swapvar( f, y, x ), y ), y, x );
453}
454//}}}
455
456//{{{ CanonicalForm vcontent ( const CanonicalForm & f, const Variable & x )
457//{{{ docu
458//
459// vcontent() - return content of f with repect to variables >= x.
460//
461// The content is recursively calculated over all coefficients in
462// f having level less than x.  x should be a polynomial
463// variable.
464//
465//}}}
466CanonicalForm
467vcontent ( const CanonicalForm & f, const Variable & x )
468{
469    ASSERT( x.level() > 0, "cannot calculate vcontent with respect to algebraic variable" );
470
471    if ( f.mvar() <= x )
472        return content( f, x );
473    else {
474        CFIterator i;
475        CanonicalForm d = 0;
476        for ( i = f; i.hasTerms() && ! d.isOne(); i++ )
477            d = gcd( d, vcontent( i.coeff(), x ) );
478        return d;
479    }
480}
481//}}}
482
483//{{{ CanonicalForm pp ( const CanonicalForm & f )
484//{{{ docu
485//
486// pp() - return primitive part of f.
487//
488// Returns zero if f equals zero, otherwise f / content(f).
489//
490//}}}
491CanonicalForm
492pp ( const CanonicalForm & f )
493{
494    if ( f.isZero() )
495        return f;
496    else
497        return f / content( f );
498}
499//}}}
500
501CanonicalForm
502gcd ( const CanonicalForm & f, const CanonicalForm & g )
503{
504    if ( f.isZero() )
505        if ( g.lc().sign() < 0 )
506            return -g;
507        else
508            return g;
509    else  if ( g.isZero() )
510        if ( f.lc().sign() < 0 )
511            return -f;
512        else
513            return f;
514    else  if ( f.inBaseDomain() )
515        if ( g.inBaseDomain() )
516            return bgcd( f, g );
517        else
518            return cf_content( g, f );
519    else  if ( g.inBaseDomain() )
520        return cf_content( f, g );
521    else  if ( f.mvar() == g.mvar() )
522        if ( f.inExtension() && getReduce( f.mvar() ) )
523            return 1;
524        else {
525            if ( divides( f, g ) )
526                if ( f.lc().sign() < 0 )
527                    return -f;
528                else
529                    return f;
530            else  if ( divides( g, f ) )
531                if ( g.lc().sign() < 0 )
532                    return -g;
533                else
534                    return g;
535            if ( getCharacteristic() == 0 && isOn( SW_RATIONAL ) ) {
536                CanonicalForm cdF = bCommonDen( f );
537                CanonicalForm cdG = bCommonDen( g );
538                Off( SW_RATIONAL );
539                CanonicalForm l = lcm( cdF, cdG );
540                On( SW_RATIONAL );
541                CanonicalForm F = f * l, G = g * l;
542                Off( SW_RATIONAL );
543                l = gcd_poly( F, G, true );
544                if ((F%l!=0) || (G % l !=0))
545                  l = gcd_poly( F, G, true );
546                On( SW_RATIONAL );
547                if ( l.lc().sign() < 0 )
548                    return -l;
549                else
550                    return l;
551            }
552            else {
553                CanonicalForm d = gcd_poly( f, g, getCharacteristic()==0 );
554                if ((f%d!=0) || (g % d !=0))
555                  d = gcd_poly( f, g, getCharacteristic()==0  );
556                if ( d.lc().sign() < 0 )
557                    return -d;
558                else
559                    return d;
560            }
561        }
562    else  if ( f.mvar() > g.mvar() )
563        return cf_content( f, g );
564    else
565        return cf_content( g, f );
566}
567
568//{{{ CanonicalForm lcm ( const CanonicalForm & f, const CanonicalForm & g )
569//{{{ docu
570//
571// lcm() - return least common multiple of f and g.
572//
573// The lcm is calculated using the formula lcm(f, g) = f * g / gcd(f, g).
574//
575// Returns zero if one of f or g equals zero.
576//
577//}}}
578CanonicalForm
579lcm ( const CanonicalForm & f, const CanonicalForm & g )
580{
581    if ( f.isZero() || g.isZero() )
582        return f;
583    else
584        return ( f / gcd( f, g ) ) * g;
585}
586//}}}
Note: See TracBrowser for help on using the repository browser.