source: git/factory/fac_ezgcd.cc @ 19fc57b

fieker-DuValspielwiese
Last change on this file since 19fc57b was 2667bc8, checked in by Hans Schönemann <hannes@…>, 19 years ago
*hannes: git-svn-id: file:///usr/local/Singular/svn/trunk@8012 2c84dea3-7e68-4137-9b89-c4e89433aadc
  • Property mode set to 100644
File size: 12.0 KB
Line 
1/* emacs edit mode for this file is -*- C++ -*- */
2/* $Id: fac_ezgcd.cc,v 1.19 2005-05-03 09:35:34 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 "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, 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    {
57        DEBDECLEVEL( cerr, "ezgcd" );
58        return d * gcd( F, G );
59    }
60    else  if ( gcd_test_one( F, G, false ) )
61    {
62        DEBDECLEVEL( cerr, "ezgcd" );
63        return d;
64    }
65    lcF = LC( F, x ); lcG = LC( G, x );
66    lcD = gcd( lcF, lcG );
67    delta = 0;
68    degF = degree( F, x ); degG = degree( G, x );
69    t = tmax( F.level(), G.level() );
70    bound = findBound( F, G, lcF, lcG, degF, degG );
71    if ( ! internal )
72        b = REvaluation( 2, t, IntRandom( 50 ) );
73    while ( ! gcdfound ) {
74        /// ---> A2
75        DEBOUTLN( cerr, "search for evaluation, delta = " << delta );
76        DEBOUTLN( cerr, "F = " << F );
77        DEBOUTLN( cerr, "G = " << G );
78        findeval( F, G, Fb, Gb, Db, b, delta, degF, degG );
79        DEBOUTLN( cerr, "found evaluation b = " << b );
80        DEBOUTLN( cerr, "F(b) = " << Fb );
81        DEBOUTLN( cerr, "G(b) = " << Gb );
82        DEBOUTLN( cerr, "D(b) = " << Db );
83        delta = degree( Db );
84        /// ---> A3
85        if ( delta == 0 )
86        {
87          DEBDECLEVEL( cerr, "ezgcd" );
88          return d;
89        }
90        /// ---> A4
91        //deltaold = delta;
92        while ( 1 ) {
93            bt = b;
94            findeval( F, G, Fbt, Gbt, Dbt, bt, delta + 1, degF, degG );
95            int dd=degree( Dbt );
96            if ( dd /*degree( Dbt )*/ == 0 )
97            {
98                DEBDECLEVEL( cerr, "ezgcd" );
99                return d;
100            }
101            if ( dd /*degree( Dbt )*/ == delta )
102                break;
103            else  if ( dd /*degree( Dbt )*/ < delta ) {
104                delta = dd /*degree( Dbt )*/;
105                b = bt;
106                Db = Dbt; Fb = Fbt; Gb = Gbt;
107            }
108        }
109        DEBOUTLN( cerr, "now after A4, delta = " << delta );
110        /// ---> A5
111        if ( degF <= degG && delta == degF && divides( F, G ) )
112        {
113            DEBDECLEVEL( cerr, "ezgcd" );
114            return d*F;
115        }
116        if ( degG < degF && delta == degG && divides( G, F ) )
117        {
118            DEBDECLEVEL( cerr, "ezgcd" );
119            return d*G;
120        }
121        if ( delta != degF && delta != degG ) {
122            /// ---> A6
123            CanonicalForm xxx;
124            //if ( gcd( (DD[1] = Fb / Db), Db ) == 1 ) {
125            DD[1] = Fb / Db;
126            xxx= gcd( DD[1], Db );
127            DEBOUTLN( cerr, "gcd((Fb/Db),Db) = " << xxx );
128            DEBOUTLN( cerr, "Fb/Db = " << DD[1] );
129            DEBOUTLN( cerr, "Db = " << Db );
130            if (xxx.inCoeffDomain()) {
131                B = F;
132                DD[2] = Db;
133                lcDD[1] = lcF;
134                lcDD[2] = lcF;
135                B *= lcF;
136            }
137            //else  if ( gcd( (DD[1] = Gb / Db), Db ) == 1 ) {
138            else
139            { 
140              DD[1] = Gb / Db;
141              xxx=gcd( DD[1], Db );
142              DEBOUTLN( cerr, "gcd((Gb/Db),Db) = " << xxx );
143              DEBOUTLN( cerr, "Gb/Db = " << DD[1] );
144              DEBOUTLN( cerr, "Db = " << Db );
145              if (xxx.inCoeffDomain())
146              {
147                B = G;
148                DD[2] = Db;
149                lcDD[1] = lcG;
150                lcDD[2] = lcG;
151                B *= lcG;
152              }
153              else
154              {
155#ifdef DEBUGOUTPUT
156                CanonicalForm dummyres = d * ezgcd_specialcase( F, G, b, bound );
157                DEBDECLEVEL( cerr, "ezgcd" );
158                return dummyres;
159#else
160                return d * ezgcd_specialcase( F, G, b, bound );
161#endif
162              }
163            }
164            /// ---> A7
165            DD[2] = DD[2] * ( b( lcDD[2] ) / lc( DD[2] ) );
166            DD[1] = DD[1] * ( b( lcDD[1] ) / lc( DD[1] ) );
167            bound = enlargeBound( B, DD[2], DD[1], bound );
168            DEBOUTLN( cerr, "(hensel) B    = " << B );
169            DEBOUTLN( cerr, "(hensel) lcB  = " << LC( B, Variable(1) ) );
170            DEBOUTLN( cerr, "(hensel) b(B) = " << b(B) );
171            DEBOUTLN( cerr, "(hensel) DD   = " << DD );
172            DEBOUTLN( cerr, "(hensel) lcDD = " << lcDD );
173            gcdfound = Hensel( B, DD, lcDD, b, bound, x );
174            DEBOUTLN( cerr, "(hensel finished) DD   = " << DD );
175           
176            if (gcdfound)
177            {
178              CanonicalForm cand=DD[2] / content(DD[2],Variable(1));
179              if (B==F)
180              {
181                DEBOUTLN( cerr, "(test) G: "<<G<<" % gcd:"<<cand<<" -> " << G%cand );
182                gcdfound= divides(cand,G);
183              }
184              else
185              { 
186                DEBOUTLN( cerr, "(test) F: "<<F<<" % gcd:"<<cand<<" -> " << F%cand);
187                gcdfound= divides(cand,F);
188              }
189            }
190            /// ---> A8 (gcdfound)
191        }
192        delta++;
193    }
194    /// ---> A9
195    DEBDECLEVEL( cerr, "ezgcd" );
196    return d * DD[2] / content(DD[2],Variable(1));
197}
198
199static CanonicalForm
200ezgcd_specialcase ( const CanonicalForm & F, const CanonicalForm & G, REvaluation & b, const modpk & bound )
201{
202    DEBOUTLN( cerr, "ezgcd: special case" );
203    CanonicalForm Ft, Gt, L, LL, Fb, Gb, Db, Lb, D, Ds, Dt, d;
204    CFArray DD( 1, 2 ), lcDD( 1, 2 );
205    Variable x = Variable( 1 );
206    bool gcdfound;
207
208    Dt = 1;
209    /// ---> S1
210    DEBOUTLN( cerr, "ezgcdspec: (S1)" );
211    Ft = ezgcd( F, F.deriv( x ) );
212    if ( Ft.isOne() ) {
213        // In this case F is squarefree and we came here by bad chance
214        // (means: bad evaluation point).  Try again with another
215        // evaluation point.  Bug fix (?) by JS.  The bad example was
216        // gcd.debug -ocr /+USE_EZGCD/@12/CB \
217        //     '(16*B^8-208*B^6*C+927*B^4*C^2-1512*B^2*C^3+432*C^4)' \
218        //     '(4*B^7*C^2-50*B^5*C^3+208*B^3*C^4-288*B*C^5)'
219        b.nextpoint();
220        return ezgcd( F, G, b, true );
221    }
222#if 1
223    Off(SW_USE_EZGCD);
224    //bool ntl0=isOn(SW_USE_NTL_GCD_0);
225    //Off(SW_USE_NTL_GCD_0);
226    //bool ntlp=isOn(SW_USE_NTL_GCD_P);
227    //Off(SW_USE_NTL_GCD_P);
228    d=gcd( F, G );
229    //if (ntl0) On(SW_USE_NTL_GCD_0);
230    //if (ntlp) On(SW_USE_NTL_GCD_P);
231    On(SW_USE_EZGCD);
232    return d;
233#else
234
235    DEBOUTLN( cerr, "ezgcdspec: (S1) done, Ft = " << Ft );
236    L = F / Ft;
237    /// ---> S2
238    DEBOUTLN( cerr, "ezgcdspec: (S2)" );
239
240    L = ezgcd( L, G, b, true );
241    DEBOUTLN( cerr, "ezgcdspec: (S2) done, Ds = " << Ds );
242    Gt = G / L;
243    Ds = L; Dt = L;
244    Lb = b( L ); Gb = b( Gt ); Fb = b( Ft );
245    d = gcd( LC( L, x ), gcd( LC( Ft, x ), LC( Gt, x ) ) );
246    while ( true ) {
247        /// ---> S3
248        DEBOUTLN( cerr, "ezgcdspec: (S3)" );
249        DEBOUTLN( cerr, "ezgcdspec: Fb = " << Fb );
250        DEBOUTLN( cerr, "ezgcdspec: Gb = " << Gb );
251        DD[1] = gcd( Lb, gcd( Fb, Gb ) );
252        /// ---> S4
253        DEBOUTLN( cerr, "ezgcdspec: (S4)" );
254        if ( degree( DD[1] ) == 0 )
255            return Ds;
256        /// ---> S5
257        DEBOUTLN( cerr, "ezgcdspec: (S5)" );
258        DEBOUTLN( cerr, "ezgcdspec: Lb = " << Lb );
259        DEBOUTLN( cerr, "ezgcdspec: Db = " << DD[1] );
260        Db = DD[1];
261        if ( ! (DD[2] = Lb/DD[1]).isOne() ) {
262            DEBOUTLN( cerr, "ezgcdspec: (S55)" );
263            lcDD[2] = LC( L, x );
264            lcDD[1] = d;
265            DD[1] = ( DD[1] * b( d ) ) / lc( DD[1] );
266            DD[2] = ( DD[2] * b( lcDD[2] ) ) / lc( DD[2] );
267            LL = L * d;
268            modpk newbound = enlargeBound( LL, DD[2], DD[1], bound );
269            DEBOUTLN( cerr, "ezgcdspec: begin with lift." );
270            DEBOUTLN( cerr, "ezgcdspec: B     = " << LL );
271            DEBOUTLN( cerr, "ezgcdspec: DD    = " << DD );
272            DEBOUTLN( cerr, "ezgcdspec: lcDD  = " << lcDD );
273            DEBOUTLN( cerr, "ezgcdspec: b     = " << b );
274            DEBOUTLN( cerr, "ezgcdspec: bound = " << bound.getpk() );
275            DEBOUTLN( cerr, "ezgcdspec: lc(B) = " << LC( LL, x ) );
276            DEBOUTLN( cerr, "ezgcdspec: test  = " << b( LL ) - DD[1] * DD[2] );
277            gcdfound = Hensel( LL, DD, lcDD, b, newbound, x );
278            ASSERT( gcdfound, "fatal error in algorithm" );
279            Dt = pp( DD[1] );
280        }
281        DEBOUTLN( cerr, "ezgcdspec: (S7)" );
282        Ds = Ds * Dt;
283        Fb = Fb / Db;
284        Gb = Gb / Db;
285    }
286    // this point is never reached, but to avoid compiler warnings let's return a value
287    return 0;
288#endif
289}
290
291static void
292findeval( const CanonicalForm & F, const CanonicalForm & G, CanonicalForm & Fb, CanonicalForm & Gb, CanonicalForm & Db, REvaluation & b, int delta, int degF, int degG )
293{
294    int i;
295    bool ok;
296    if ( delta != 0 )
297        b.nextpoint();
298    DEBOUTLN( cerr, "ezgcd: (findeval) F = " << F  <<", G="<< G);
299    DEBOUTLN( cerr, "ezgcd: (findeval) degF = " << degF << ", degG="<<degG );
300    do {
301        DEBOUTLN( cerr, "ezgcd: (findeval) b = " << b );
302        Fb = b( F );
303        ok = degree( Fb ) == degF;
304        if ( ok ) {
305            Gb = b( G );
306            ok = degree( Gb ) == degG;
307        }
308       
309        if ( ok )
310        {   
311            Db = gcd( Fb, Gb );
312            if ( delta > 0 )
313              ok = degree( Db ) < delta;
314        }
315        if ( ! ok )
316            b.nextpoint();
317    } while ( ! ok );
318}
319
320static modpk
321enlargeBound ( const CanonicalForm & F, const CanonicalForm & Lb, const CanonicalForm & Db, const modpk & pk )
322{
323    DEBOUTLN( cerr, "ezgcd: (enlarge bound) F      = " << F );
324    DEBOUTLN( cerr, "ezgcd: (enlarge bound) Lb     = " << Lb );
325    DEBOUTLN( cerr, "ezgcd: (enlarge bound) Db     = " << Db );
326    DEBOUTLN( cerr, "ezgcd: (enlarge bound) Lb % p = " << mod( Lb, pk.getp() ) );
327    DEBOUTLN( cerr, "ezgcd: (enlarge bound) Db % p = " << mod( Db, pk.getp() ) );
328
329    CanonicalForm limit = power( CanonicalForm(2), degree( Db ) ) *
330        tmax( maxNorm( Lb ), tmax( maxNorm( Db ), maxNorm( F ) ) );
331    int p = pk.getp();
332    int k = pk.getk();
333    CanonicalForm bound = pk.getpk();
334    while ( bound < limit ) {
335        k++;
336        bound *= p;
337    }
338    k *= 2;
339    DEBOUTLN( cerr, "ezgcd: (enlarge bound) newbound = " << power( CanonicalForm( p ), k ) );
340    return modpk( p, k );
341}
342
343static modpk
344findBound ( const CanonicalForm & F, const CanonicalForm & G, const CanonicalForm & lcF, const CanonicalForm & lcG, int degF, int degG )
345{
346    CanonicalForm limit = power( CanonicalForm(2), tmin( degF, degG ) ) *
347        gcd( icontent( lcF ), icontent( lcG ) ) * tmin( maxNorm( F ), maxNorm( G ) );
348    int p, i = 0;
349    do {
350        p = cf_getBigPrime( i );
351        i++;
352    } while ( mod( lcF, p ).isZero() && mod( lcG, p ).isZero() );
353    CanonicalForm bound = p;
354    i = 1;
355    while ( bound < limit ) {
356        i++;
357        bound *= p;
358    }
359    return modpk( p, i );
360}
Note: See TracBrowser for help on using the repository browser.