source: git/factory/fac_ezgcd.cc @ a70441f

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