source: git/factory/cfEzgcd.cc @ 1950d7

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