source: git/factory/fac_ezgcd.cc @ 76a7114

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