source: git/factory/cfEzgcd.cc @ 03c742

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