source: git/factory/fac_ezgcd.cc @ f7a4e9

spielwiese
Last change on this file since f7a4e9 was f7a4e9, checked in by Martin Lee <martinlee84@…>, 11 years ago
chg: towards better EZGCD
  • Property mode set to 100644
File size: 20.6 KB
Line 
1/* emacs edit mode for this file is -*- C++ -*- */
2/* $Id$ */
3
4#include "config.h"
5
6#include "cf_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 "templates/ftmpl_functions.h"
18#include "cf_map.h"
19#include "facHensel.h"
20
21static
22int compress4EZGCD (const CanonicalForm& F, const CanonicalForm& G, CFMap & M,
23                    CFMap & N, int& both_non_zero)
24{
25  int n= tmax (F.level(), G.level());
26  int * degsf= new int [n + 1];
27  int * degsg= new int [n + 1];
28
29  for (int i = 0; i <= n; i++)
30    degsf[i]= degsg[i]= 0;
31
32  degsf= degrees (F, degsf);
33  degsg= degrees (G, degsg);
34
35  both_non_zero= 0;
36  int f_zero= 0;
37  int g_zero= 0;
38
39  for (int i= 1; i <= n; i++)
40  {
41    if (degsf[i] != 0 && degsg[i] != 0)
42    {
43      both_non_zero++;
44      continue;
45    }
46    if (degsf[i] == 0 && degsg[i] != 0 && i <= G.level())
47    {
48      f_zero++;
49      continue;
50    }
51    if (degsg[i] == 0 && degsf[i] && i <= F.level())
52    {
53      g_zero++;
54      continue;
55    }
56  }
57
58  if (both_non_zero == 0)
59  {
60    delete [] degsf;
61    delete [] degsg;
62    return 0;
63  }
64
65  // map Variables which do not occur in both polynomials to higher levels
66  int k= 1;
67  int l= 1;
68  for (int i= 1; i <= n; i++)
69  {
70    if (degsf[i] != 0 && degsg[i] == 0 && i <= F.level())
71    {
72      if (k + both_non_zero != i)
73      {
74        M.newpair (Variable (i), Variable (k + both_non_zero));
75        N.newpair (Variable (k + both_non_zero), Variable (i));
76      }
77      k++;
78    }
79    if (degsf[i] == 0 && degsg[i] != 0 && i <= G.level())
80    {
81      if (l + g_zero + both_non_zero != i)
82      {
83        M.newpair (Variable (i), Variable (l + g_zero + both_non_zero));
84        N.newpair (Variable (l + g_zero + both_non_zero), Variable (i));
85      }
86      l++;
87    }
88  }
89
90  // sort Variables x_{i} in decreasing order of
91  // min(deg_{x_{i}}(f),deg_{x_{i}}(g))
92  int m= tmin (F.level(), G.level());
93  int max_min_deg;
94  k= both_non_zero;
95  l= 0;
96  int i= 1;
97  while (k > 0)
98  {
99    max_min_deg= tmin (degsf[i], degsg[i]);
100    while (max_min_deg == 0)
101    {
102      i++;
103      max_min_deg= tmin (degsf[i], degsg[i]);
104    }
105    for (int j= i + 1; j <= m; j++)
106    {
107      if ((tmin (degsf[j],degsg[j]) < max_min_deg) &&
108          (tmin (degsf[j], degsg[j]) != 0))
109      {
110        max_min_deg= tmin (degsf[j], degsg[j]);
111        l= j;
112      }
113    }
114
115    if (l != 0)
116    {
117      if (l != k)
118      {
119        M.newpair (Variable (l), Variable(k));
120        N.newpair (Variable (k), Variable(l));
121        degsf[l]= 0;
122        degsg[l]= 0;
123        l= 0;
124      }
125      else
126      {
127        degsf[l]= 0;
128        degsg[l]= 0;
129        l= 0;
130      }
131    }
132    else if (l == 0)
133    {
134      if (i != k)
135      {
136        M.newpair (Variable (i), Variable (k));
137        N.newpair (Variable (k), Variable (i));
138        degsf[i]= 0;
139        degsg[i]= 0;
140      }
141      else
142      {
143        degsf[i]= 0;
144        degsg[i]= 0;
145      }
146      i++;
147    }
148    k--;
149  }
150
151  delete [] degsf;
152  delete [] degsg;
153
154  return both_non_zero;
155}
156
157static inline
158CanonicalForm myShift2Zero (const CanonicalForm& F, CFList& Feval,
159                            const CFList& evaluation)
160{
161  CanonicalForm A= F;
162  int k= 2;
163  for (CFListIterator i= evaluation; i.hasItem(); i++, k++)
164    A= A (Variable (k) + i.getItem(), k);
165
166  CanonicalForm buf= A;
167  Feval= CFList();
168  Feval.append (buf);
169  for (k= evaluation.length() + 1; k > 2; k--)
170  {
171    buf= mod (buf, Variable (k));
172    Feval.insert (buf);
173  }
174  return A;
175}
176
177static inline
178CanonicalForm myReverseShift (const CanonicalForm& F, const CFList& evaluation)
179{
180  int l= evaluation.length() + 1;
181  CanonicalForm result= F;
182  CFListIterator j= evaluation;
183  for (int i= 2; i < l + 1; i++, j++)
184  {
185    if (F.level() < i)
186      continue;
187    result= result (Variable (i) - j.getItem(), i);
188  }
189  return result;
190}
191
192static inline
193Evaluation optimize4Lift (const CanonicalForm& F, CFMap & M,
194                    CFMap & N, const Evaluation& A)
195{
196  int n= F.level();
197  int * degsf= new int [n + 1];
198
199  for (int i = 0; i <= n; i++)
200    degsf[i]= 0;
201
202  degsf= degrees (F, degsf);
203
204  Evaluation result= Evaluation (A.min(), A.max());
205  ASSERT (A.min() == 2, "expected A.min() == 2");
206  ASSERT (A.max() == n, "expected A.max() == n");
207  int max_deg;
208  int k= n;
209  int l= 1;
210  int i= 2;
211  int pos= 2;
212  while (k > 1)
213  {
214    max_deg= degsf [i];
215    while (max_deg == 0)
216    {
217      i++;
218      max_deg= degsf [i];
219    }
220    l= i;
221    for (int j= i + 1; j <=  n; j++)
222    {
223      if (degsf[j] > max_deg)
224      {
225        max_deg= degsf[j];
226        l= j;
227      }
228    }
229
230    if (l <= n)
231    {
232      if (l != pos)
233      {
234        result.setValue (pos, A [l]);
235        M.newpair (Variable (l), Variable (pos));
236        N.newpair (Variable (pos), Variable (l));
237        degsf[l]= 0;
238        l= 2;
239        if (k == 2 && n == 3)
240        {
241          result.setValue (l, A [pos]);
242          M.newpair (Variable (pos), Variable (l));
243          N.newpair (Variable (l), Variable (pos));
244          degsf[pos]= 0;
245        }
246      }
247      else
248      {
249        result.setValue (l, A [l]);
250        degsf [l]= 0;
251      }
252    }
253    pos++;
254    k--;
255    l= 2;
256  }
257
258  delete [] degsf;
259
260  return result;
261}
262
263static inline
264int Hensel (const CanonicalForm & UU, CFArray & G, const Evaluation & AA,
265            const CFArray& LeadCoeffs )
266{
267  CFList factors;
268  factors.append (G[1]);
269  factors.append (G[2]);
270
271  CFMap NN, MM;
272  Evaluation A= optimize4Lift (UU, MM, NN, AA);
273
274  CanonicalForm U= MM (UU);
275  CFArray LCs= CFArray (1,2);
276  LCs [1]= MM (LeadCoeffs [1]);
277  LCs [2]= MM (LeadCoeffs [2]);
278
279  CFList evaluation;
280  for (int i= A.min(); i <= A.max(); i++)
281    evaluation.append (A [i]);
282  CFList UEval;
283  CanonicalForm shiftedU= myShift2Zero (U, UEval, evaluation);
284
285  if (size (shiftedU)/getNumVars (U) > 500)
286    return -1;
287
288  CFArray shiftedLCs= CFArray (2);
289  CFList shiftedLCsEval1, shiftedLCsEval2;
290  shiftedLCs[0]= myShift2Zero (LCs[1], shiftedLCsEval1, evaluation);
291  shiftedLCs[1]= myShift2Zero (LCs[2], shiftedLCsEval2, evaluation);
292  factors.insert (1);
293  int liftBound= degree (UEval.getLast(), 2) + 1;
294  CFArray Pi;
295  CFMatrix M= CFMatrix (liftBound, factors.length() - 1);
296  CFList diophant;
297  CFArray lcs= CFArray (2);
298  lcs [0]= shiftedLCsEval1.getFirst();
299  lcs [1]= shiftedLCsEval2.getFirst();
300  nonMonicHenselLift12 (UEval.getFirst(), factors, liftBound, Pi, diophant, M,
301                        lcs, false);
302
303  for (CFListIterator i= factors; i.hasItem(); i++)
304  {
305    if (!fdivides (i.getItem(), UEval.getFirst()))
306      return 0;
307  }
308
309  int * liftBounds;
310  bool noOneToOne= false;
311  if (U.level() > 2)
312  {
313    liftBounds= new int [U.level() - 1]; /* index: 0.. U.level()-2 */
314    liftBounds[0]= liftBound;
315    for (int i= 1; i < U.level() - 1; i++)
316      liftBounds[i]= degree (shiftedU, Variable (i + 2)) + 1;
317    factors= nonMonicHenselLift2 (UEval, factors, liftBounds, U.level() - 1,
318                                  false, shiftedLCsEval1, shiftedLCsEval2, Pi,
319                                  diophant, noOneToOne);
320    delete [] liftBounds;
321    if (noOneToOne)
322      return 0;
323  }
324  G[1]= factors.getFirst();
325  G[2]= factors.getLast();
326  G[1]= myReverseShift (G[1], evaluation);
327  G[2]= myReverseShift (G[2], evaluation);
328  G[1]= NN (G[1]);
329  G[2]= NN (G[2]);
330  return 1;
331}
332
333static void findeval( const CanonicalForm & F, const CanonicalForm & G, CanonicalForm & Fb, CanonicalForm & Gb, CanonicalForm & Db, REvaluation & b, int delta, int degF, int degG );
334
335static CanonicalForm ezgcd ( const CanonicalForm & FF, const CanonicalForm & GG, REvaluation & b, bool internal );
336
337static CanonicalForm ezgcd_specialcase ( const CanonicalForm & F, const CanonicalForm & G, REvaluation & b, const modpk & bound );
338
339static modpk findBound ( const CanonicalForm & F, const CanonicalForm & G, const CanonicalForm & lcF, const CanonicalForm & lcG, int degF, int degG );
340
341static modpk enlargeBound ( const CanonicalForm & F, const CanonicalForm & Lb, const CanonicalForm & Db, const modpk & pk );
342
343CanonicalForm
344ezgcd ( const CanonicalForm & FF, const CanonicalForm & GG )
345{
346    REvaluation b;
347    return ezgcd( FF, GG, b, false );
348}
349
350static CanonicalForm
351ezgcd ( const CanonicalForm & FF, const CanonicalForm & GG, REvaluation & b, bool internal )
352{
353    CanonicalForm F, G, f, g, d, Fb, Gb, Db, Fbt, Gbt, Dbt, B0, B, D0, lcF, lcG, lcD, cand, result;
354    CFArray DD( 1, 2 ), lcDD( 1, 2 );
355    int degF, degG, delta, t;
356    REvaluation bt;
357    bool gcdfound = false;
358    Variable x = Variable(1);
359    modpk bound;
360
361    F= FF;
362    G= GG;
363
364    CFMap M,N;
365    int smallestDegLev;
366    int best_level= compress4EZGCD (F, G, M, N, smallestDegLev);
367
368    if (best_level == 0) return G.genOne();
369
370    F= M (F);
371    G= M (G);
372
373
374    DEBINCLEVEL( cerr, "ezgcd" );
375    DEBOUTLN( cerr, "FF = " << FF );
376    DEBOUTLN( cerr, "GG = " << GG );
377    f = content( F, x ); g = content( G, x ); d = gcd( f, g );
378    DEBOUTLN( cerr, "f = " << f );
379    DEBOUTLN( cerr, "g = " << g );
380    F /= f; G /= g;
381    if ( F.isUnivariate() && G.isUnivariate() )
382    {
383        DEBDECLEVEL( cerr, "ezgcd" );
384        if(F.mvar()==G.mvar())
385          d*=gcd(F,G);
386        return N (d);
387    }
388    else  if ( gcd_test_one( F, G, false ) )
389    {
390        DEBDECLEVEL( cerr, "ezgcd" );
391        return N (d);
392    }
393    lcF = LC( F, x ); lcG = LC( G, x );
394    lcD = gcd( lcF, lcG );
395    delta = 0;
396    degF = degree( F, x ); degG = degree( G, x );
397    t = tmax( F.level(), G.level() );
398    //bound = findBound( F, G, lcF, lcG, degF, degG );
399    if ( ! internal )
400        b = REvaluation( 2, t, IntRandom( 25 ) );
401    while ( ! gcdfound )
402    {
403        /// ---> A2
404        DEBOUTLN( cerr, "search for evaluation, delta = " << delta );
405        DEBOUTLN( cerr, "F = " << F );
406        DEBOUTLN( cerr, "G = " << G );
407        findeval( F, G, Fb, Gb, Db, b, delta, degF, degG );
408        DEBOUTLN( cerr, "found evaluation b = " << b );
409        DEBOUTLN( cerr, "F(b) = " << Fb );
410        DEBOUTLN( cerr, "G(b) = " << Gb );
411        DEBOUTLN( cerr, "D(b) = " << Db );
412        delta = degree( Db );
413        /// ---> A3
414        if ( delta == 0 )
415        {
416          DEBDECLEVEL( cerr, "ezgcd" );
417          return N (d);
418        }
419        /// ---> A4
420        //deltaold = delta;
421        while ( 1 ) {
422            bt = b;
423            findeval( F, G, Fbt, Gbt, Dbt, bt, delta + 1, degF, degG );
424            int dd=degree( Dbt );
425            if ( dd /*degree( Dbt )*/ == 0 )
426            {
427                DEBDECLEVEL( cerr, "ezgcd" );
428                return N (d);
429            }
430            if ( dd /*degree( Dbt )*/ == delta )
431                break;
432            else  if ( dd /*degree( Dbt )*/ < delta ) {
433                delta = dd /*degree( Dbt )*/;
434                b = bt;
435                Db = Dbt; Fb = Fbt; Gb = Gbt;
436            }
437        }
438        DEBOUTLN( cerr, "now after A4, delta = " << delta );
439        /// ---> A5
440        if (delta == degF)
441        {
442          if (degF <= degG  && fdivides (F, G))
443          {
444            DEBDECLEVEL( cerr, "ezgcd" );
445            return N (d*F);
446          }
447          else
448            delta--;
449        }
450        else if (delta == degG)
451        {
452          if (degG <= degF && fdivides( G, F ))
453          {
454            DEBDECLEVEL( cerr, "ezgcd" );
455            return N (d*G);
456          }
457          else
458            delta--;
459        }
460        if ( delta != degF && delta != degG )
461        {
462            /// ---> A6
463            bool B_is_F;
464            CanonicalForm xxx1, xxx2;
465            CanonicalForm buf;
466            DD[1] = Fb / Db;
467            buf= Gb/Db;
468            xxx1 = gcd( DD[1], Db );
469            xxx2 = gcd( buf, Db );
470            if (((xxx1.inCoeffDomain() && xxx2.inCoeffDomain()) &&
471                (size (F) <= size (G)))
472                 || (xxx1.inCoeffDomain() && !xxx2.inCoeffDomain()))
473            {
474              B = F;
475              DD[2] = Db;
476              lcDD[1] = lcF;
477              lcDD[2] = lcD;
478              B_is_F = true;
479            }
480            else if (((xxx1.inCoeffDomain() && xxx2.inCoeffDomain()) &&
481                     (size (G) < size (F)))
482                     || (!xxx1.inCoeffDomain() && xxx2.inCoeffDomain()))
483            {
484              DD[1] = buf;
485              B = G;
486              DD[2] = Db;
487              lcDD[1] = lcG;
488              lcDD[2] = lcD;
489              B_is_F = false;
490            }
491            else
492            {
493              Off(SW_USE_EZGCD);
494              result=gcd( F, G );
495              On(SW_USE_EZGCD);
496              return N (d*result);
497            }
498            /*CanonicalForm xxx;
499            //if ( gcd( (DD[1] = Fb / Db), Db ) == 1 ) {
500            DD[1] = Fb / Db;
501            xxx= gcd( DD[1], Db );
502            DEBOUTLN( cerr, "gcd((Fb/Db),Db) = " << xxx );
503            DEBOUTLN( cerr, "Fb/Db = " << DD[1] );
504            DEBOUTLN( cerr, "Db = " << Db );
505            if (xxx.inCoeffDomain()) {
506                B = F;
507                DD[2] = Db;
508                lcDD[1] = lcF;
509                lcDD[2] = lcF;
510                B *= lcF;
511                B_is_F=true;
512            }
513            //else  if ( gcd( (DD[1] = Gb / Db), Db ) == 1 ) {
514            else
515            {
516              DD[1] = Gb / Db;
517              xxx=gcd( DD[1], Db );
518              DEBOUTLN( cerr, "gcd((Gb/Db),Db) = " << xxx );
519              DEBOUTLN( cerr, "Gb/Db = " << DD[1] );
520              DEBOUTLN( cerr, "Db = " << Db );
521              if (xxx.inCoeffDomain())
522              {
523                B = G;
524                DD[2] = Db;
525                lcDD[1] = lcG;
526                lcDD[2] = lcG;
527                B *= lcG;
528                B_is_F=false;
529              }
530              else
531              {
532#ifdef DEBUGOUTPUT
533                CanonicalForm dummyres = d * ezgcd_specialcase( F, G, b, bound );
534                DEBDECLEVEL( cerr, "ezgcd" );
535                return dummyres;
536#else
537                return d * ezgcd_specialcase( F, G, b, bound );
538#endif
539              }
540            }*/
541            /// ---> A7
542            DD[2] = DD[2] * ( b( lcDD[2] ) / lc( DD[2] ) );
543            DD[1] = DD[1] * ( b( lcDD[1] ) / lc( DD[1] ) );
544            //bound = enlargeBound( B, DD[2], DD[1], bound );
545            DEBOUTLN( cerr, "(hensel) B    = " << B );
546            DEBOUTLN( cerr, "(hensel) lcB  = " << LC( B, Variable(1) ) );
547            DEBOUTLN( cerr, "(hensel) b(B) = " << b(B) );
548            DEBOUTLN( cerr, "(hensel) DD   = " << DD );
549            DEBOUTLN( cerr, "(hensel) lcDD = " << lcDD );
550            gcdfound= Hensel (B*lcD, DD, b, lcDD);
551            DEBOUTLN( cerr, "(hensel finished) DD   = " << DD );
552
553            if (gcdfound == -1)
554            {
555              Off (SW_USE_EZGCD);
556              result= gcd (F,G);
557              On (SW_USE_EZGCD_P);
558              return N (d*result);
559            }
560
561            if (gcdfound)
562            {
563              cand = DD[2] / content( DD[2], Variable(1) );
564              gcdfound = fdivides( cand, G ) && fdivides ( cand, F );
565            }
566            /// ---> A8 (gcdfound)
567        }
568        delta--;
569    }
570    /// ---> A9
571    DEBDECLEVEL( cerr, "ezgcd" );
572    return N (d*cand);
573}
574
575static CanonicalForm
576ezgcd_specialcase ( const CanonicalForm & F, const CanonicalForm & G, REvaluation & /*b*/, const modpk & /*bound*/ )
577{
578    CanonicalForm d;
579#if 1
580    Off(SW_USE_EZGCD);
581    //bool ntl0=isOn(SW_USE_NTL_GCD_0);
582    //Off(SW_USE_NTL_GCD_0);
583    //bool ntlp=isOn(SW_USE_NTL_GCD_P);
584    //Off(SW_USE_NTL_GCD_P);
585    d=gcd( F, G );
586    //if (ntl0) On(SW_USE_NTL_GCD_0);
587    //if (ntlp) On(SW_USE_NTL_GCD_P);
588    On(SW_USE_EZGCD);
589    return d;
590#else
591    /*
592     * I am not sure, if these lines are much of use.
593     * Remove?
594     */
595    DEBOUTLN( cerr, "ezgcd: special case" );
596    CanonicalForm Ft, Gt, L, LL, Fb, Gb, Db, Lb, D, Ds, Dt;
597    CFArray DD( 1, 2 ), lcDD( 1, 2 );
598    Variable x = Variable( 1 );
599    bool gcdfound;
600
601    Dt = 1;
602    /// ---> S1
603    DEBOUTLN( cerr, "ezgcdspec: (S1)" );
604    Ft = ezgcd( F, F.deriv( x ) );
605    if ( Ft.isOne() )
606    {
607        // In this case F is squarefree and we came here by bad chance
608        // (means: bad evaluation point).  Try again with another
609        // evaluation point.  Bug fix (?) by JS.  The bad example was
610        // gcd.debug -ocr /+USE_EZGCD/@12/CB
611        //     '(16*B^8-208*B^6*C+927*B^4*C^2-1512*B^2*C^3+432*C^4)'
612        //     '(4*B^7*C^2-50*B^5*C^3+208*B^3*C^4-288*B*C^5)'
613        b.nextpoint();
614        return ezgcd( F, G, b, true );
615    }
616
617    DEBOUTLN( cerr, "ezgcdspec: (S1) done, Ft = " << Ft );
618    L = F / Ft;
619    /// ---> S2
620    DEBOUTLN( cerr, "ezgcdspec: (S2)" );
621
622    L = ezgcd( L, G, b, true );
623    DEBOUTLN( cerr, "ezgcdspec: (S2) done, Ds = " << Ds );
624    Gt = G / L;
625    Ds = L; Dt = L;
626    Lb = b( L ); Gb = b( Gt ); Fb = b( Ft );
627    d = gcd( LC( L, x ), gcd( LC( Ft, x ), LC( Gt, x ) ) );
628    while ( true ) {
629        /// ---> S3
630        DEBOUTLN( cerr, "ezgcdspec: (S3)" );
631        DEBOUTLN( cerr, "ezgcdspec: Fb = " << Fb );
632        DEBOUTLN( cerr, "ezgcdspec: Gb = " << Gb );
633        DD[1] = gcd( Lb, gcd( Fb, Gb ) );
634        /// ---> S4
635        DEBOUTLN( cerr, "ezgcdspec: (S4)" );
636        if ( degree( DD[1] ) == 0 )
637            return Ds;
638        /// ---> S5
639        DEBOUTLN( cerr, "ezgcdspec: (S5)" );
640        DEBOUTLN( cerr, "ezgcdspec: Lb = " << Lb );
641        DEBOUTLN( cerr, "ezgcdspec: Db = " << DD[1] );
642        Db = DD[1];
643        if ( ! (DD[2] = Lb/DD[1]).isOne() ) {
644            DEBOUTLN( cerr, "ezgcdspec: (S55)" );
645            lcDD[2] = LC( L, x );
646            lcDD[1] = d;
647            DD[1] = ( DD[1] * b( d ) ) / lc( DD[1] );
648            DD[2] = ( DD[2] * b( lcDD[2] ) ) / lc( DD[2] );
649            LL = L * d;
650            modpk newbound = enlargeBound( LL, DD[2], DD[1], bound );
651            DEBOUTLN( cerr, "ezgcdspec: begin with lift." );
652            DEBOUTLN( cerr, "ezgcdspec: B     = " << LL );
653            DEBOUTLN( cerr, "ezgcdspec: DD    = " << DD );
654            DEBOUTLN( cerr, "ezgcdspec: lcDD  = " << lcDD );
655            DEBOUTLN( cerr, "ezgcdspec: b     = " << b );
656            DEBOUTLN( cerr, "ezgcdspec: bound = " << bound.getpk() );
657            DEBOUTLN( cerr, "ezgcdspec: lc(B) = " << LC( LL, x ) );
658            DEBOUTLN( cerr, "ezgcdspec: test  = " << b( LL ) - DD[1] * DD[2] );
659            gcdfound = Hensel( LL, DD, lcDD, b, newbound, x );
660            ASSERT( gcdfound, "fatal error in algorithm" );
661            Dt = pp( DD[1] );
662        }
663        DEBOUTLN( cerr, "ezgcdspec: (S7)" );
664        Ds = Ds * Dt;
665        Fb = Fb / Db;
666        Gb = Gb / Db;
667    }
668    // this point is never reached, but to avoid compiler warnings let's return a value
669    return 0;
670#endif
671}
672
673static void
674findeval( const CanonicalForm & F, const CanonicalForm & G, CanonicalForm & Fb, CanonicalForm & Gb, CanonicalForm & Db, REvaluation & b, int delta, int degF, int degG )
675{
676    bool ok;
677    if ( delta != 0 )
678        b.nextpoint();
679    DEBOUTLN( cerr, "ezgcd: (findeval) F = " << F  <<", G="<< G);
680    DEBOUTLN( cerr, "ezgcd: (findeval) degF = " << degF << ", degG="<<degG );
681    do {
682        DEBOUTLN( cerr, "ezgcd: (findeval) b = " << b );
683        Fb = b( F );
684        ok = degree( Fb ) == degF;
685        if ( ok ) {
686            Gb = b( G );
687            ok = degree( Gb ) == degG;
688        }
689
690        if ( ok )
691        {
692            Db = gcd( Fb, Gb );
693            if ( delta > 0 )
694              ok = degree( Db ) < delta;
695        }
696        if ( ! ok )
697        {
698            b.nextpoint();
699        }
700    } while ( ! ok );
701}
702
703static modpk
704enlargeBound ( const CanonicalForm & F, const CanonicalForm & Lb, const CanonicalForm & Db, const modpk & pk )
705{
706    DEBOUTLN( cerr, "ezgcd: (enlarge bound) F      = " << F );
707    DEBOUTLN( cerr, "ezgcd: (enlarge bound) Lb     = " << Lb );
708    DEBOUTLN( cerr, "ezgcd: (enlarge bound) Db     = " << Db );
709    DEBOUTLN( cerr, "ezgcd: (enlarge bound) Lb % p = " << mod( Lb, pk.getp() ) );
710    DEBOUTLN( cerr, "ezgcd: (enlarge bound) Db % p = " << mod( Db, pk.getp() ) );
711
712    CanonicalForm limit = power( CanonicalForm(2), degree( Db ) ) *
713        tmax( maxNorm( Lb ), tmax( maxNorm( Db ), maxNorm( F ) ) );
714    int p = pk.getp();
715    int k = pk.getk();
716    CanonicalForm bound = pk.getpk();
717    while ( bound < limit ) {
718        k++;
719        bound *= p;
720    }
721    k *= 2;
722    DEBOUTLN( cerr, "ezgcd: (enlarge bound) newbound = " << power( CanonicalForm( p ), k ) );
723    return modpk( p, k );
724}
725
726static modpk
727findBound ( const CanonicalForm & F, const CanonicalForm & G, const CanonicalForm & lcF, const CanonicalForm & lcG, int degF, int degG )
728{
729    CanonicalForm limit = power( CanonicalForm(2), tmin( degF, degG ) ) *
730        gcd( icontent( lcF ), icontent( lcG ) ) * tmin( maxNorm( F ), maxNorm( G ) );
731    int p, i = 0;
732    do {
733        p = cf_getBigPrime( i );
734        i++;
735    } while ( mod( lcF, p ).isZero() && mod( lcG, p ).isZero() );
736    CanonicalForm bound = p;
737    i = 1;
738    while ( bound < limit ) {
739        i++;
740        bound *= p;
741    }
742    return modpk( p, i );
743}
Note: See TracBrowser for help on using the repository browser.