source: git/factory/cf_gcd.cc @ edb4893

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