source: git/factory/cf_gcd.cc @ 1986a2

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