source: git/factory/cf_gcd.cc @ 71da5e

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