source: git/factory/cf_gcd.cc @ 7befae

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