source: git/factory/cf_gcd.cc @ f11d7b

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