source: git/factory/cf_gcd.cc @ 398b15

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