source: git/factory/fac_ezgcd.cc @ b973c0

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