source: git/factory/fac_ezgcd.cc @ 72dd6e

spielwiese
Last change on this file since 72dd6e was 2c4671, checked in by Jens Schmidt <schmidt@…>, 27 years ago
* fac_ezgcd.cc: include paths for mac added git-svn-id: file:///usr/local/Singular/svn/trunk@449 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 8.8 KB
Line 
1/* emacs edit mode for this file is -*- C++ -*- */
2/* $Id: fac_ezgcd.cc,v 1.7 1997-06-30 15:26:31 schmidt Exp $ */
3
4#include <config.h>
5
6#include "assert.h"
7#include "debug.h"
8
9#include "cf_defs.h"
10#ifdef macintosh
11#include <::templates:ftmpl_functions.h>
12#else
13#include "templates/ftmpl_functions.h"
14#endif
15#include "canonicalform.h"
16#include "fac_util.h"
17#include "cf_reval.h"
18#include "cf_random.h"
19#include "cf_primes.h"
20#include "fac_distrib.h"
21
22
23static void findeval( const CanonicalForm & F, const CanonicalForm & G, CanonicalForm & Fb, CanonicalForm & Gb, CanonicalForm & Db, REvaluation & b, int delta, int degF, int degG );
24
25static CanonicalForm ezgcd ( const CanonicalForm & FF, const CanonicalForm & GG, REvaluation & b, bool internal );
26
27static CanonicalForm ezgcd_specialcase ( const CanonicalForm & F, const CanonicalForm & G, REvaluation & b, const modpk & bound );
28
29static modpk findBound ( const CanonicalForm & F, const CanonicalForm & G, const CanonicalForm & lcF, const CanonicalForm & lcG, int degF, int degG );
30
31static modpk enlargeBound ( const CanonicalForm & F, const CanonicalForm & Lb, const CanonicalForm & Db, const modpk & pk );
32
33CanonicalForm
34ezgcd ( const CanonicalForm & FF, const CanonicalForm & GG )
35{
36    REvaluation b;
37    return ezgcd( FF, GG, b, false );
38}
39
40static CanonicalForm
41ezgcd ( const CanonicalForm & FF, const CanonicalForm & GG, REvaluation & b, bool internal )
42{
43    CanonicalForm F, G, f, g, d, Fb, Gb, Db, Fbt, Gbt, Dbt, B0, B, D0, lcF, lcG, lcD;
44    CFArray DD( 1, 2 ), lcDD( 1, 2 );
45    int degF, degG, delta, deltaold, t;
46    REvaluation bt;
47    bool gcdfound = false;
48    Variable x = Variable(1);
49    modpk bound;
50
51    DEBINCLEVEL( cerr, "ezgcd" );
52    DEBOUTLN( cerr, "FF = " << FF );
53    DEBOUTLN( cerr, "GG = " << GG );
54    f = content( FF, x ); g = content( GG, x ); d = gcd( f, g );
55    DEBOUTLN( cerr, "f = " << f );
56    DEBOUTLN( cerr, "g = " << g );
57    F = FF / f; G = GG / g;
58    if ( F.isUnivariate() && G.isUnivariate() )
59        return d * gcd( F, G );
60    else  if ( gcd_test_one( F, G, false ) )
61        return d;
62    lcF = LC( F, x ); lcG = LC( G, x );
63    lcD = gcd( lcF, lcG );
64    delta = 0;
65    degF = degree( F, x ); degG = degree( G, x );
66    t = tmax( F.level(), G.level() );
67    bound = findBound( F, G, lcF, lcG, degF, degG );
68    if ( ! internal )
69        b = REvaluation( 2, t, IntRandom( 100 ) );
70    while ( ! gcdfound ) {
71        /// ---> A2
72        DEBOUTLN( cerr, "search for evaluation, delta = " << delta );
73        DEBOUTLN( cerr, "F = " << F );
74        DEBOUTLN( cerr, "G = " << G );
75        findeval( F, G, Fb, Gb, Db, b, delta, degF, degG );
76        DEBOUTLN( cerr, "found evaluation b = " << b );
77        DEBOUTLN( cerr, "F(b) = " << Fb );
78        DEBOUTLN( cerr, "G(b) = " << Gb );
79        DEBOUTLN( cerr, "D(b) = " << Db );
80        delta = degree( Db );
81        /// ---> A3
82        if ( delta == 0 )
83            return d;
84        /// ---> A4
85        deltaold = delta + 1;
86        while ( deltaold != delta ) {
87            bt = b;
88            findeval( F, G, Fbt, Gbt, Dbt, bt, delta + 1, degF, degG );
89            if ( degree( Dbt ) == 0 )
90                return d;
91            if ( degree( Dbt ) == delta )
92                deltaold = delta;
93            else  if ( degree( Dbt ) < delta ) {
94                deltaold = delta;
95                delta = degree( Dbt );
96                b = bt;
97                Db = Dbt; Fb = Fbt; Gb = Gbt;
98            }
99        }
100        DEBOUTLN( cerr, "now after A4, delta = " << delta );
101        /// ---> A5
102        if ( degF <= degG && delta == degF && divides( F, G ) )
103            return d*F;
104        if ( degG < degF && delta == degG && divides( G, F ) )
105            return d*G;
106        if ( delta != degF && delta != degG ) {
107            /// ---> A6
108            CanonicalForm xxx;
109            if ( gcd( (DD[1] = Fb / Db), Db ) == 1 ) {
110                B = F;
111                DD[2] = Db;
112                lcDD[1] = lcF;
113                lcDD[2] = lcF;
114                B *= lcF;
115            }
116            else  if ( gcd( (DD[1] = Gb / Db), Db ) == 1 ) {
117                B = G;
118                DD[2] = Db;
119                lcDD[1] = lcG;
120                lcDD[2] = lcG;
121                B *= lcG;
122            }
123            else {
124#ifdef DEBUGOUTPUT
125                CanonicalForm dummyres = d * ezgcd_specialcase( F, G, b, bound );
126                DEBDECLEVEL( cerr, "ezgcd" );
127                return dummyres;
128#else
129                return d * ezgcd_specialcase( F, G, b, bound );
130#endif
131            }
132            /// ---> A7
133            DD[2] = DD[2] * ( b( lcDD[2] ) / lc( DD[2] ) );
134            DD[1] = DD[1] * ( b( lcDD[1] ) / lc( DD[1] ) );
135            bound = enlargeBound( B, DD[2], DD[1], bound );
136            DEBOUTLN( cerr, "(hensel) B    = " << B );
137            DEBOUTLN( cerr, "(hensel) lcB  = " << LC( B, Variable(1) ) );
138            DEBOUTLN( cerr, "(hensel) b(B) = " << b(B) );
139            DEBOUTLN( cerr, "(hensel) DD   = " << DD );
140            DEBOUTLN( cerr, "(hensel) lcDD = " << lcDD );
141            gcdfound = Hensel( B, DD, lcDD, b, bound, x );
142            DEBOUTLN( cerr, "(hensel finished) DD   = " << DD );
143            /// ---> A8 (gcdfound)
144        }
145    }
146    /// ---> A9
147    DEBDECLEVEL( cerr, "ezgcd" );
148    return d * DD[2] / content(DD[2],Variable(1));
149}
150
151static CanonicalForm
152ezgcd_specialcase ( const CanonicalForm & F, const CanonicalForm & G, REvaluation & b, const modpk & bound )
153{
154    DEBOUTLN( cerr, "ezgcd: special case" );
155    CanonicalForm Ft, Gt, L, LL, Fb, Gb, Db, Lb, D, Ds, Dt, d;
156    CFArray DD( 1, 2 ), lcDD( 1, 2 );
157    Variable x = Variable( 1 );
158    bool gcdfound;
159
160    Dt = 1;
161    /// ---> S1
162    DEBOUTLN( cerr, "ezgcdspec: (S1)" );
163    Ft = ezgcd( F, F.deriv( x ) );
164    DEBOUTLN( cerr, "ezgcdspec: (S1) done, Ft = " << Ft );
165    L = F / Ft;
166    /// ---> S2
167    DEBOUTLN( cerr, "ezgcdspec: (S2)" );
168    L = ezgcd( L, G, b, true );
169    DEBOUTLN( cerr, "ezgcdspec: (S2) done, Ds = " << Ds );
170    Gt = G / L;
171    Ds = L; Dt = L;
172    Lb = b( L ); Gb = b( Gt ); Fb = b( Ft );
173    d = gcd( LC( L, x ), gcd( LC( Ft, x ), LC( Gt, x ) ) );
174    while ( true ) {
175        /// ---> S3
176        DEBOUTLN( cerr, "ezgcdspec: (S3)" );
177        DEBOUTLN( cerr, "ezgcdspec: Fb = " << Fb );
178        DEBOUTLN( cerr, "ezgcdspec: Gb = " << Gb );
179        DD[1] = gcd( Lb, gcd( Fb, Gb ) );
180        /// ---> S4
181        DEBOUTLN( cerr, "ezgcdspec: (S4)" );
182        if ( degree( DD[1] ) == 0 )
183            return Ds;
184        /// ---> S5
185        DEBOUTLN( cerr, "ezgcdspec: (S5)" );
186        DEBOUTLN( cerr, "ezgcdspec: Lb = " << Lb );
187        DEBOUTLN( cerr, "ezgcdspec: Db = " << DD[1] );
188        Db = DD[1];
189        if ( ! (DD[2] = Lb/DD[1]).isOne() ) {
190            DEBOUTLN( cerr, "ezgcdspec: (S55)" );
191            lcDD[2] = LC( L, x );
192            lcDD[1] = d;
193            DD[1] = ( DD[1] * b( d ) ) / lc( DD[1] );
194            DD[2] = ( DD[2] * b( lcDD[2] ) ) / lc( DD[2] );
195            LL = L * d;
196            modpk newbound = enlargeBound( LL, DD[2], DD[1], bound );
197            DEBOUTLN( cerr, "ezgcdspec: begin with lift." );
198            DEBOUTLN( cerr, "ezgcdspec: B     = " << LL );
199            DEBOUTLN( cerr, "ezgcdspec: DD    = " << DD );
200            DEBOUTLN( cerr, "ezgcdspec: lcDD  = " << lcDD );
201            DEBOUTLN( cerr, "ezgcdspec: b     = " << b );
202            DEBOUTLN( cerr, "ezgcdspec: bound = " << bound.getpk() );
203            DEBOUTLN( cerr, "ezgcdspec: lc(B) = " << LC( LL, x ) );
204            DEBOUTLN( cerr, "ezgcdspec: test  = " << b( LL ) - DD[1] * DD[2] );
205            gcdfound = Hensel( LL, DD, lcDD, b, newbound, x );
206            ASSERT( gcdfound, "fatal error in algorithm" );
207            Dt = pp( DD[1] );
208        }
209        DEBOUTLN( cerr, "ezgcdspec: (S7)" );
210        Ds = Ds * Dt;
211        Fb = Fb / Db;
212        Gb = Gb / Db;
213    }
214    // this point is never reached, but to avoid compiler warnings let's return a value
215    return 0;
216}
217
218static void
219findeval( const CanonicalForm & F, const CanonicalForm & G, CanonicalForm & Fb, CanonicalForm & Gb, CanonicalForm & Db, REvaluation & b, int delta, int degF, int degG )
220{
221    bool ok;
222    if ( delta != 0 )
223        b.nextpoint();
224    do {
225        DEBOUTLN( cerr, "ezgcd: (findeval) b = " << b );
226        Fb = b( F );
227        ok = degree( Fb ) == degF;
228        if ( ok ) {
229            Gb = b( G );
230            ok = degree( Gb ) == degG;
231        }
232        if ( ok ) {
233            Db = gcd( Fb, Gb );
234            if ( delta > 0 )
235                ok = degree( Db ) < delta;
236        }
237        if ( ! ok )
238            b.nextpoint();
239    } while ( ! ok );
240}
241
242static modpk
243enlargeBound ( const CanonicalForm & F, const CanonicalForm & Lb, const CanonicalForm & Db, const modpk & pk )
244{
245    DEBOUTLN( cerr, "ezgcd: (enlarge bound) F      = " << F );
246    DEBOUTLN( cerr, "ezgcd: (enlarge bound) Lb     = " << Lb );
247    DEBOUTLN( cerr, "ezgcd: (enlarge bound) Db     = " << Db );
248    DEBOUTLN( cerr, "ezgcd: (enlarge bound) Lb % p = " << mod( Lb, pk.getp() ) );
249    DEBOUTLN( cerr, "ezgcd: (enlarge bound) Db % p = " << mod( Db, pk.getp() ) );
250
251    CanonicalForm limit = power( CanonicalForm(2), degree( Db ) ) *
252        tmax( maxCoeff( Lb ), tmax( maxCoeff( Db ), maxCoeff( F ) ) );
253    int p = pk.getp();
254    int k = pk.getk();
255    CanonicalForm bound = pk.getpk();
256    while ( bound < limit ) {
257        k++;
258        bound *= p;
259    }
260    k *= 2;
261    DEBOUTLN( cerr, "ezgcd: (enlarge bound) newbound = " << power( CanonicalForm( p ), k ) );
262    return modpk( p, k );
263}
264
265static modpk
266findBound ( const CanonicalForm & F, const CanonicalForm & G, const CanonicalForm & lcF, const CanonicalForm & lcG, int degF, int degG )
267{
268    CanonicalForm limit = power( CanonicalForm(2), tmin( degF, degG ) ) *
269        gcd( icontent( lcF ), icontent( lcG ) ) * tmin( maxCoeff( F ), maxCoeff( G ) );
270    int p, i = 0;
271    do {
272        p = cf_getBigPrime( i );
273        i++;
274    } while ( mod( lcF, p ).isZero() && mod( lcG, p ).isZero() );
275    CanonicalForm bound = p;
276    i = 1;
277    while ( bound < limit ) {
278        i++;
279        bound *= p;
280    }
281    return modpk( p, i );
282}
Note: See TracBrowser for help on using the repository browser.