source: git/factory/cf_gcd.cc @ cd86ac

spielwiese
Last change on this file since cd86ac was 150dc8, checked in by Hans Schönemann <hannes@…>, 23 years ago
*** empty log message *** git-svn-id: file:///usr/local/Singular/svn/trunk@5140 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 13.9 KB
Line 
1/* emacs edit mode for this file is -*- C++ -*- */
2/* $Id: cf_gcd.cc,v 1.21 2001-01-30 13:41:12 Singular 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_algorithm.h"
15#include "cf_map.h"
16#include "sm_sparsemod.h"
17#include "fac_util.h"
18#include "ftmpl_functions.h"
19
20static CanonicalForm gcd_poly( const CanonicalForm & f, const CanonicalForm& g, bool modularflag );
21
22bool
23gcd_test_one ( const CanonicalForm & f, const CanonicalForm & g, bool swap )
24{
25    int count = 0;
26    // assume polys have same level;
27    CFRandom * sample = CFRandomFactory::generate();
28    REvaluation e( 2, tmax( f.level(), g.level() ), *sample );
29    delete sample;
30    CanonicalForm lcf, lcg;
31    if ( swap ) {
32        lcf = swapvar( LC( f ), Variable(1), f.mvar() );
33        lcg = swapvar( LC( g ), Variable(1), f.mvar() );
34    }
35    else {
36        lcf = LC( f, Variable(1) );
37        lcg = LC( g, Variable(1) );
38    }
39    while ( ( e( lcf ).isZero() || e( lcg ).isZero() ) && count < 100 ) {
40        e.nextpoint();
41        count++;
42    }
43    if ( count == 100 )
44        return false;
45    CanonicalForm F, G;
46    if ( swap ) {
47        F=swapvar( f, Variable(1), f.mvar() );
48        G=swapvar( g, Variable(1), g.mvar() );
49    }
50    else {
51        F = f;
52        G = g;
53    }
54    if ( e(F).taildegree() > 0 && e(G).taildegree() > 0 )
55        return false;
56    return gcd( e( F ), e( G ) ).degree() < 1;
57}
58
59//{{{ static CanonicalForm balance ( const CanonicalForm & f, const CanonicalForm & q )
60//{{{ docu
61//
62// balance() - map f from positive to symmetric representation
63//   mod q.
64//
65// This makes sense for univariate polynomials over Z only.
66// q should be an integer.
67//
68// Used by gcd_poly_univar0().
69//
70//}}}
71static CanonicalForm
72balance ( const CanonicalForm & f, const CanonicalForm & q )
73{
74    Variable x = f.mvar();
75    CanonicalForm result = 0, qh = q / 2;
76    CanonicalForm c;
77    CFIterator i;
78    for ( i = f; i.hasTerms(); i++ ) {
79        c = mod( i.coeff(), q );
80        if ( c > qh )
81            result += power( x, i.exp() ) * (c - q);
82        else
83            result += power( x, i.exp() ) * c;
84    }
85    return result;
86}
87//}}}
88
89//{{{ static CanonicalForm icontent ( const CanonicalForm & f, const CanonicalForm & c )
90//{{{ docu
91//
92// icontent() - return gcd of c and all coefficients of f which
93//   are in a coefficient domain.
94//
95// Used by icontent().
96//
97//}}}
98static CanonicalForm
99icontent ( const CanonicalForm & f, const CanonicalForm & c )
100{
101    if ( f.inCoeffDomain() )
102        return gcd( f, c );
103    else {
104        CanonicalForm g = c;
105        for ( CFIterator i = f; i.hasTerms() && ! g.isOne(); i++ )
106            g = icontent( i.coeff(), g );
107        return g;
108    }
109}
110//}}}
111
112//{{{ CanonicalForm icontent ( const CanonicalForm & f )
113//{{{ docu
114//
115// icontent() - return gcd over all coefficients of f which are
116//   in a coefficient domain.
117//
118//}}}
119CanonicalForm
120icontent ( const CanonicalForm & f )
121{
122    return icontent( f, 0 );
123}
124//}}}
125
126//{{{ CanonicalForm extgcd ( const CanonicalForm & f, const CanonicalForm & g, CanonicalForm & a, CanonicalForm & b )
127//{{{ docu
128//
129// extgcd() - returns polynomial extended gcd of f and g.
130//
131// Returns gcd(f, g) and a and b sucht that f*a+g*b=gcd(f, g).
132// The gcd is calculated using an extended euclidean polynomial
133// remainder sequence, so f and g should be polynomials over an
134// euclidean domain.  Normalizes result.
135//
136// Note: be sure that f and g have the same level!
137//
138//}}}
139CanonicalForm
140extgcd ( const CanonicalForm & f, const CanonicalForm & g, CanonicalForm & a, CanonicalForm & b )
141{
142    CanonicalForm contf = content( f );
143    CanonicalForm contg = content( g );
144
145    CanonicalForm p0 = f / contf, p1 = g / contg;
146    CanonicalForm f0 = 1, f1 = 0, g0 = 0, g1 = 1, q, r;
147
148    while ( ! p1.isZero() ) {
149        divrem( p0, p1, q, r );
150        p0 = p1; p1 = r;
151        r = g0 - g1 * q;
152        g0 = g1; g1 = r;
153        r = f0 - f1 * q;
154        f0 = f1; f1 = r;
155    }
156    CanonicalForm contp0 = content( p0 );
157    a = f0 / ( contf * contp0 );
158    b = g0 / ( contg * contp0 );
159    p0 /= contp0;
160    if ( p0.sign() < 0 ) {
161        p0 = -p0;
162        a = -a;
163        b = -b;
164    }
165    return p0;
166}
167//}}}
168
169static CanonicalForm
170gcd_poly_univar0( const CanonicalForm & F, const CanonicalForm & G, bool primitive )
171{
172    CanonicalForm f, g, c, cg, cl, BB, B, M, q, Dp, newD, D, newq;
173    int p, i, n;
174
175    if ( primitive ) {
176        f = F;
177        g = G;
178        c = 1;
179    }
180    else {
181        CanonicalForm cF = content( F ), cG = content( G );
182        f = F / cF;
183        g = G / cG;
184        c = bgcd( cF, cG );
185    }
186    cg = gcd( f.lc(), g.lc() );
187    cl = ( f.lc() / cg ) * g.lc();
188//     B = 2 * cg * tmin(
189//         maxnorm(f)*power(CanonicalForm(2),f.degree())*isqrt(f.degree()+1),
190//         maxnorm(g)*power(CanonicalForm(2),g.degree())*isqrt(g.degree()+1)
191//         )+1;
192    M = tmin( maxNorm(f), maxNorm(g) );
193    BB = power(CanonicalForm(2),tmin(f.degree(),g.degree()))*M;
194    q = 0;
195    i = cf_getNumSmallPrimes() - 1;
196    while ( true ) {
197        B = BB;
198        while ( i >= 0 && q < B ) {
199            p = cf_getSmallPrime( i );
200            i--;
201            while ( i >= 0 && mod( cl, p ) == 0 ) {
202                p = cf_getSmallPrime( i );
203                i--;
204            }
205            setCharacteristic( p );
206            Dp = gcd( mapinto( f ), mapinto( g ) );
207            Dp = ( Dp / Dp.lc() ) * mapinto( cg );
208            setCharacteristic( 0 );
209            if ( Dp.degree() == 0 ) return c;
210            if ( q.isZero() ) {
211                D = mapinto( Dp );
212                q = p;
213                B = power(CanonicalForm(2),D.degree())*M+1;
214            }
215            else {
216                if ( Dp.degree() == D.degree() ) {
217                    chineseRemainder( D, q, mapinto( Dp ), p, newD, newq );
218                    q = newq;
219                    D = newD;
220                }
221                else if ( Dp.degree() < D.degree() ) {
222                    // all previous p's are bad primes
223                    q = p;
224                    D = mapinto( Dp );
225                    B = power(CanonicalForm(2),D.degree())*M+1;
226                }
227                // else p is a bad prime
228            }
229        }
230        if ( i >= 0 ) {
231            // now balance D mod q
232            D = pp( balance( D, q ) );
233            if ( divides( D, f ) && divides( D, g ) )
234                return D * c;
235            else
236                q = 0;
237        }
238        else {
239            return gcd_poly( F, G, false );
240        }
241        DEBOUTLN( cerr, "another try ..." );
242    }
243}
244
245static CanonicalForm
246gcd_poly1( const CanonicalForm & f, const CanonicalForm & g, bool modularflag )
247{
248    CanonicalForm C, Ci, Ci1, Hi, bi, pi, pi1, pi2;
249    int delta;
250    Variable v = f.mvar();
251
252    if ( f.degree( v ) >= g.degree( v ) ) {
253        pi = f; pi1 = g;
254    }
255    else {
256        pi = g; pi1 = f;
257    }
258    Ci = content( pi ); Ci1 = content( pi1 );
259    C = gcd( Ci, Ci1 );
260    pi1 = pi1 / Ci1; pi = pi / Ci;
261    if ( pi.isUnivariate() && pi1.isUnivariate() ) {
262        if ( modularflag )
263            return gcd_poly_univar0( pi, pi1, true ) * C;
264    }
265    else
266        if ( gcd_test_one( pi1, pi, true ) )
267            return C;
268    delta = degree( pi, v ) - degree( pi1, v );
269    Hi = power( LC( pi1, v ), delta );
270    if ( (delta+1) % 2 )
271        bi = 1;
272    else
273        bi = -1;
274    while ( degree( pi1, v ) > 0 ) {
275        pi2 = psr( pi, pi1, v );
276        pi2 = pi2 / bi;
277        pi = pi1; pi1 = pi2;
278        if ( degree( pi1, v ) > 0 ) {
279            delta = degree( pi, v ) - degree( pi1, v );
280            if ( (delta+1) % 2 )
281                bi = LC( pi, v ) * power( Hi, delta );
282            else
283                bi = -LC( pi, v ) * power( Hi, delta );
284            Hi = power( LC( pi1, v ), delta ) / power( Hi, delta-1 );
285        }
286    }
287    if ( degree( pi1, v ) == 0 )
288        return C;
289    else {
290        return C * pp( pi );
291    }
292}
293
294//{{{ static CanonicalForm gcd_poly ( const CanonicalForm & f, const CanonicalForm & g, bool modularflag )
295//{{{ docu
296//
297// gcd_poly() - calculate polynomial gcd.
298//
299// This is the dispatcher for polynomial gcd calculation.  We call either
300// ezgcd(), sparsemod() or gcd_poly1() in dependecy on the current
301// characteristic and settings of SW_USE_EZGCD and SW_USE_SPARSEMOD, resp.
302//
303// modularflag is reached down to gcd_poly1() without change in case of
304// zero characteristic.
305//
306// Used by gcd() and gcd_poly_univar0().
307//
308//}}}
309static CanonicalForm
310gcd_poly ( const CanonicalForm & f, const CanonicalForm & g, bool modularflag )
311{
312    if ( getCharacteristic() != 0 ) {
313        return gcd_poly1( f, g, false );
314    }
315    else if ( isOn( SW_USE_EZGCD ) && ! ( f.isUnivariate() && g.isUnivariate() ) ) {
316        CFMap M, N;
317        compress( f, g, M, N );
318        return N( ezgcd( M(f), M(g) ) );
319    }
320    else if ( isOn( SW_USE_SPARSEMOD ) && ! ( f.isUnivariate() && g.isUnivariate() ) ) {
321        return sparsemod( f, g );
322    }
323    else {
324        return gcd_poly1( f, g, modularflag );
325    }
326}
327//}}}
328
329//{{{ static CanonicalForm cf_content ( const CanonicalForm & f, const CanonicalForm & g )
330//{{{ docu
331//
332// cf_content() - return gcd(g, content(f)).
333//
334// content(f) is calculated with respect to f's main variable.
335//
336// Used by gcd(), content(), content( CF, Variable ).
337//
338//}}}
339static CanonicalForm
340cf_content ( const CanonicalForm & f, const CanonicalForm & g )
341{
342    if ( f.inPolyDomain() || ( f.inExtension() && ! getReduce( f.mvar() ) ) ) {
343        CFIterator i = f;
344        CanonicalForm result = g;
345        while ( i.hasTerms() && ! result.isOne() ) {
346            result = gcd( result, i.coeff() );
347            i++;
348        }
349        return result;
350    }
351    else
352        if ( f.sign() < 0 )
353            return -f;
354        else
355            return f;
356}
357//}}}
358
359//{{{ CanonicalForm content ( const CanonicalForm & f )
360//{{{ docu
361//
362// content() - return content(f) with respect to main variable.
363//
364// Normalizes result.
365//
366//}}}
367CanonicalForm
368content ( const CanonicalForm & f )
369{
370    return cf_content( f, 0 );
371}
372//}}}
373
374//{{{ CanonicalForm content ( const CanonicalForm & f, const Variable & x )
375//{{{ docu
376//
377// content() - return content(f) with respect to x.
378//
379// x should be a polynomial variable.
380//
381//}}}
382CanonicalForm
383content ( const CanonicalForm & f, const Variable & x )
384{
385    ASSERT( x.level() > 0, "cannot calculate content with respect to algebraic variable" );
386    Variable y = f.mvar();
387
388    if ( y == x )
389        return cf_content( f, 0 );
390    else  if ( y < x )
391        return f;
392    else
393        return swapvar( content( swapvar( f, y, x ), y ), y, x );
394}
395//}}}
396
397//{{{ CanonicalForm vcontent ( const CanonicalForm & f, const Variable & x )
398//{{{ docu
399//
400// vcontent() - return content of f with repect to variables >= x.
401//
402// The content is recursively calculated over all coefficients in
403// f having level less than x.  x should be a polynomial
404// variable.
405//
406//}}}
407CanonicalForm
408vcontent ( const CanonicalForm & f, const Variable & x )
409{
410    ASSERT( x.level() > 0, "cannot calculate vcontent with respect to algebraic variable" );
411
412    if ( f.mvar() <= x )
413        return content( f, x );
414    else {
415        CFIterator i;
416        CanonicalForm d = 0;
417        for ( i = f; i.hasTerms() && ! d.isOne(); i++ )
418            d = gcd( d, vcontent( i.coeff(), x ) );
419        return d;
420    }
421}
422//}}}
423
424//{{{ CanonicalForm pp ( const CanonicalForm & f )
425//{{{ docu
426//
427// pp() - return primitive part of f.
428//
429// Returns zero if f equals zero, otherwise f / content(f).
430//
431//}}}
432CanonicalForm
433pp ( const CanonicalForm & f )
434{
435    if ( f.isZero() )
436        return f;
437    else
438        return f / content( f );
439}
440//}}}
441
442CanonicalForm
443gcd ( const CanonicalForm & f, const CanonicalForm & g )
444{
445    if ( f.isZero() )
446        if ( g.lc().sign() < 0 )
447            return -g;
448        else
449            return g;
450    else  if ( g.isZero() )
451        if ( f.lc().sign() < 0 )
452            return -f;
453        else
454            return f;
455    else  if ( f.inBaseDomain() )
456        if ( g.inBaseDomain() )
457            return bgcd( f, g );
458        else
459            return cf_content( g, f );
460    else  if ( g.inBaseDomain() )
461        return cf_content( f, g );
462    else  if ( f.mvar() == g.mvar() )
463        if ( f.inExtension() && getReduce( f.mvar() ) )
464            return 1;
465        else {
466            if ( divides( f, g ) )
467                if ( f.lc().sign() < 0 )
468                    return -f;
469                else
470                    return f;
471            else  if ( divides( g, f ) )
472                if ( g.lc().sign() < 0 )
473                    return -g;
474                else
475                    return g;
476            if ( getCharacteristic() == 0 && isOn( SW_RATIONAL ) ) {
477                CanonicalForm cdF = bCommonDen( f );
478                CanonicalForm cdG = bCommonDen( g );
479                Off( SW_RATIONAL );
480                CanonicalForm l = lcm( cdF, cdG );
481                On( SW_RATIONAL );
482                CanonicalForm F = f * l, G = g * l;
483                Off( SW_RATIONAL );
484                l = gcd_poly( F, G, true );
485                On( SW_RATIONAL );
486                if ( l.lc().sign() < 0 )
487                    return -l;
488                else
489                    return l;
490            }
491            else {
492                CanonicalForm d = gcd_poly( f, g, getCharacteristic()==0 );
493                if ( d.lc().sign() < 0 )
494                    return -d;
495                else
496                    return d;
497            }
498        }
499    else  if ( f.mvar() > g.mvar() )
500        return cf_content( f, g );
501    else
502        return cf_content( g, f );
503}
504
505//{{{ CanonicalForm lcm ( const CanonicalForm & f, const CanonicalForm & g )
506//{{{ docu
507//
508// lcm() - return least common multiple of f and g.
509//
510// The lcm is calculated using the formula lcm(f, g) = f * g / gcd(f, g).
511//
512// Returns zero if one of f or g equals zero.
513//
514//}}}
515CanonicalForm
516lcm ( const CanonicalForm & f, const CanonicalForm & g )
517{
518    if ( f.isZero() || g.isZero() )
519        return f;
520    else
521        return ( f / gcd( f, g ) ) * g;
522}
523//}}}
Note: See TracBrowser for help on using the repository browser.