source: git/factory/cfEzgcd.cc @ 52a933f

fieker-DuValspielwiese
Last change on this file since 52a933f was 52a933f, checked in by Martin Lee <martinlee84@…>, 10 years ago
chg: renamed modular GCD variants in cf_gcd_smallp to modGCDFp, modGCDFq, modGCDGF, modGCDZ
  • Property mode set to 100644
File size: 36.8 KB
Line 
1/*****************************************************************************\
2 * Computer Algebra System SINGULAR
3\*****************************************************************************/
4/** @file cfEzgcd.cc
5 *
6 * This file implements the GCD of two multivariate polynomials over Q or F_q
7 * using EZ-GCD as described in "Algorithms for Computer Algebra" by Geddes,
8 * Czapor, Labahnn
9 *
10 * @author Martin Lee
11 *
12 **/
13/*****************************************************************************/
14
15
16#include "config.h"
17
18
19#include "timing.h"
20#include "cf_assert.h"
21#include "debug.h"
22
23#include "cf_defs.h"
24#include "canonicalform.h"
25#include "cfEzgcd.h"
26#include "cf_gcd_smallp.h"
27#include "cf_util.h"
28#include "cf_map_ext.h"
29#include "cf_algorithm.h"
30#include "cf_reval.h"
31#include "cf_random.h"
32#include "cf_primes.h"
33#include "templates/ftmpl_functions.h"
34#include "cf_map.h"
35#include "facHensel.h"
36
37#ifdef HAVE_NTL
38#include "NTLconvert.h"
39
40static const double log2exp= 1.442695041;
41
42TIMING_DEFINE_PRINT(ez_eval)
43TIMING_DEFINE_PRINT(ez_compress)
44TIMING_DEFINE_PRINT(ez_hensel_lift)
45TIMING_DEFINE_PRINT(ez_content)
46TIMING_DEFINE_PRINT(ez_termination)
47
48static
49int compress4EZGCD (const CanonicalForm& F, const CanonicalForm& G, CFMap & M,
50                    CFMap & N, int& both_non_zero)
51{
52  int n= tmax (F.level(), G.level());
53  int * degsf= new int [n + 1];
54  int * degsg= new int [n + 1];
55
56  for (int i = 0; i <= n; i++)
57    degsf[i]= degsg[i]= 0;
58
59  degsf= degrees (F, degsf);
60  degsg= degrees (G, degsg);
61
62  both_non_zero= 0;
63  int f_zero= 0;
64  int g_zero= 0;
65
66  for (int i= 1; i <= n; i++)
67  {
68    if (degsf[i] != 0 && degsg[i] != 0)
69    {
70      both_non_zero++;
71      continue;
72    }
73    if (degsf[i] == 0 && degsg[i] != 0 && i <= G.level())
74    {
75      f_zero++;
76      continue;
77    }
78    if (degsg[i] == 0 && degsf[i] && i <= F.level())
79    {
80      g_zero++;
81      continue;
82    }
83  }
84
85  if (both_non_zero == 0)
86  {
87    delete [] degsf;
88    delete [] degsg;
89    return 0;
90  }
91
92  // map Variables which do not occur in both polynomials to higher levels
93  int k= 1;
94  int l= 1;
95  for (int i= 1; i <= n; i++)
96  {
97    if (degsf[i] != 0 && degsg[i] == 0 && i <= F.level())
98    {
99      if (k + both_non_zero != i)
100      {
101        M.newpair (Variable (i), Variable (k + both_non_zero));
102        N.newpair (Variable (k + both_non_zero), Variable (i));
103      }
104      k++;
105    }
106    if (degsf[i] == 0 && degsg[i] != 0 && i <= G.level())
107    {
108      if (l + g_zero + both_non_zero != i)
109      {
110        M.newpair (Variable (i), Variable (l + g_zero + both_non_zero));
111        N.newpair (Variable (l + g_zero + both_non_zero), Variable (i));
112      }
113      l++;
114    }
115  }
116
117  // sort Variables x_{i} in decreasing order of
118  // min(deg_{x_{i}}(f),deg_{x_{i}}(g))
119  int m= tmin (F.level(), G.level());
120  int max_min_deg;
121  k= both_non_zero;
122  l= 0;
123  int i= 1;
124  while (k > 0)
125  {
126    max_min_deg= tmin (degsf[i], degsg[i]);
127    while (max_min_deg == 0)
128    {
129      i++;
130      max_min_deg= tmin (degsf[i], degsg[i]);
131    }
132    for (int j= i + 1; j <= m; j++)
133    {
134      if ((tmin (degsf[j],degsg[j]) < max_min_deg) &&
135          (tmin (degsf[j], degsg[j]) != 0))
136      {
137        max_min_deg= tmin (degsf[j], degsg[j]);
138        l= j;
139      }
140    }
141
142    if (l != 0)
143    {
144      if (l != k)
145      {
146        M.newpair (Variable (l), Variable(k));
147        N.newpair (Variable (k), Variable(l));
148        degsf[l]= 0;
149        degsg[l]= 0;
150        l= 0;
151      }
152      else
153      {
154        degsf[l]= 0;
155        degsg[l]= 0;
156        l= 0;
157      }
158    }
159    else if (l == 0)
160    {
161      if (i != k)
162      {
163        M.newpair (Variable (i), Variable (k));
164        N.newpair (Variable (k), Variable (i));
165        degsf[i]= 0;
166        degsg[i]= 0;
167      }
168      else
169      {
170        degsf[i]= 0;
171        degsg[i]= 0;
172      }
173      i++;
174    }
175    k--;
176  }
177
178  delete [] degsf;
179  delete [] degsg;
180
181  return both_non_zero;
182}
183
184static inline
185CanonicalForm myShift2Zero (const CanonicalForm& F, CFList& Feval,
186                            const CFList& evaluation)
187{
188  CanonicalForm A= F;
189  int k= 2;
190  for (CFListIterator i= evaluation; i.hasItem(); i++, k++)
191    A= A (Variable (k) + i.getItem(), k);
192
193  CanonicalForm buf= A;
194  Feval= CFList();
195  Feval.append (buf);
196  for (k= evaluation.length() + 1; k > 2; k--)
197  {
198    buf= mod (buf, Variable (k));
199    Feval.insert (buf);
200  }
201  return A;
202}
203
204static inline
205CanonicalForm myReverseShift (const CanonicalForm& F, const CFList& evaluation)
206{
207  int l= evaluation.length() + 1;
208  CanonicalForm result= F;
209  CFListIterator j= evaluation;
210  for (int i= 2; i < l + 1; i++, j++)
211  {
212    if (F.level() < i)
213      continue;
214    result= result (Variable (i) - j.getItem(), i);
215  }
216  return result;
217}
218
219static inline
220Evaluation optimize4Lift (const CanonicalForm& F, CFMap & M,
221                    CFMap & N, const Evaluation& A)
222{
223  int n= F.level();
224  int * degsf= new int [n + 1];
225
226  for (int i = 0; i <= n; i++)
227    degsf[i]= 0;
228
229  degsf= degrees (F, degsf);
230
231  Evaluation result= Evaluation (A.min(), A.max());
232  ASSERT (A.min() == 2, "expected A.min() == 2");
233  int max_deg;
234  int k= n;
235  int l= 1;
236  int i= 2;
237  int pos= 2;
238  while (k > 1)
239  {
240    max_deg= degsf [i];
241    while (max_deg == 0)
242    {
243      i++;
244      max_deg= degsf [i];
245    }
246    l= i;
247    for (int j= i + 1; j <=  n; j++)
248    {
249      if (degsf[j] > max_deg)
250      {
251        max_deg= degsf[j];
252        l= j;
253      }
254    }
255
256    if (l <= n)
257    {
258      if (l != pos)
259      {
260        result.setValue (pos, A [l]);
261        M.newpair (Variable (l), Variable (pos));
262        N.newpair (Variable (pos), Variable (l));
263        degsf[l]= 0;
264        l= 2;
265        if (k == 2 && n == 3)
266        {
267          result.setValue (l, A [pos]);
268          M.newpair (Variable (pos), Variable (l));
269          N.newpair (Variable (l), Variable (pos));
270          degsf[pos]= 0;
271        }
272      }
273      else
274      {
275        result.setValue (l, A [l]);
276        degsf [l]= 0;
277      }
278    }
279    pos++;
280    k--;
281    l= 2;
282  }
283
284  delete [] degsf;
285
286  return result;
287}
288
289static inline
290int Hensel (const CanonicalForm & UU, CFArray & G, const Evaluation & AA,
291            const CFArray& LeadCoeffs )
292{
293  CFList factors;
294  factors.append (G[1]);
295  factors.append (G[2]);
296
297  CFMap NN, MM;
298  Evaluation A= optimize4Lift (UU, MM, NN, AA);
299
300  CanonicalForm U= MM (UU);
301  CFArray LCs= CFArray (1,2);
302  LCs [1]= MM (LeadCoeffs [1]);
303  LCs [2]= MM (LeadCoeffs [2]);
304
305  CFList evaluation;
306  long termEstimate= size (U);
307  for (int i= A.min(); i <= A.max(); i++)
308  {
309    if (!A[i].isZero())
310    {
311      termEstimate *= degree (U,i)*2;
312      termEstimate /= 3;
313    }
314    evaluation.append (A [i]);
315  }
316  if (termEstimate/getNumVars(U) > 500)
317    return -1;
318  CFList UEval;
319  CanonicalForm shiftedU= myShift2Zero (U, UEval, evaluation);
320
321  if (size (shiftedU)/getNumVars (U) > 500)
322    return -1;
323
324  CFArray shiftedLCs= CFArray (2);
325  CFList shiftedLCsEval1, shiftedLCsEval2;
326  shiftedLCs[0]= myShift2Zero (LCs[1], shiftedLCsEval1, evaluation);
327  shiftedLCs[1]= myShift2Zero (LCs[2], shiftedLCsEval2, evaluation);
328  factors.insert (1);
329  int liftBound= degree (UEval.getLast(), 2) + 1;
330  CFArray Pi;
331  CFMatrix M= CFMatrix (liftBound, factors.length() - 1);
332  CFList diophant;
333  CFArray lcs= CFArray (2);
334  lcs [0]= shiftedLCsEval1.getFirst();
335  lcs [1]= shiftedLCsEval2.getFirst();
336  nonMonicHenselLift12 (UEval.getFirst(), factors, liftBound, Pi, diophant, M,
337                        lcs, false);
338
339  for (CFListIterator i= factors; i.hasItem(); i++)
340  {
341    if (!fdivides (i.getItem(), UEval.getFirst()))
342      return 0;
343  }
344
345  int * liftBounds;
346  bool noOneToOne= false;
347  if (U.level() > 2)
348  {
349    liftBounds= new int [U.level() - 1]; /* index: 0.. U.level()-2 */
350    liftBounds[0]= liftBound;
351    for (int i= 1; i < U.level() - 1; i++)
352      liftBounds[i]= degree (shiftedU, Variable (i + 2)) + 1;
353    factors= nonMonicHenselLift2 (UEval, factors, liftBounds, U.level() - 1,
354                                  false, shiftedLCsEval1, shiftedLCsEval2, Pi,
355                                  diophant, noOneToOne);
356    delete [] liftBounds;
357    if (noOneToOne)
358      return 0;
359  }
360  G[1]= factors.getFirst();
361  G[2]= factors.getLast();
362  G[1]= myReverseShift (G[1], evaluation);
363  G[2]= myReverseShift (G[2], evaluation);
364  G[1]= NN (G[1]);
365  G[2]= NN (G[2]);
366  return 1;
367}
368
369static
370bool findeval (const CanonicalForm & F, const CanonicalForm & G,
371               CanonicalForm & Fb, CanonicalForm & Gb, CanonicalForm & Db,
372               REvaluation & b, int delta, int degF, int degG, int maxeval,
373               int & count, int& k, int bound, int& l)
374{
375  if( count == 0 && delta != 0)
376  {
377    if( count++ > maxeval )
378      return false;
379  }
380  if (count > 0)
381  {
382    b.nextpoint(k);
383    if (k == 0)
384      k++;
385    l++;
386    if (l > bound)
387    {
388      l= 1;
389      k++;
390      if (k > tmax (F.level(), G.level()) - 1)
391        return false;
392      b.nextpoint (k);
393    }
394    if (count++ > maxeval)
395      return false;
396  }
397  while( true )
398  {
399    Fb = b( F );
400    if( degree( Fb, 1 ) == degF )
401    {
402      Gb = b( G );
403      if( degree( Gb, 1 ) == degG )
404      {
405        Db = gcd( Fb, Gb );
406        if( delta > 0 )
407        {
408          if( degree( Db, 1 ) <= delta )
409            return true;
410        }
411        else
412        {
413          k++;
414          return true;
415        }
416      }
417    }
418    if (k == 0)
419      k++;
420    b.nextpoint(k);
421    l++;
422    if (l > bound)
423    {
424      l= 1;
425      k++;
426      if (k > tmax (F.level(), G.level()) - 1)
427        return false;
428      b.nextpoint (k);
429    }
430    if( count++ > maxeval )
431      return false;
432  }
433}
434
435static CanonicalForm
436ezgcd ( const CanonicalForm & FF, const CanonicalForm & GG, REvaluation & b,
437        bool internal )
438{
439  bool isRat= isOn (SW_RATIONAL);
440
441  int maxNumVars= tmax (getNumVars (FF), getNumVars (GG));
442  int sizeF= size (FF);
443  int sizeG= size (GG);
444
445
446  if (!isRat)
447    On (SW_RATIONAL);
448  if (sizeF/maxNumVars > 500 && sizeG/maxNumVars > 500)
449  {
450    Off(SW_USE_EZGCD);
451    CanonicalForm result=gcd( FF, GG );
452    On(SW_USE_EZGCD);
453    if (!isRat)
454      Off (SW_RATIONAL);
455    result /= icontent (result);
456    DEBDECLEVEL( cerr, "ezgcd" );
457    return result;
458  }
459
460
461  CanonicalForm F, G, f, g, d, Fb, Gb, Db, Fbt, Gbt, Dbt, B0, B, D0, lcF, lcG,
462                lcD, cand, contcand, result;
463  CFArray DD( 1, 2 ), lcDD( 1, 2 );
464  int degF, degG, delta, t, count, maxeval;
465  REvaluation bt;
466  int gcdfound = 0;
467  Variable x = Variable(1);
468  count= 0;
469  maxeval= 200;
470  int o, l;
471  o= 0;
472  l= 1;
473
474  if (!isRat)
475    On (SW_RATIONAL);
476  F= FF*bCommonDen (FF);
477  G= GG*bCommonDen (GG);
478  if (!isRat)
479    Off (SW_RATIONAL);
480
481  TIMING_START (ez_compress)
482  CFMap M,N;
483  int smallestDegLev;
484  int best_level= compress4EZGCD (F, G, M, N, smallestDegLev);
485
486  if (best_level == 0)
487  {
488    DEBDECLEVEL( cerr, "ezgcd" );
489    return G.genOne();
490  }
491
492  F= M (F);
493  G= M (G);
494  TIMING_END_AND_PRINT (ez_compress, "time for compression in EZ: ")
495
496  DEBINCLEVEL( cerr, "ezgcd" );
497  DEBOUTLN( cerr, "FF = " << FF );
498  DEBOUTLN( cerr, "GG = " << GG );
499  TIMING_START (ez_content)
500  f = content( F, x ); g = content( G, x ); d = gcd( f, g );
501  DEBOUTLN( cerr, "f = " << f );
502  DEBOUTLN( cerr, "g = " << g );
503  F /= f; G /= g;
504  TIMING_END_AND_PRINT (ez_content, "time to extract content in EZ: ")
505  if ( F.isUnivariate() && G.isUnivariate() )
506  {
507    DEBDECLEVEL( cerr, "ezgcd" );
508    if(F.mvar()==G.mvar())
509      d*=gcd(F,G);
510    else
511      return N (d);
512    return N (d);
513  }
514  if ( F.isUnivariate())
515  {
516    g= content (G,G.mvar());
517    return N(d*gcd(F,g));
518  }
519  if ( G.isUnivariate())
520  {
521    f= content (F,F.mvar());
522    return N(d*gcd(G,f));
523  }
524
525  maxNumVars= tmax (getNumVars (F), getNumVars (G));
526  sizeF= size (F);
527  sizeG= size (G);
528
529  if (!isRat)
530    On (SW_RATIONAL);
531  if (sizeF/maxNumVars > 500 && sizeG/maxNumVars > 500)
532  {
533    Off(SW_USE_EZGCD);
534    result=gcd( F, G );
535    On(SW_USE_EZGCD);
536    if (!isRat)
537      Off (SW_RATIONAL);
538    result /= icontent (result);
539    DEBDECLEVEL( cerr, "ezgcd" );
540    return N (d*result);
541  }
542
543  int dummy= 0;
544  if ( gcd_test_one( F, G, false, dummy ) )
545  {
546    DEBDECLEVEL( cerr, "ezgcd" );
547    if (!isRat)
548      Off (SW_RATIONAL);
549    return N (d);
550  }
551  lcF = LC( F, x ); lcG = LC( G, x );
552  lcD = gcd( lcF, lcG );
553  delta = 0;
554  degF = degree( F, x ); degG = degree( G, x );
555  t = tmax( F.level(), G.level() );
556  if ( ! internal )
557    b = REvaluation( 2, t, IntRandom( 25 ) );
558  while ( ! gcdfound )
559  {
560    /// ---> A2
561    DEBOUTLN( cerr, "search for evaluation, delta = " << delta );
562    DEBOUTLN( cerr, "F = " << F );
563    DEBOUTLN( cerr, "G = " << G );
564    TIMING_START (ez_eval)
565    if (!findeval( F, G, Fb, Gb, Db, b, delta, degF, degG, maxeval, count,
566                   o, 25, l))
567    {
568      Off(SW_USE_EZGCD);
569      result=gcd( F, G );
570      On(SW_USE_EZGCD);
571      if (!isRat)
572        Off (SW_RATIONAL);
573      DEBDECLEVEL( cerr, "ezgcd" );
574      result /= icontent (result);
575      return N (d*result);
576    }
577    TIMING_END_AND_PRINT (ez_eval, "time to find eval point in EZ1: ")
578    DEBOUTLN( cerr, "found evaluation b = " << b );
579    DEBOUTLN( cerr, "F(b) = " << Fb );
580    DEBOUTLN( cerr, "G(b) = " << Gb );
581    DEBOUTLN( cerr, "D(b) = " << Db );
582    delta = degree( Db );
583    /// ---> A3
584    if (delta == degF)
585    {
586      if (degF <= degG  && fdivides (F, G))
587      {
588        DEBDECLEVEL( cerr, "ezgcd" );
589        if (!isRat)
590          Off (SW_RATIONAL);
591        return N (d*F);
592      }
593      else
594        delta--;
595    }
596    else if (delta == degG)
597    {
598      if (degG <= degF && fdivides( G, F ))
599      {
600        DEBDECLEVEL( cerr, "ezgcd" );
601        if (!isRat)
602          Off (SW_RATIONAL);
603        return N (d*G);
604      }
605      else
606        delta--;
607    }
608    if ( delta == 0 )
609    {
610      DEBDECLEVEL( cerr, "ezgcd" );
611      if (!isRat)
612        Off (SW_RATIONAL);
613      return N (d);
614    }
615    /// ---> A4
616    //deltaold = delta;
617    while ( 1 )
618    {
619      bt = b;
620      TIMING_START (ez_eval)
621      if (!findeval( F, G, Fbt, Gbt, Dbt, bt, delta, degF, degG, maxeval, count,
622                     o, 25,l ))
623      {
624        Off(SW_USE_EZGCD);
625        result=gcd( F, G );
626        On(SW_USE_EZGCD);
627        if (!isRat)
628          Off (SW_RATIONAL);
629        DEBDECLEVEL( cerr, "ezgcd" );
630        result /= icontent (result);
631        return N (d*result);
632      }
633      TIMING_END_AND_PRINT (ez_eval, "time to find eval point in EZ2: ")
634      int dd=degree( Dbt );
635      if ( dd /*degree( Dbt )*/ == 0 )
636      {
637        DEBDECLEVEL( cerr, "ezgcd" );
638        if (!isRat)
639          Off (SW_RATIONAL);
640        return N (d);
641      }
642      if ( dd /*degree( Dbt )*/ == delta )
643        break;
644      else  if ( dd /*degree( Dbt )*/ < delta )
645      {
646        delta = dd /*degree( Dbt )*/;
647        b = bt;
648        Db = Dbt; Fb = Fbt; Gb = Gbt;
649      }
650      DEBOUTLN( cerr, "now after A4, delta = " << delta );
651      /// ---> A5
652      if (delta == degF)
653      {
654        if (degF <= degG  && fdivides (F, G))
655        {
656          DEBDECLEVEL( cerr, "ezgcd" );
657          if (!isRat)
658            Off (SW_RATIONAL);
659          return N (d*F);
660        }
661        else
662          delta--;
663      }
664      else if (delta == degG)
665      {
666        if (degG <= degF && fdivides( G, F ))
667        {
668          DEBDECLEVEL( cerr, "ezgcd" );
669          if (!isRat)
670            Off (SW_RATIONAL);
671          return N (d*G);
672        }
673        else
674          delta--;
675      }
676      if ( delta == 0 )
677      {
678        DEBDECLEVEL( cerr, "ezgcd" );
679        if (!isRat)
680          Off (SW_RATIONAL);
681        return N (d);
682      }
683    }
684    if ( delta != degF && delta != degG )
685    {
686      /// ---> A6
687      bool B_is_F;
688      CanonicalForm xxx1, xxx2;
689      CanonicalForm buf;
690      DD[1] = Fb / Db;
691      buf= Gb/Db;
692      xxx1 = gcd( DD[1], Db );
693      xxx2 = gcd( buf, Db );
694      if (((xxx1.inCoeffDomain() && xxx2.inCoeffDomain()) &&
695          (size (F) <= size (G)))
696            || (xxx1.inCoeffDomain() && !xxx2.inCoeffDomain()))
697      {
698        B = F;
699        DD[2] = Db;
700        lcDD[1] = lcF;
701        lcDD[2] = lcD;
702        B_is_F = true;
703      }
704      else if (((xxx1.inCoeffDomain() && xxx2.inCoeffDomain()) &&
705                (size (G) < size (F)))
706                || (!xxx1.inCoeffDomain() && xxx2.inCoeffDomain()))
707      {
708        DD[1] = buf;
709        B = G;
710        DD[2] = Db;
711        lcDD[1] = lcG;
712        lcDD[2] = lcD;
713        B_is_F = false;
714      }
715      else
716      {
717        //special case
718        Off(SW_USE_EZGCD);
719        result=gcd( F, G );
720        On(SW_USE_EZGCD);
721        if (!isRat)
722          Off (SW_RATIONAL);
723        DEBDECLEVEL( cerr, "ezgcd" );
724        result /= icontent (result);
725        return N (d*result);
726      }
727      /// ---> A7
728      DD[2] = DD[2] * ( b( lcDD[2] ) / lc( DD[2] ) );
729      DD[1] = DD[1] * ( b( lcDD[1] ) / lc( DD[1] ) );
730      DEBOUTLN( cerr, "(hensel) B    = " << B );
731      DEBOUTLN( cerr, "(hensel) lcB  = " << LC( B, Variable(1) ) );
732      DEBOUTLN( cerr, "(hensel) b(B) = " << b(B) );
733      DEBOUTLN( cerr, "(hensel) DD   = " << DD );
734      DEBOUTLN( cerr, "(hensel) lcDD = " << lcDD );
735      TIMING_START (ez_hensel_lift)
736      gcdfound= Hensel (B*lcD, DD, b, lcDD);
737      TIMING_END_AND_PRINT (ez_hensel_lift, "time to hensel lift in EZ: ")
738      DEBOUTLN( cerr, "(hensel finished) DD   = " << DD );
739
740      if (gcdfound == -1)
741      {
742        Off (SW_USE_EZGCD);
743        result= gcd (F,G);
744        On (SW_USE_EZGCD);
745        if (!isRat)
746          Off (SW_RATIONAL);
747        DEBDECLEVEL( cerr, "ezgcd" );
748        result /= icontent (result);
749        return N (d*result);
750      }
751
752      if (gcdfound)
753      {
754        TIMING_START (ez_termination)
755        contcand= content (DD[2], Variable (1));
756        cand = DD[2] / contcand;
757        if (B_is_F)
758          gcdfound = fdivides( cand, G ) && cand*(DD[1]/(lcD/contcand)) == F;
759        else
760          gcdfound = fdivides( cand, F ) && cand*(DD[1]/(lcD/contcand)) == G;
761        TIMING_END_AND_PRINT (ez_termination,
762                              "time for termination test in EZ: ")
763      }
764      /// ---> A8 (gcdfound)
765    }
766    delta--;
767  }
768  /// ---> A9
769  DEBDECLEVEL( cerr, "ezgcd" );
770  cand *= bCommonDen (cand);
771  if (!isRat)
772    Off (SW_RATIONAL);
773  cand /= icontent (cand);
774  return N (d*cand);
775}
776#endif
777
778CanonicalForm
779ezgcd ( const CanonicalForm & FF, const CanonicalForm & GG )
780{
781#ifdef HAVE_NTL
782  REvaluation b;
783  return ezgcd( FF, GG, b, false );
784#else
785  Off (SW_USE_EZGCD);
786  return gcd (FF, GG);
787  On (SW_USE_EZGCD);
788#endif
789}
790
791static inline
792int Hensel_P (const CanonicalForm & UU, CFArray & G, const Evaluation & AA,
793              const CFArray& LeadCoeffs )
794{
795  CFList factors;
796  factors.append (G[1]);
797  factors.append (G[2]);
798
799  CFMap NN, MM;
800  Evaluation A= optimize4Lift (UU, MM, NN, AA);
801
802  CanonicalForm U= MM (UU);
803  CFArray LCs= CFArray (1,2);
804  LCs [1]= MM (LeadCoeffs [1]);
805  LCs [2]= MM (LeadCoeffs [2]);
806
807  CFList evaluation;
808  long termEstimate= size (U);
809  for (int i= A.min(); i <= A.max(); i++)
810  {
811    if (!A[i].isZero() && (getCharacteristic() > degree (U,i))) //TODO find a good estimate for getCharacteristic() <= degree (U,i)
812    {
813      termEstimate *= degree (U,i)*2;
814      termEstimate /= 3;
815    }
816    evaluation.append (A [i]);
817  }
818  if (termEstimate/getNumVars(U) > 500)
819    return -1;
820  CFList UEval;
821  CanonicalForm shiftedU= myShift2Zero (U, UEval, evaluation);
822
823  if (size (shiftedU)/getNumVars (U) > 500)
824    return -1;
825
826  CFArray shiftedLCs= CFArray (2);
827  CFList shiftedLCsEval1, shiftedLCsEval2;
828  shiftedLCs[0]= myShift2Zero (LCs[1], shiftedLCsEval1, evaluation);
829  shiftedLCs[1]= myShift2Zero (LCs[2], shiftedLCsEval2, evaluation);
830  factors.insert (1);
831  int liftBound= degree (UEval.getLast(), 2) + 1;
832  CFArray Pi;
833  CFMatrix M= CFMatrix (liftBound, factors.length() - 1);
834  CFList diophant;
835  CFArray lcs= CFArray (2);
836  lcs [0]= shiftedLCsEval1.getFirst();
837  lcs [1]= shiftedLCsEval2.getFirst();
838  nonMonicHenselLift12 (UEval.getFirst(), factors, liftBound, Pi, diophant, M,
839                        lcs, false);
840
841  for (CFListIterator i= factors; i.hasItem(); i++)
842  {
843    if (!fdivides (i.getItem(), UEval.getFirst()))
844      return 0;
845  }
846
847  int * liftBounds;
848  bool noOneToOne= false;
849  if (U.level() > 2)
850  {
851    liftBounds= new int [U.level() - 1]; /* index: 0.. U.level()-2 */
852    liftBounds[0]= liftBound;
853    for (int i= 1; i < U.level() - 1; i++)
854      liftBounds[i]= degree (shiftedU, Variable (i + 2)) + 1;
855    factors= nonMonicHenselLift2 (UEval, factors, liftBounds, U.level() - 1,
856                                  false, shiftedLCsEval1, shiftedLCsEval2, Pi,
857                                  diophant, noOneToOne);
858    delete [] liftBounds;
859    if (noOneToOne)
860      return 0;
861  }
862  G[1]= factors.getFirst();
863  G[2]= factors.getLast();
864  G[1]= myReverseShift (G[1], evaluation);
865  G[2]= myReverseShift (G[2], evaluation);
866  G[1]= NN (G[1]);
867  G[2]= NN (G[2]);
868  return 1;
869}
870
871static inline
872bool findeval_P (const CanonicalForm & F, const CanonicalForm & G,
873                 CanonicalForm & Fb, CanonicalForm & Gb, CanonicalForm & Db,
874                 REvaluation & b, int delta, int degF, int degG, int maxeval,
875                 int & count, int& k, int bound, int& l)
876{
877  if( count == 0 && delta != 0)
878  {
879    if( count++ > maxeval )
880      return false;
881  }
882  if (count > 0)
883  {
884    b.nextpoint(k);
885    if (k == 0)
886      k++;
887    l++;
888    if (l > bound)
889    {
890      l= 1;
891      k++;
892      if (k > tmax (F.level(), G.level()) - 1)
893        return false;
894      b.nextpoint (k);
895    }
896    if (count++ > maxeval)
897      return false;
898  }
899  while( true )
900  {
901    Fb = b( F );
902    if( degree( Fb, 1 ) == degF )
903    {
904      Gb = b( G );
905      if( degree( Gb, 1 ) == degG )
906      {
907        Db = gcd( Fb, Gb );
908        if( delta > 0 )
909        {
910          if( degree( Db, 1 ) <= delta )
911            return true;
912        }
913        else
914          return true;
915      }
916    }
917    if (k == 0)
918      k++;
919    b.nextpoint(k);
920    l++;
921    if (l > bound)
922    {
923      l= 1;
924      k++;
925      if (k > tmax (F.level(), G.level()) - 1)
926        return false;
927      b.nextpoint (k);
928    }
929    if( count++ > maxeval )
930      return false;
931  }
932}
933
934// parameters for heuristic
935static int maxNumEval= 200;
936static int sizePerVars1= 500; //try dense gcd if size/#variables is bigger
937
938/// Extended Zassenhaus GCD for finite fields.
939/// In case things become too dense we switch to a modular algorithm
940CanonicalForm EZGCD_P( const CanonicalForm & FF, const CanonicalForm & GG )
941{
942  if (FF.isZero() && degree(GG) > 0) return GG/Lc(GG);
943  else if (GG.isZero() && degree (FF) > 0) return FF/Lc(FF);
944  else if (FF.isZero() && GG.isZero()) return FF.genOne();
945  if (FF.inBaseDomain() || GG.inBaseDomain()) return FF.genOne();
946  if (FF.isUnivariate() && fdivides(FF, GG)) return FF/Lc(FF);
947  if (GG.isUnivariate() && fdivides(GG, FF)) return GG/Lc(GG);
948  if (FF == GG) return FF/Lc(FF);
949
950  int maxNumVars= tmax (getNumVars (FF), getNumVars (GG));
951  Variable a, oldA;
952  int sizeF= size (FF);
953  int sizeG= size (GG);
954
955  if (sizeF/maxNumVars > sizePerVars1 && sizeG/maxNumVars > sizePerVars1)
956  {
957    if (hasFirstAlgVar (FF, a) || hasFirstAlgVar (GG, a))
958      return modGCDFq (FF, GG, a);
959    else if (CFFactory::gettype() == GaloisFieldDomain)
960      return modGCDGF (FF, GG);
961    else
962      return modGCDFp (FF, GG);
963  }
964
965  CanonicalForm F, G, f, g, d, Fb, Gb, Db, Fbt, Gbt, Dbt, B0, B, D0, lcF, lcG,
966                lcD;
967  CFArray DD( 1, 2 ), lcDD( 1, 2 );
968  int degF, degG, delta, count;
969  int maxeval;
970  maxeval= tmin((getCharacteristic()/
971                (int)(ilog2(getCharacteristic())*log2exp))*2, maxNumEval);
972  count= 0; // number of eval. used
973  REvaluation b, bt;
974  int gcdfound = 0;
975  Variable x = Variable(1);
976
977  F= FF;
978  G= GG;
979
980  CFMap M,N;
981  int smallestDegLev;
982  TIMING_START (ez_p_compress)
983  int best_level= compress4EZGCD (F, G, M, N, smallestDegLev);
984
985  if (best_level == 0) return G.genOne();
986
987  F= M (F);
988  G= M (G);
989  TIMING_END_AND_PRINT (ez_p_compress, "time for compression in EZ_P: ")
990
991  TIMING_START (ez_p_content)
992  f = content( F, x ); g = content( G, x );
993  d = gcd( f, g );
994  F /= f; G /= g;
995  TIMING_END_AND_PRINT (ez_p_content, "time to extract content in EZ_P: ")
996
997  if( F.isUnivariate() && G.isUnivariate() )
998  {
999    if( F.mvar() == G.mvar() )
1000      d *= gcd( F, G );
1001    else
1002      return N (d);
1003    return N (d);
1004  }
1005  if ( F.isUnivariate())
1006  {
1007    g= content (G,G.mvar());
1008    return N(d*gcd(F,g));
1009  }
1010  if ( G.isUnivariate())
1011  {
1012    f= content (F,F.mvar());
1013    return N(d*gcd(G,f));
1014  }
1015
1016  maxNumVars= tmax (getNumVars (F), getNumVars (G));
1017  sizeF= size (F);
1018  sizeG= size (G);
1019
1020  if (sizeF/maxNumVars > sizePerVars1 && sizeG/maxNumVars > sizePerVars1)
1021  {
1022    if (hasFirstAlgVar (F, a) || hasFirstAlgVar (G, a))
1023      return N (d*modGCDFq (F, G, a));
1024    else if (CFFactory::gettype() == GaloisFieldDomain)
1025      return N (d*modGCDGF (F, G));
1026    else
1027      return N (d*modGCDFp (F, G));
1028  }
1029
1030  int dummy= 0;
1031  if( gcd_test_one( F, G, false, dummy ) )
1032  {
1033    return N (d);
1034  }
1035
1036  bool passToGF= false;
1037  bool extOfExt= false;
1038  int p= getCharacteristic();
1039  bool algExtension= (hasFirstAlgVar(F,a) || hasFirstAlgVar(G,a));
1040  int k= 1;
1041  CanonicalForm primElem, imPrimElem;
1042  CFList source, dest;
1043  if (p < 50 && CFFactory::gettype() != GaloisFieldDomain && !algExtension)
1044  {
1045    if (p == 2)
1046      setCharacteristic (2, 12, 'Z');
1047    else if (p == 3)
1048      setCharacteristic (3, 4, 'Z');
1049    else if (p == 5 || p == 7)
1050      setCharacteristic (p, 3, 'Z');
1051    else
1052      setCharacteristic (p, 2, 'Z');
1053    passToGF= true;
1054    F= F.mapinto();
1055    G= G.mapinto();
1056    maxeval= 2*ipower (p, getGFDegree());
1057  }
1058  else if (CFFactory::gettype() == GaloisFieldDomain &&
1059           ipower (p , getGFDegree()) < 50)
1060  {
1061    k= getGFDegree();
1062    if (ipower (p, 2*k) > 50)
1063      setCharacteristic (p, 2*k, gf_name);
1064    else
1065      setCharacteristic (p, 3*k, gf_name);
1066    F= GFMapUp (F, k);
1067    G= GFMapUp (G, k);
1068    maxeval= tmin (2*ipower (p, getGFDegree()), maxNumEval);
1069  }
1070  else if (p < 50 && algExtension && CFFactory::gettype() != GaloisFieldDomain)
1071  {
1072    int d= degree (getMipo (a));
1073    oldA= a;
1074    Variable v2;
1075    if (p == 2 && d < 6)
1076    {
1077      if (fac_NTL_char != p)
1078      {
1079        fac_NTL_char= p;
1080        zz_p::init (p);
1081      }
1082      bool primFail= false;
1083      Variable vBuf;
1084      primElem= primitiveElement (a, vBuf, primFail);
1085      ASSERT (!primFail, "failure in integer factorizer");
1086      if (d < 3)
1087      {
1088        zz_pX NTLIrredpoly;
1089        BuildIrred (NTLIrredpoly, d*3);
1090        CanonicalForm newMipo= convertNTLzzpX2CF (NTLIrredpoly, Variable (1));
1091        v2= rootOf (newMipo);
1092      }
1093      else
1094      {
1095        zz_pX NTLIrredpoly;
1096        BuildIrred (NTLIrredpoly, d*2);
1097        CanonicalForm newMipo= convertNTLzzpX2CF (NTLIrredpoly, Variable (1));
1098        v2= rootOf (newMipo);
1099      }
1100      imPrimElem= mapPrimElem (primElem, a, v2);
1101      extOfExt= true;
1102    }
1103    else if ((p == 3 && d < 4) || ((p == 5 || p == 7) && d < 3))
1104    {
1105      if (fac_NTL_char != p)
1106      {
1107        fac_NTL_char= p;
1108        zz_p::init (p);
1109      }
1110      bool primFail= false;
1111      Variable vBuf;
1112      primElem= primitiveElement (a, vBuf, primFail);
1113      ASSERT (!primFail, "failure in integer factorizer");
1114      zz_pX NTLIrredpoly;
1115      BuildIrred (NTLIrredpoly, d*2);
1116      CanonicalForm newMipo= convertNTLzzpX2CF (NTLIrredpoly, Variable (1));
1117      v2= rootOf (newMipo);
1118      imPrimElem= mapPrimElem (primElem, a, v2);
1119      extOfExt= true;
1120    }
1121    if (extOfExt)
1122    {
1123      maxeval= tmin (2*ipower (p, degree (getMipo (v2))), maxNumEval);
1124      F= mapUp (F, a, v2, primElem, imPrimElem, source, dest);
1125      G= mapUp (G, a, v2, primElem, imPrimElem, source, dest);
1126      a= v2;
1127    }
1128  }
1129
1130  lcF = LC( F, x ); lcG = LC( G, x );
1131  lcD = gcd( lcF, lcG );
1132
1133  delta = 0;
1134  degF = degree( F, x ); degG = degree( G, x );
1135
1136  if(hasFirstAlgVar(G,a))
1137    b = REvaluation( 2, tmax(F.level(), G.level()), AlgExtRandomF( a ) );
1138  else
1139  { // both not in extension given by algebraic variable
1140    if (CFFactory::gettype() != GaloisFieldDomain)
1141      b = REvaluation( 2, tmax(F.level(), G.level()), FFRandom() );
1142    else
1143      b = REvaluation( 2, tmax(F.level(), G.level()), GFRandom() );
1144  }
1145
1146  CanonicalForm cand, contcand;
1147  CanonicalForm result;
1148  int o, t;
1149  o= 0;
1150  t= 1;
1151  int goodPointCount= 0;
1152  while( !gcdfound )
1153  {
1154    TIMING_START (ez_p_eval);
1155    if( !findeval_P( F, G, Fb, Gb, Db, b, delta, degF, degG, maxeval, count, o,
1156         maxeval/maxNumVars, t ))
1157    { // too many eval. used --> try another method
1158      Off (SW_USE_EZGCD_P);
1159      result= gcd (F,G);
1160      On (SW_USE_EZGCD_P);
1161      if (passToGF)
1162      {
1163        CanonicalForm mipo= gf_mipo;
1164        setCharacteristic (p);
1165        Variable alpha= rootOf (mipo.mapinto());
1166        result= GF2FalphaRep (result, alpha);
1167      }
1168      if (k > 1)
1169      {
1170        result= GFMapDown (result, k);
1171        setCharacteristic (p, k, gf_name);
1172      }
1173      if (extOfExt)
1174        result= mapDown (result, primElem, imPrimElem, oldA, dest, source);
1175      return N (d*result);
1176    }
1177    TIMING_END_AND_PRINT (ez_p_eval, "time for eval point search in EZ_P1: ");
1178    delta = degree( Db );
1179    if (delta == degF)
1180    {
1181      if (degF <= degG && fdivides (F, G))
1182      {
1183        if (passToGF)
1184        {
1185          CanonicalForm mipo= gf_mipo;
1186          setCharacteristic (p);
1187          Variable alpha= rootOf (mipo.mapinto());
1188          F= GF2FalphaRep (F, alpha);
1189        }
1190        if (k > 1)
1191        {
1192          F= GFMapDown (F, k);
1193          setCharacteristic (p, k, gf_name);
1194        }
1195        if (extOfExt)
1196          F= mapDown (F, primElem, imPrimElem, oldA, dest, source);
1197        return N (d*F);
1198      }
1199      else
1200        delta--;
1201    }
1202    else if (delta == degG)
1203    {
1204      if (degG <= degF && fdivides (G, F))
1205      {
1206        if (passToGF)
1207        {
1208          CanonicalForm mipo= gf_mipo;
1209          setCharacteristic (p);
1210          Variable alpha= rootOf (mipo.mapinto());
1211          G= GF2FalphaRep (G, alpha);
1212        }
1213        if (k > 1)
1214        {
1215          G= GFMapDown (G, k);
1216          setCharacteristic (p, k, gf_name);
1217        }
1218        if (extOfExt)
1219          G= mapDown (G, primElem, imPrimElem, oldA, dest, source);
1220        return N (d*G);
1221      }
1222      else
1223        delta--;
1224    }
1225    if( delta == 0 )
1226    {
1227      if (passToGF)
1228        setCharacteristic (p);
1229      if (k > 1)
1230        setCharacteristic (p, k, gf_name);
1231      return N (d);
1232    }
1233    while( true )
1234    {
1235      bt = b;
1236      TIMING_START (ez_p_eval);
1237      if( !findeval_P(F,G,Fbt,Gbt,Dbt, bt, delta, degF, degG, maxeval, count, o,
1238           maxeval/maxNumVars, t ))
1239      { // too many eval. used --> try another method
1240        Off (SW_USE_EZGCD_P);
1241        result= gcd (F,G);
1242        On (SW_USE_EZGCD_P);
1243        if (passToGF)
1244        {
1245          CanonicalForm mipo= gf_mipo;
1246          setCharacteristic (p);
1247          Variable alpha= rootOf (mipo.mapinto());
1248          result= GF2FalphaRep (result, alpha);
1249        }
1250        if (k > 1)
1251        {
1252          result= GFMapDown (result, k);
1253          setCharacteristic (p, k, gf_name);
1254        }
1255        if (extOfExt)
1256          result= mapDown (result, primElem, imPrimElem, oldA, dest, source);
1257        return N (d*result);
1258      }
1259      TIMING_END_AND_PRINT (ez_p_eval, "time for eval point search in EZ_P2: ");
1260      int dd = degree( Dbt );
1261      if( dd == 0 )
1262      {
1263        if (passToGF)
1264          setCharacteristic (p);
1265        if (k > 1)
1266          setCharacteristic (p, k, gf_name);
1267        return N (d);
1268      }
1269      if( dd == delta )
1270      {
1271        goodPointCount++;
1272        if (goodPointCount == 5)
1273          break;
1274      }
1275      if( dd < delta )
1276      {
1277        goodPointCount= 0;
1278        delta = dd;
1279        b = bt;
1280        Db = Dbt; Fb = Fbt; Gb = Gbt;
1281      }
1282      if (delta == degF)
1283      {
1284        if (degF <= degG && fdivides (F, G))
1285        {
1286          if (passToGF)
1287          {
1288            CanonicalForm mipo= gf_mipo;
1289            setCharacteristic (p);
1290            Variable alpha= rootOf (mipo.mapinto());
1291            F= GF2FalphaRep (F, alpha);
1292          }
1293          if (k > 1)
1294          {
1295            F= GFMapDown (F, k);
1296            setCharacteristic (p, k, gf_name);
1297          }
1298          if (extOfExt)
1299            F= mapDown (F, primElem, imPrimElem, oldA, dest, source);
1300          return N (d*F);
1301        }
1302        else
1303          delta--;
1304      }
1305      else if (delta == degG)
1306      {
1307        if (degG <= degF && fdivides (G, F))
1308        {
1309          if (passToGF)
1310          {
1311            CanonicalForm mipo= gf_mipo;
1312            setCharacteristic (p);
1313            Variable alpha= rootOf (mipo.mapinto());
1314            G= GF2FalphaRep (G, alpha);
1315          }
1316          if (k > 1)
1317          {
1318            G= GFMapDown (G, k);
1319            setCharacteristic (p, k, gf_name);
1320          }
1321          if (extOfExt)
1322            G= mapDown (G, primElem, imPrimElem, oldA, dest, source);
1323          return N (d*G);
1324        }
1325        else
1326          delta--;
1327      }
1328      if( delta == 0 )
1329      {
1330        if (passToGF)
1331          setCharacteristic (p);
1332        if (k > 1)
1333          setCharacteristic (p, k, gf_name);
1334        return N (d);
1335      }
1336    }
1337    if( delta != degF && delta != degG )
1338    {
1339      bool B_is_F;
1340      CanonicalForm xxx1, xxx2;
1341      CanonicalForm buf;
1342      DD[1] = Fb / Db;
1343      buf= Gb/Db;
1344      xxx1 = gcd( DD[1], Db );
1345      xxx2 = gcd( buf, Db );
1346      if (((xxx1.inCoeffDomain() && xxx2.inCoeffDomain()) &&
1347          (size (F) <= size (G)))
1348          || (xxx1.inCoeffDomain() && !xxx2.inCoeffDomain()))
1349      {
1350        B = F;
1351        DD[2] = Db;
1352        lcDD[1] = lcF;
1353        lcDD[2] = lcD;
1354        B_is_F = true;
1355      }
1356      else if (((xxx1.inCoeffDomain() && xxx2.inCoeffDomain()) &&
1357               (size (G) < size (F)))
1358               || (!xxx1.inCoeffDomain() && xxx2.inCoeffDomain()))
1359      {
1360        DD[1] = buf;
1361        B = G;
1362        DD[2] = Db;
1363        lcDD[1] = lcG;
1364        lcDD[2] = lcD;
1365        B_is_F = false;
1366      }
1367      else // special case handling
1368      {
1369        Off (SW_USE_EZGCD_P);
1370        result= gcd (F,G);
1371        On (SW_USE_EZGCD_P);
1372        if (passToGF)
1373        {
1374          CanonicalForm mipo= gf_mipo;
1375          setCharacteristic (p);
1376          Variable alpha= rootOf (mipo.mapinto());
1377          result= GF2FalphaRep (result, alpha);
1378        }
1379        if (k > 1)
1380        {
1381          result= GFMapDown (result, k);
1382          setCharacteristic (p, k, gf_name);
1383        }
1384        if (extOfExt)
1385          result= mapDown (result, primElem, imPrimElem, oldA, dest, source);
1386        return N (d*result);
1387      }
1388      DD[2] = DD[2] * ( b( lcDD[2] ) / lc( DD[2] ) );
1389      DD[1] = DD[1] * ( b( lcDD[1] ) / lc( DD[1] ) );
1390
1391      if (size (B*lcDD[2])/maxNumVars > sizePerVars1)
1392      {
1393        if (algExtension)
1394        {
1395          result= modGCDFq (F, G, a);
1396          if (extOfExt)
1397            result= mapDown (result, primElem, imPrimElem, oldA, dest, source);
1398          return N (d*result);
1399        }
1400        if (CFFactory::gettype() == GaloisFieldDomain)
1401        {
1402          result= modGCDGF (F, G);
1403          if (passToGF)
1404          {
1405            CanonicalForm mipo= gf_mipo;
1406            setCharacteristic (p);
1407            Variable alpha= rootOf (mipo.mapinto());
1408            result= GF2FalphaRep (result, alpha);
1409          }
1410          if (k > 1)
1411          {
1412            result= GFMapDown (result, k);
1413            setCharacteristic (p, k, gf_name);
1414          }
1415          return N (d*result);
1416        }
1417        else
1418          return N (d*modGCDFp (F,G));
1419      }
1420
1421      TIMING_START (ez_p_hensel_lift);
1422      gcdfound= Hensel_P (B*lcD, DD, b, lcDD);
1423      TIMING_END_AND_PRINT (ez_p_hensel_lift, "time for Hensel lift in EZ_P: ");
1424
1425      if (gcdfound == -1) //things became dense
1426      {
1427        if (algExtension)
1428        {
1429          result= modGCDFq (F, G, a);
1430          if (extOfExt)
1431            result= mapDown (result, primElem, imPrimElem, oldA, dest, source);
1432          return N (d*result);
1433        }
1434        if (CFFactory::gettype() == GaloisFieldDomain)
1435        {
1436          result= modGCDGF (F, G);
1437          if (passToGF)
1438          {
1439            CanonicalForm mipo= gf_mipo;
1440            setCharacteristic (p);
1441            Variable alpha= rootOf (mipo.mapinto());
1442            result= GF2FalphaRep (result, alpha);
1443          }
1444          if (k > 1)
1445          {
1446            result= GFMapDown (result, k);
1447            setCharacteristic (p, k, gf_name);
1448          }
1449          return N (d*result);
1450        }
1451        else
1452        {
1453          if (p >= cf_getBigPrime(0))
1454            return N (d*sparseGCDFp (F,G));
1455          else
1456            return N (d*modGCDFp (F,G));
1457        }
1458      }
1459
1460      if (gcdfound == 1)
1461      {
1462        TIMING_START (termination_test);
1463        contcand= content (DD[2], Variable (1));
1464        cand = DD[2] / contcand;
1465        if (B_is_F)
1466          gcdfound = fdivides( cand, G ) && cand*(DD[1]/(lcD/contcand)) == F;
1467        else
1468          gcdfound = fdivides( cand, F ) && cand*(DD[1]/(lcD/contcand)) == G;
1469        TIMING_END_AND_PRINT (termination_test,
1470                              "time for termination test EZ_P: ");
1471
1472        if (passToGF && gcdfound)
1473        {
1474          CanonicalForm mipo= gf_mipo;
1475          setCharacteristic (p);
1476          Variable alpha= rootOf (mipo.mapinto());
1477          cand= GF2FalphaRep (cand, alpha);
1478        }
1479        if (k > 1 && gcdfound)
1480        {
1481          cand= GFMapDown (cand, k);
1482          setCharacteristic (p, k, gf_name);
1483        }
1484        if (extOfExt && gcdfound)
1485          cand= mapDown (cand, primElem, imPrimElem, oldA, dest, source);
1486      }
1487    }
1488    delta--;
1489    goodPointCount= 0;
1490  }
1491  return N (d*cand);
1492}
1493
Note: See TracBrowser for help on using the repository browser.