source: git/factory/fac_ezgcd.cc @ 597bd3

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