source: git/factory/cf_gcd.cc @ 4bca3e

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