source: git/factory/cf_gcd.cc @ 84250a6

spielwiese
Last change on this file since 84250a6 was ab4548f, checked in by Jens Schmidt <schmidt@…>, 27 years ago
#include <config.h> added git-svn-id: file:///usr/local/Singular/svn/trunk@133 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 12.0 KB
Line 
1// emacs edit mode for this file is -*- C++ -*-
2// $Id: cf_gcd.cc,v 1.7 1997-04-07 15:05:38 schmidt Exp $
3
4/*
5$Log: not supported by cvs2svn $
6Revision 1.6  1997/03/26 16:39:06  schmidt
7debug output changed to DEBOUT
8
9Revision 1.5  1996/12/05 18:24:53  schmidt
10``Unconditional'' check-in.
11Now it is my turn to develop factory.
12
13Revision 1.4  1996/07/08 08:21:10  stobbe
14"gcd_poly: now uses ezgcd if the switch SW_USE_EZGCD is on.
15"
16
17Revision 1.3  1996/06/18 12:22:54  stobbe
18"gcd_poly_univar0: now uses getSmallPrimes (due to changes in the handling
19                  of prime numbers in cf_primes.
20"
21
22Revision 1.2  1996/06/13 08:18:34  stobbe
23"balance: Now balaces polynomials even if the coefficient sizes are higher
24         than the bound.
25gcd: Now returns the results with positive leading coefficient.
26     The isOne test is now performed if pi or pi1 is multivariate.
27"
28
29Revision 1.1  1996/06/03 08:32:56  stobbe
30"gcd_poly: now uses new function gcd_poly_univar0 to compute univariate
31          polynomial gcd's over Z.
32gcd_poly_univar0: computes univariate polynomial gcd's in characteristic 0
33                  via chinese remaindering.
34maxnorm: computes the maximum norm of all coefficients of a polynomial.
35"
36
37Revision 1.0  1996/05/17 11:56:37  stobbe
38Initial revision
39
40*/
41
42#include <config.h>
43
44#include "assert.h"
45#include "debug.h"
46
47#include "cf_defs.h"
48#include "canonicalform.h"
49#include "cf_iter.h"
50#include "cf_reval.h"
51#include "cf_primes.h"
52#include "cf_chinese.h"
53#include "cf_map.h"
54#include "fac_util.h"
55#include "templates/functions.h"
56
57static CanonicalForm gcd_poly( const CanonicalForm & f, const CanonicalForm& g, bool modularflag );
58
59
60static int
61isqrt ( int a )
62{
63    int h, x0, x1 = a;
64    do {
65        x0 = x1;
66        h = x0 * x0 + a - 1;
67        if ( h % (2 * x0) == 0 )
68            x1 = h / (2 * x0);
69        else
70            x1 = (h - 1)  / (2 * x0);
71    } while ( x1 < x0 );
72    return x1;
73}
74
75bool
76gcd_test_one ( const CanonicalForm & f, const CanonicalForm & g, bool swap )
77{
78    int count = 0;
79    // assume polys have same level;
80    CFRandom * sample = CFRandomFactory::generate();
81    REvaluation e( 2, tmax( f.level(), g.level() ), *sample );
82    delete sample;
83    CanonicalForm lcf, lcg;
84    if ( swap ) {
85        lcf = swapvar( LC( f ), Variable(1), f.mvar() );
86        lcg = swapvar( LC( g ), Variable(1), f.mvar() );
87    }
88    else {
89        lcf = LC( f, Variable(1) );
90        lcg = LC( g, Variable(1) );
91    }
92    while ( ( e( lcf ).isZero() || e( lcg ).isZero() ) && count < 100 ) {
93        e.nextpoint();
94        count++;
95    }
96    if ( count == 100 )
97        return false;
98    CanonicalForm F, G;
99    if ( swap ) {
100        F=swapvar( f, Variable(1), f.mvar() );
101        G=swapvar( g, Variable(1), g.mvar() );
102    }
103    else {
104        F = f;
105        G = g;
106    }
107    if ( e(F).taildegree() > 0 && e(G).taildegree() > 0 )
108        return false;
109    return gcd( e( F ), e( G ) ).degree() < 1;
110}
111
112static CanonicalForm
113maxnorm ( const CanonicalForm & f )
114{
115    CanonicalForm m = 0, h;
116    CFIterator i;
117    for ( i = f; i.hasTerms(); i++ )
118        m = tmax( m, abs( i.coeff() ) );
119    return m;
120}
121
122static void
123chinesePoly ( const CanonicalForm & f1, const CanonicalForm & q1, const CanonicalForm & f2, const CanonicalForm & q2, CanonicalForm & f, CanonicalForm & q )
124{
125    CFIterator i1 = f1, i2 = f2;
126    CanonicalForm c;
127    Variable x = f1.mvar();
128    f = 0;
129    while ( i1.hasTerms() && i2.hasTerms() ) {
130        if ( i1.exp() == i2.exp() ) {
131            chineseRemainder( i1.coeff(), q1, i2.coeff(), q2, c, q );
132            f += power( x, i1.exp() ) * c;
133            i1++; i2++;
134        }
135        else if ( i1.exp() > i2.exp() ) {
136            chineseRemainder( 0, q1, i2.coeff(), q2, c, q );
137            f += power( x, i2.exp() ) * c;
138            i2++;
139        }
140        else {
141            chineseRemainder( i1.coeff(), q1, 0, q2, c, q );
142            f += power( x, i1.exp() ) * c;
143            i1++;
144        }
145    }
146    while ( i1.hasTerms() ) {
147        chineseRemainder( i1.coeff(), q1, 0, q2, c, q );
148        f += power( x, i1.exp() ) * c;
149        i1++;
150    }
151    while ( i2.hasTerms() ) {
152        chineseRemainder( 0, q1, i2.coeff(), q2, c, q );
153        f += power( x, i2.exp() ) * c;
154        i2++;
155    }
156}
157
158static CanonicalForm
159balance ( const CanonicalForm & f, const CanonicalForm & q )
160{
161    CFIterator i;
162    CanonicalForm result = 0, qh = q / 2;
163    CanonicalForm c;
164    Variable x = f.mvar();
165    for ( i = f; i.hasTerms(); i++ ) {
166        c = i.coeff() % q;
167        if ( c > qh )
168            result += power( x, i.exp() ) * (c - q);
169        else
170            result += power( x, i.exp() ) * c;
171    }
172    return result;
173}
174
175CanonicalForm
176igcd ( const CanonicalForm & f, const CanonicalForm & g )
177{
178    CanonicalForm a, b, c, dummy;
179
180    if ( f.inZ() && g.inZ() && ! isOn( SW_RATIONAL ) ) {
181        if ( f.sign() < 0 ) a = -f; else a = f;
182        if ( g.sign() < 0 ) b = -g; else b = g;
183        while ( ! b.isZero() ) {
184            divrem( a, b, dummy, c );
185            a = b;
186            b = c;
187        }
188        return a;
189    }
190    else
191        return 1;
192}
193
194static CanonicalForm
195icontent ( const CanonicalForm & f, const CanonicalForm & c )
196{
197    if ( f.inCoeffDomain() )
198        return gcd( f, c );
199    else {
200        CanonicalForm g = c;
201        for ( CFIterator i = f; i.hasTerms() && ! g.isOne(); i++ )
202            g = icontent( i.coeff(), g );
203        return g;
204    }
205}
206
207CanonicalForm
208icontent ( const CanonicalForm & f )
209{
210    return icontent( f, 0 );
211}
212
213CanonicalForm
214iextgcd ( const CanonicalForm & f, const CanonicalForm & g, CanonicalForm & a, CanonicalForm & b )
215{
216    CanonicalForm p0 = f, p1 = g;
217    CanonicalForm f0 = 1, f1 = 0, g0 = 0, g1 = 1, q, r;
218
219    while ( ! p1.isZero() ) {
220        divrem( p0, p1, q, r );
221        p0 = p1; p1 = r;
222        r = g0 - g1 * q;
223        g0 = g1; g1 = r;
224        r = f0 - f1 * q;
225        f0 = f1; f1 = r;
226    }
227    a = f0;
228    b = g0;
229    if ( p0.sign() < 0 ) {
230        p0 = -p0;
231        a = -a;
232        b = -b;
233    }
234    return p0;
235}
236
237CanonicalForm
238extgcd ( const CanonicalForm & f, const CanonicalForm & g, CanonicalForm & a, CanonicalForm & b )
239{
240    CanonicalForm contf = content( f );
241    CanonicalForm contg = content( g );
242
243    CanonicalForm p0 = f / contf, p1 = g / contg;
244    CanonicalForm f0 = 1, f1 = 0, g0 = 0, g1 = 1, q, r;
245
246    while ( ! p1.isZero() ) {
247        divrem( p0, p1, q, r );
248        p0 = p1; p1 = r;
249        r = g0 - g1 * q;
250        g0 = g1; g1 = r;
251        r = f0 - f1 * q;
252        f0 = f1; f1 = r;
253    }
254    CanonicalForm contp0 = content( p0 );
255    a = f0 / ( contf * contp0 );
256    b = g0 / ( contg * contp0 );
257    p0 /= contp0;
258    if ( p0.sign() < 0 ) {
259        p0 = -p0;
260        a = -a;
261        b = -b;
262    }
263    return p0;
264}
265
266static CanonicalForm
267gcd_poly_univar0( const CanonicalForm & F, const CanonicalForm & G, bool primitive )
268{
269    CanonicalForm f, g, c, cg, cl, BB, B, M, q, Dp, newD, D, newq;
270    int p, i, n;
271
272    if ( primitive ) {
273        f = F;
274        g = G;
275        c = 1;
276    }
277    else {
278        CanonicalForm cF = content( F ), cG = content( G );
279        f = F / cF;
280        g = G / cG;
281        c = igcd( cF, cG );
282    }
283    cg = gcd( f.lc(), g.lc() );
284    cl = ( f.lc() / cg ) * g.lc();
285//     B = 2 * cg * tmin(
286//      maxnorm(f)*power(CanonicalForm(2),f.degree())*isqrt(f.degree()+1),
287//      maxnorm(g)*power(CanonicalForm(2),g.degree())*isqrt(g.degree()+1)
288//      )+1;
289    M = tmin( maxnorm(f), maxnorm(g) );
290    BB = power(CanonicalForm(2),tmin(f.degree(),g.degree()))*M;
291    q = 0;
292    i = cf_getNumSmallPrimes() - 1;
293    while ( true ) {
294        B = BB;
295        while ( i >= 0 && q < B ) {
296            p = cf_getSmallPrime( i );
297            i--;
298            while ( i >= 0 && mod( cl, p ) == 0 ) {
299                p = cf_getSmallPrime( i );
300                i--;
301            }
302            setCharacteristic( p );
303            Dp = gcd( mapinto( f ), mapinto( g ) );
304            Dp = ( Dp / Dp.lc() ) * mapinto( cg );
305            setCharacteristic( 0 );
306            if ( Dp.degree() == 0 ) return c;
307            if ( q.isZero() ) {
308                D = mapinto( Dp );
309                q = p;
310                B = power(CanonicalForm(2),D.degree())*M+1;
311            }
312            else {
313                if ( Dp.degree() == D.degree() ) {
314                    chinesePoly( D, q, mapinto( Dp ), p, newD, newq );
315                    q = newq;
316                    D = newD;
317                }
318                else if ( Dp.degree() < D.degree() ) {
319                    // all previous p's are bad primes
320                    q = p;
321                    D = mapinto( Dp );
322                    B = power(CanonicalForm(2),D.degree())*M+1;
323                }
324                // else p is a bad prime
325            }
326        }
327        if ( i >= 0 ) {
328            // now balance D mod q
329            D = pp( balance( D, q ) );
330            if ( divides( D, f ) && divides( D, g ) )
331                return D * c;
332            else
333                q = 0;
334        }
335        else {
336            return gcd_poly( F, G, false );
337        }
338        DEBOUTLN( cerr, "another try ...", ' ' );
339    }
340}
341
342
343static CanonicalForm
344gcd_poly1( const CanonicalForm & f, const CanonicalForm & g, bool modularflag )
345{
346    CanonicalForm C, Ci, Ci1, Hi, bi, pi, pi1, pi2;
347    int delta;
348    Variable v = f.mvar();
349
350    if ( f.degree( v ) >= g.degree( v ) ) {
351        pi = f; pi1 = g;
352    }
353    else {
354        pi = g; pi1 = f;
355    }
356    Ci = content( pi ); Ci1 = content( pi1 );
357    C = gcd( Ci, Ci1 );
358    pi1 = pi1 / Ci1; pi = pi / Ci;
359    if ( pi.isUnivariate() && pi1.isUnivariate() ) {
360        if ( modularflag )
361            return gcd_poly_univar0( pi, pi1, true ) * C;
362    }
363    else
364        if ( gcd_test_one( pi1, pi, true ) )
365            return C;
366    delta = degree( pi, v ) - degree( pi1, v );
367    Hi = power( LC( pi1, v ), delta );
368    if ( (delta+1) % 2 )
369        bi = 1;
370    else
371        bi = -1;
372    while ( degree( pi1, v ) > 0 ) {
373        pi2 = psr( pi, pi1, v );
374        pi2 = pi2 / bi;
375        pi = pi1; pi1 = pi2;
376        if ( degree( pi1, v ) > 0 ) {
377            delta = degree( pi, v ) - degree( pi1, v );
378            if ( (delta+1) % 2 )
379                bi = LC( pi, v ) * power( Hi, delta );
380            else
381                bi = -LC( pi, v ) * power( Hi, delta );
382            Hi = power( LC( pi1, v ), delta ) / power( Hi, delta-1 );
383        }
384    }
385    if ( degree( pi1, v ) == 0 )
386        return C;
387    else {
388        return C * pp( pi );
389    }
390}
391
392
393static CanonicalForm
394gcd_poly( const CanonicalForm & f, const CanonicalForm & g, bool modularflag )
395{
396    if ( getCharacteristic() != 0 ) {
397        return gcd_poly1( f, g, false );
398    }
399    else if ( isOn( SW_USE_EZGCD ) && ! ( f.isUnivariate() && g.isUnivariate() ) ) {
400        CFMap M, N;
401        compress( f, g, M, N );
402        return N( ezgcd( M(f), M(g) ) );
403    }
404    else {
405        return gcd_poly1( f, g, modularflag );
406    }
407}
408
409
410static CanonicalForm
411cf_content ( const CanonicalForm & f, const CanonicalForm & g )
412{
413    if ( f.inPolyDomain() || ( f.inExtension() && ! getReduce( f.mvar() ) ) ) {
414        CFIterator i = f;
415        CanonicalForm result = g;
416        while ( i.hasTerms() && ! result.isOne() ) {
417            result = gcd( result, i.coeff() );
418            i++;
419        }
420        return result;
421    }
422    else
423        if ( f.sign() < 0 )
424            return -f;
425        else
426            return f;
427}
428
429CanonicalForm
430content ( const CanonicalForm & f )
431{
432    return cf_content( f, 0 );
433}
434
435CanonicalForm
436content ( const CanonicalForm & f, const Variable & x )
437{
438    if ( f.mvar() == x )
439        return cf_content( f, 0 );
440    else  if ( f.mvar() < x )
441        return f;
442    else
443        return swapvar( content( swapvar( f, f.mvar(), x ), f.mvar() ), f.mvar(), x );
444}
445
446CanonicalForm
447vcontent ( const CanonicalForm & f, const Variable & x )
448{
449    if ( f.mvar() <= x )
450        return content( f, x );
451    else {
452        CFIterator i;
453        CanonicalForm d = 0;
454        for ( i = f; i.hasTerms() && ! d.isOne(); i++ )
455            d = gcd( d, vcontent( i.coeff(), x ) );
456        return d;
457    }
458}
459
460CanonicalForm
461pp ( const CanonicalForm & f )
462{
463    if ( f.isZero() )
464        return f;
465    else
466        return f / content( f );
467}
468
469CanonicalForm
470gcd ( const CanonicalForm & f, const CanonicalForm & g )
471{
472    if ( f.isZero() )
473        if ( g.lc().sign() < 0 )
474            return -g;
475        else
476            return g;
477    else  if ( g.isZero() )
478        if ( f.lc().sign() < 0 )
479            return -f;
480        else
481            return f;
482    else  if ( f.inBaseDomain() )
483        if ( g.inBaseDomain() )
484            return igcd( f, g );
485        else
486            return cf_content( g, f );
487    else  if ( g.inBaseDomain() )
488        return cf_content( f, g );
489    else  if ( f.mvar() == g.mvar() )
490        if ( f.inExtension() && getReduce( f.mvar() ) )
491            return 1;
492        else {
493            if ( divides( f, g ) )
494                if ( f.lc().sign() < 0 )
495                    return -f;
496                else
497                    return f;
498            else  if ( divides( g, f ) )
499                if ( g.lc().sign() < 0 )
500                    return -g;
501                else
502                    return g;
503            if ( getCharacteristic() == 0 && isOn( SW_RATIONAL ) ) {
504                Off( SW_RATIONAL );
505                CanonicalForm l = lcm( common_den( f ), common_den( g ) );
506                On( SW_RATIONAL );
507                CanonicalForm F = f * l, G = g * l;
508                Off( SW_RATIONAL );
509                l = gcd_poly( F, G, true );
510                On( SW_RATIONAL );
511                if ( l.lc().sign() < 0 )
512                    return -l;
513                else
514                    return l;
515            }
516            else {
517                CanonicalForm d = gcd_poly( f, g, getCharacteristic()==0 );
518                if ( d.lc().sign() < 0 )
519                    return -d;
520                else
521                    return d;
522            }
523        }
524    else  if ( f.mvar() > g.mvar() )
525        return cf_content( f, g );
526    else
527        return cf_content( g, f );
528}
529
530CanonicalForm
531lcm ( const CanonicalForm & f, const CanonicalForm & g )
532{
533    return ( f / gcd( f, g ) ) * g;
534}
Note: See TracBrowser for help on using the repository browser.