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

spielwiese
Last change on this file since 03f640 was 03f640, checked in by Martin Lee <martinlee84@…>, 10 years ago
chg: added pruning of alg extensions
  • Property mode set to 100644
File size: 37.5 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 "cfModGcd.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
791#ifdef HAVE_NTL
792static inline
793int Hensel_P (const CanonicalForm & UU, CFArray & G, const Evaluation & AA,
794              const CFArray& LeadCoeffs )
795{
796  CFList factors;
797  factors.append (G[1]);
798  factors.append (G[2]);
799
800  CFMap NN, MM;
801  Evaluation A= optimize4Lift (UU, MM, NN, AA);
802
803  CanonicalForm U= MM (UU);
804  CFArray LCs= CFArray (1,2);
805  LCs [1]= MM (LeadCoeffs [1]);
806  LCs [2]= MM (LeadCoeffs [2]);
807
808  CFList evaluation;
809  long termEstimate= size (U);
810  for (int i= A.min(); i <= A.max(); i++)
811  {
812    if (!A[i].isZero() && (getCharacteristic() > degree (U,i))) //TODO find a good estimate for getCharacteristic() <= degree (U,i)
813    {
814      termEstimate *= degree (U,i)*2;
815      termEstimate /= 3;
816    }
817    evaluation.append (A [i]);
818  }
819  if (termEstimate/getNumVars(U) > 500)
820    return -1;
821  CFList UEval;
822  CanonicalForm shiftedU= myShift2Zero (U, UEval, evaluation);
823
824  if (size (shiftedU)/getNumVars (U) > 500)
825    return -1;
826
827  CFArray shiftedLCs= CFArray (2);
828  CFList shiftedLCsEval1, shiftedLCsEval2;
829  shiftedLCs[0]= myShift2Zero (LCs[1], shiftedLCsEval1, evaluation);
830  shiftedLCs[1]= myShift2Zero (LCs[2], shiftedLCsEval2, evaluation);
831  factors.insert (1);
832  int liftBound= degree (UEval.getLast(), 2) + 1;
833  CFArray Pi;
834  CFMatrix M= CFMatrix (liftBound, factors.length() - 1);
835  CFList diophant;
836  CFArray lcs= CFArray (2);
837  lcs [0]= shiftedLCsEval1.getFirst();
838  lcs [1]= shiftedLCsEval2.getFirst();
839  nonMonicHenselLift12 (UEval.getFirst(), factors, liftBound, Pi, diophant, M,
840                        lcs, false);
841
842  for (CFListIterator i= factors; i.hasItem(); i++)
843  {
844    if (!fdivides (i.getItem(), UEval.getFirst()))
845      return 0;
846  }
847
848  int * liftBounds;
849  bool noOneToOne= false;
850  if (U.level() > 2)
851  {
852    liftBounds= new int [U.level() - 1]; /* index: 0.. U.level()-2 */
853    liftBounds[0]= liftBound;
854    for (int i= 1; i < U.level() - 1; i++)
855      liftBounds[i]= degree (shiftedU, Variable (i + 2)) + 1;
856    factors= nonMonicHenselLift2 (UEval, factors, liftBounds, U.level() - 1,
857                                  false, shiftedLCsEval1, shiftedLCsEval2, Pi,
858                                  diophant, noOneToOne);
859    delete [] liftBounds;
860    if (noOneToOne)
861      return 0;
862  }
863  G[1]= factors.getFirst();
864  G[2]= factors.getLast();
865  G[1]= myReverseShift (G[1], evaluation);
866  G[2]= myReverseShift (G[2], evaluation);
867  G[1]= NN (G[1]);
868  G[2]= NN (G[2]);
869  return 1;
870}
871
872static inline
873bool findeval_P (const CanonicalForm & F, const CanonicalForm & G,
874                 CanonicalForm & Fb, CanonicalForm & Gb, CanonicalForm & Db,
875                 REvaluation & b, int delta, int degF, int degG, int maxeval,
876                 int & count, int& k, int bound, int& l)
877{
878  if( count == 0 && delta != 0)
879  {
880    if( count++ > maxeval )
881      return false;
882  }
883  if (count > 0)
884  {
885    b.nextpoint(k);
886    if (k == 0)
887      k++;
888    l++;
889    if (l > bound)
890    {
891      l= 1;
892      k++;
893      if (k > tmax (F.level(), G.level()) - 1)
894        return false;
895      b.nextpoint (k);
896    }
897    if (count++ > maxeval)
898      return false;
899  }
900  while( true )
901  {
902    Fb = b( F );
903    if( degree( Fb, 1 ) == degF )
904    {
905      Gb = b( G );
906      if( degree( Gb, 1 ) == degG )
907      {
908        Db = gcd( Fb, Gb );
909        if( delta > 0 )
910        {
911          if( degree( Db, 1 ) <= delta )
912            return true;
913        }
914        else
915          return true;
916      }
917    }
918    if (k == 0)
919      k++;
920    b.nextpoint(k);
921    l++;
922    if (l > bound)
923    {
924      l= 1;
925      k++;
926      if (k > tmax (F.level(), G.level()) - 1)
927        return false;
928      b.nextpoint (k);
929    }
930    if( count++ > maxeval )
931      return false;
932  }
933}
934
935// parameters for heuristic
936static int maxNumEval= 200;
937static int sizePerVars1= 500; //try dense gcd if size/#variables is bigger
938
939/// Extended Zassenhaus GCD for finite fields.
940/// In case things become too dense we switch to a modular algorithm
941CanonicalForm EZGCD_P( const CanonicalForm & FF, const CanonicalForm & GG )
942{
943  if (FF.isZero() && degree(GG) > 0) return GG/Lc(GG);
944  else if (GG.isZero() && degree (FF) > 0) return FF/Lc(FF);
945  else if (FF.isZero() && GG.isZero()) return FF.genOne();
946  if (FF.inBaseDomain() || GG.inBaseDomain()) return FF.genOne();
947  if (FF.isUnivariate() && fdivides(FF, GG)) return FF/Lc(FF);
948  if (GG.isUnivariate() && fdivides(GG, FF)) return GG/Lc(GG);
949  if (FF == GG) return FF/Lc(FF);
950
951  int maxNumVars= tmax (getNumVars (FF), getNumVars (GG));
952  Variable a, oldA;
953  int sizeF= size (FF);
954  int sizeG= size (GG);
955
956  if (sizeF/maxNumVars > sizePerVars1 && sizeG/maxNumVars > sizePerVars1)
957  {
958    if (hasFirstAlgVar (FF, a) || hasFirstAlgVar (GG, a))
959      return modGCDFq (FF, GG, a);
960    else if (CFFactory::gettype() == GaloisFieldDomain)
961      return modGCDGF (FF, GG);
962    else
963      return modGCDFp (FF, GG);
964  }
965
966  CanonicalForm F, G, f, g, d, Fb, Gb, Db, Fbt, Gbt, Dbt, B0, B, D0, lcF, lcG,
967                lcD;
968  CFArray DD( 1, 2 ), lcDD( 1, 2 );
969  int degF, degG, delta, count;
970  int maxeval;
971  maxeval= tmin((getCharacteristic()/
972                (int)(ilog2(getCharacteristic())*log2exp))*2, maxNumEval);
973  count= 0; // number of eval. used
974  REvaluation b, bt;
975  int gcdfound = 0;
976  Variable x = Variable(1);
977
978  F= FF;
979  G= GG;
980
981  CFMap M,N;
982  int smallestDegLev;
983  TIMING_START (ez_p_compress)
984  int best_level= compress4EZGCD (F, G, M, N, smallestDegLev);
985
986  if (best_level == 0) return G.genOne();
987
988  F= M (F);
989  G= M (G);
990  TIMING_END_AND_PRINT (ez_p_compress, "time for compression in EZ_P: ")
991
992  TIMING_START (ez_p_content)
993  f = content( F, x ); g = content( G, x );
994  d = gcd( f, g );
995  F /= f; G /= g;
996  TIMING_END_AND_PRINT (ez_p_content, "time to extract content in EZ_P: ")
997
998  if( F.isUnivariate() && G.isUnivariate() )
999  {
1000    if( F.mvar() == G.mvar() )
1001      d *= gcd( F, G );
1002    else
1003      return N (d);
1004    return N (d);
1005  }
1006  if ( F.isUnivariate())
1007  {
1008    g= content (G,G.mvar());
1009    return N(d*gcd(F,g));
1010  }
1011  if ( G.isUnivariate())
1012  {
1013    f= content (F,F.mvar());
1014    return N(d*gcd(G,f));
1015  }
1016
1017  maxNumVars= tmax (getNumVars (F), getNumVars (G));
1018  sizeF= size (F);
1019  sizeG= size (G);
1020
1021  if (sizeF/maxNumVars > sizePerVars1 && sizeG/maxNumVars > sizePerVars1)
1022  {
1023    if (hasFirstAlgVar (F, a) || hasFirstAlgVar (G, a))
1024      return N (d*modGCDFq (F, G, a));
1025    else if (CFFactory::gettype() == GaloisFieldDomain)
1026      return N (d*modGCDGF (F, G));
1027    else
1028      return N (d*modGCDFp (F, G));
1029  }
1030
1031  int dummy= 0;
1032  if( gcd_test_one( F, G, false, dummy ) )
1033  {
1034    return N (d);
1035  }
1036
1037  bool passToGF= false;
1038  bool extOfExt= false;
1039  int p= getCharacteristic();
1040  bool algExtension= (hasFirstAlgVar(F,a) || hasFirstAlgVar(G,a));
1041  int k= 1;
1042  CanonicalForm primElem, imPrimElem;
1043  CFList source, dest;
1044  if (p < 50 && CFFactory::gettype() != GaloisFieldDomain && !algExtension)
1045  {
1046    if (p == 2)
1047      setCharacteristic (2, 12, 'Z');
1048    else if (p == 3)
1049      setCharacteristic (3, 4, 'Z');
1050    else if (p == 5 || p == 7)
1051      setCharacteristic (p, 3, 'Z');
1052    else
1053      setCharacteristic (p, 2, 'Z');
1054    passToGF= true;
1055    F= F.mapinto();
1056    G= G.mapinto();
1057    maxeval= 2*ipower (p, getGFDegree());
1058  }
1059  else if (CFFactory::gettype() == GaloisFieldDomain &&
1060           ipower (p , getGFDegree()) < 50)
1061  {
1062    k= getGFDegree();
1063    if (ipower (p, 2*k) > 50)
1064      setCharacteristic (p, 2*k, gf_name);
1065    else
1066      setCharacteristic (p, 3*k, gf_name);
1067    F= GFMapUp (F, k);
1068    G= GFMapUp (G, k);
1069    maxeval= tmin (2*ipower (p, getGFDegree()), maxNumEval);
1070  }
1071  else if (p < 50 && algExtension && CFFactory::gettype() != GaloisFieldDomain)
1072  {
1073    int d= degree (getMipo (a));
1074    oldA= a;
1075    Variable v2;
1076    if (p == 2 && d < 6)
1077    {
1078      if (fac_NTL_char != p)
1079      {
1080        fac_NTL_char= p;
1081        zz_p::init (p);
1082      }
1083      bool primFail= false;
1084      Variable vBuf;
1085      primElem= primitiveElement (a, vBuf, primFail);
1086      ASSERT (!primFail, "failure in integer factorizer");
1087      if (d < 3)
1088      {
1089        zz_pX NTLIrredpoly;
1090        BuildIrred (NTLIrredpoly, d*3);
1091        CanonicalForm newMipo= convertNTLzzpX2CF (NTLIrredpoly, Variable (1));
1092        v2= rootOf (newMipo);
1093      }
1094      else
1095      {
1096        zz_pX NTLIrredpoly;
1097        BuildIrred (NTLIrredpoly, d*2);
1098        CanonicalForm newMipo= convertNTLzzpX2CF (NTLIrredpoly, Variable (1));
1099        v2= rootOf (newMipo);
1100      }
1101      imPrimElem= mapPrimElem (primElem, a, v2);
1102      extOfExt= true;
1103    }
1104    else if ((p == 3 && d < 4) || ((p == 5 || p == 7) && d < 3))
1105    {
1106      if (fac_NTL_char != p)
1107      {
1108        fac_NTL_char= p;
1109        zz_p::init (p);
1110      }
1111      bool primFail= false;
1112      Variable vBuf;
1113      primElem= primitiveElement (a, vBuf, primFail);
1114      ASSERT (!primFail, "failure in integer factorizer");
1115      zz_pX NTLIrredpoly;
1116      BuildIrred (NTLIrredpoly, d*2);
1117      CanonicalForm newMipo= convertNTLzzpX2CF (NTLIrredpoly, Variable (1));
1118      v2= rootOf (newMipo);
1119      imPrimElem= mapPrimElem (primElem, a, v2);
1120      extOfExt= true;
1121    }
1122    if (extOfExt)
1123    {
1124      maxeval= tmin (2*ipower (p, degree (getMipo (v2))), maxNumEval);
1125      F= mapUp (F, a, v2, primElem, imPrimElem, source, dest);
1126      G= mapUp (G, a, v2, primElem, imPrimElem, source, dest);
1127      a= v2;
1128    }
1129  }
1130
1131  lcF = LC( F, x ); lcG = LC( G, x );
1132  lcD = gcd( lcF, lcG );
1133
1134  delta = 0;
1135  degF = degree( F, x ); degG = degree( G, x );
1136
1137  if (algExtension)
1138    b = REvaluation( 2, tmax(F.level(), G.level()), AlgExtRandomF( a ) );
1139  else
1140  { // both not in extension given by algebraic variable
1141    if (CFFactory::gettype() != GaloisFieldDomain)
1142      b = REvaluation( 2, tmax(F.level(), G.level()), FFRandom() );
1143    else
1144      b = REvaluation( 2, tmax(F.level(), G.level()), GFRandom() );
1145  }
1146
1147  CanonicalForm cand, contcand;
1148  CanonicalForm result;
1149  int o, t;
1150  o= 0;
1151  t= 1;
1152  int goodPointCount= 0;
1153  while( !gcdfound )
1154  {
1155    TIMING_START (ez_p_eval);
1156    if( !findeval_P( F, G, Fb, Gb, Db, b, delta, degF, degG, maxeval, count, o,
1157         maxeval/maxNumVars, t ))
1158    { // too many eval. used --> try another method
1159      Off (SW_USE_EZGCD_P);
1160      result= gcd (F,G);
1161      On (SW_USE_EZGCD_P);
1162      if (passToGF)
1163      {
1164        CanonicalForm mipo= gf_mipo;
1165        setCharacteristic (p);
1166        Variable alpha= rootOf (mipo.mapinto());
1167        result= GF2FalphaRep (result, alpha);
1168        prune (alpha);
1169      }
1170      if (k > 1)
1171      {
1172        result= GFMapDown (result, k);
1173        setCharacteristic (p, k, gf_name);
1174      }
1175      if (extOfExt)
1176      {
1177        result= mapDown (result, primElem, imPrimElem, oldA, dest, source);
1178        prune1 (oldA);
1179      }
1180      return N (d*result);
1181    }
1182    TIMING_END_AND_PRINT (ez_p_eval, "time for eval point search in EZ_P1: ");
1183    delta = degree( Db );
1184    if (delta == degF)
1185    {
1186      if (degF <= degG && fdivides (F, G))
1187      {
1188        if (passToGF)
1189        {
1190          CanonicalForm mipo= gf_mipo;
1191          setCharacteristic (p);
1192          Variable alpha= rootOf (mipo.mapinto());
1193          F= GF2FalphaRep (F, alpha);
1194          prune (alpha);
1195        }
1196        if (k > 1)
1197        {
1198          F= GFMapDown (F, k);
1199          setCharacteristic (p, k, gf_name);
1200        }
1201        if (extOfExt)
1202        {
1203          F= mapDown (F, primElem, imPrimElem, oldA, dest, source);
1204          prune1 (oldA);
1205        }
1206        return N (d*F);
1207      }
1208      else
1209        delta--;
1210    }
1211    else if (delta == degG)
1212    {
1213      if (degG <= degF && fdivides (G, F))
1214      {
1215        if (passToGF)
1216        {
1217          CanonicalForm mipo= gf_mipo;
1218          setCharacteristic (p);
1219          Variable alpha= rootOf (mipo.mapinto());
1220          G= GF2FalphaRep (G, alpha);
1221          prune (alpha);
1222        }
1223        if (k > 1)
1224        {
1225          G= GFMapDown (G, k);
1226          setCharacteristic (p, k, gf_name);
1227        }
1228        if (extOfExt)
1229        {
1230          G= mapDown (G, primElem, imPrimElem, oldA, dest, source);
1231          prune1 (oldA);
1232        }
1233        return N (d*G);
1234      }
1235      else
1236        delta--;
1237    }
1238    if( delta == 0 )
1239    {
1240      if (passToGF)
1241        setCharacteristic (p);
1242      if (k > 1)
1243        setCharacteristic (p, k, gf_name);
1244      return N (d);
1245    }
1246    while( true )
1247    {
1248      bt = b;
1249      TIMING_START (ez_p_eval);
1250      if( !findeval_P(F,G,Fbt,Gbt,Dbt, bt, delta, degF, degG, maxeval, count, o,
1251           maxeval/maxNumVars, t ))
1252      { // too many eval. used --> try another method
1253        Off (SW_USE_EZGCD_P);
1254        result= gcd (F,G);
1255        On (SW_USE_EZGCD_P);
1256        if (passToGF)
1257        {
1258          CanonicalForm mipo= gf_mipo;
1259          setCharacteristic (p);
1260          Variable alpha= rootOf (mipo.mapinto());
1261          result= GF2FalphaRep (result, alpha);
1262          prune (alpha);
1263        }
1264        if (k > 1)
1265        {
1266          result= GFMapDown (result, k);
1267          setCharacteristic (p, k, gf_name);
1268        }
1269        if (extOfExt)
1270        {
1271          result= mapDown (result, primElem, imPrimElem, oldA, dest, source);
1272          prune1 (oldA);
1273        }
1274        return N (d*result);
1275      }
1276      TIMING_END_AND_PRINT (ez_p_eval, "time for eval point search in EZ_P2: ");
1277      int dd = degree( Dbt );
1278      if( dd == 0 )
1279      {
1280        if (passToGF)
1281          setCharacteristic (p);
1282        if (k > 1)
1283          setCharacteristic (p, k, gf_name);
1284        return N (d);
1285      }
1286      if( dd == delta )
1287      {
1288        goodPointCount++;
1289        if (goodPointCount == 5)
1290          break;
1291      }
1292      if( dd < delta )
1293      {
1294        goodPointCount= 0;
1295        delta = dd;
1296        b = bt;
1297        Db = Dbt; Fb = Fbt; Gb = Gbt;
1298      }
1299      if (delta == degF)
1300      {
1301        if (degF <= degG && fdivides (F, G))
1302        {
1303          if (passToGF)
1304          {
1305            CanonicalForm mipo= gf_mipo;
1306            setCharacteristic (p);
1307            Variable alpha= rootOf (mipo.mapinto());
1308            F= GF2FalphaRep (F, alpha);
1309            prune (alpha);
1310          }
1311          if (k > 1)
1312          {
1313            F= GFMapDown (F, k);
1314            setCharacteristic (p, k, gf_name);
1315          }
1316          if (extOfExt)
1317          {
1318            F= mapDown (F, primElem, imPrimElem, oldA, dest, source);
1319            prune1 (oldA);
1320          }
1321          return N (d*F);
1322        }
1323        else
1324          delta--;
1325      }
1326      else if (delta == degG)
1327      {
1328        if (degG <= degF && fdivides (G, F))
1329        {
1330          if (passToGF)
1331          {
1332            CanonicalForm mipo= gf_mipo;
1333            setCharacteristic (p);
1334            Variable alpha= rootOf (mipo.mapinto());
1335            G= GF2FalphaRep (G, alpha);
1336            prune (alpha);
1337          }
1338          if (k > 1)
1339          {
1340            G= GFMapDown (G, k);
1341            setCharacteristic (p, k, gf_name);
1342          }
1343          if (extOfExt)
1344          {
1345            G= mapDown (G, primElem, imPrimElem, oldA, dest, source);
1346            prune1 (oldA);
1347          }
1348          return N (d*G);
1349        }
1350        else
1351          delta--;
1352      }
1353      if( delta == 0 )
1354      {
1355        if (passToGF)
1356          setCharacteristic (p);
1357        if (k > 1)
1358          setCharacteristic (p, k, gf_name);
1359        return N (d);
1360      }
1361    }
1362    if( delta != degF && delta != degG )
1363    {
1364      bool B_is_F;
1365      CanonicalForm xxx1, xxx2;
1366      CanonicalForm buf;
1367      DD[1] = Fb / Db;
1368      buf= Gb/Db;
1369      xxx1 = gcd( DD[1], Db );
1370      xxx2 = gcd( buf, Db );
1371      if (((xxx1.inCoeffDomain() && xxx2.inCoeffDomain()) &&
1372          (size (F) <= size (G)))
1373          || (xxx1.inCoeffDomain() && !xxx2.inCoeffDomain()))
1374      {
1375        B = F;
1376        DD[2] = Db;
1377        lcDD[1] = lcF;
1378        lcDD[2] = lcD;
1379        B_is_F = true;
1380      }
1381      else if (((xxx1.inCoeffDomain() && xxx2.inCoeffDomain()) &&
1382               (size (G) < size (F)))
1383               || (!xxx1.inCoeffDomain() && xxx2.inCoeffDomain()))
1384      {
1385        DD[1] = buf;
1386        B = G;
1387        DD[2] = Db;
1388        lcDD[1] = lcG;
1389        lcDD[2] = lcD;
1390        B_is_F = false;
1391      }
1392      else // special case handling
1393      {
1394        Off (SW_USE_EZGCD_P);
1395        result= gcd (F,G);
1396        On (SW_USE_EZGCD_P);
1397        if (passToGF)
1398        {
1399          CanonicalForm mipo= gf_mipo;
1400          setCharacteristic (p);
1401          Variable alpha= rootOf (mipo.mapinto());
1402          result= GF2FalphaRep (result, alpha);
1403          prune (alpha);
1404        }
1405        if (k > 1)
1406        {
1407          result= GFMapDown (result, k);
1408          setCharacteristic (p, k, gf_name);
1409        }
1410        if (extOfExt)
1411        {
1412          result= mapDown (result, primElem, imPrimElem, oldA, dest, source);
1413          prune1 (oldA);
1414        }
1415        return N (d*result);
1416      }
1417      DD[2] = DD[2] * ( b( lcDD[2] ) / lc( DD[2] ) );
1418      DD[1] = DD[1] * ( b( lcDD[1] ) / lc( DD[1] ) );
1419
1420      if (size (B*lcDD[2])/maxNumVars > sizePerVars1)
1421      {
1422        if (algExtension)
1423        {
1424          result= modGCDFq (F, G, a);
1425          if (extOfExt)
1426          {
1427            result= mapDown (result, primElem, imPrimElem, oldA, dest, source);
1428            prune1 (oldA);
1429          }
1430          return N (d*result);
1431        }
1432        if (CFFactory::gettype() == GaloisFieldDomain)
1433        {
1434          result= modGCDGF (F, G);
1435          if (passToGF)
1436          {
1437            CanonicalForm mipo= gf_mipo;
1438            setCharacteristic (p);
1439            Variable alpha= rootOf (mipo.mapinto());
1440            result= GF2FalphaRep (result, alpha);
1441            prune (alpha);
1442          }
1443          if (k > 1)
1444          {
1445            result= GFMapDown (result, k);
1446            setCharacteristic (p, k, gf_name);
1447          }
1448          return N (d*result);
1449        }
1450        else
1451          return N (d*modGCDFp (F,G));
1452      }
1453
1454      TIMING_START (ez_p_hensel_lift);
1455      gcdfound= Hensel_P (B*lcD, DD, b, lcDD);
1456      TIMING_END_AND_PRINT (ez_p_hensel_lift, "time for Hensel lift in EZ_P: ");
1457
1458      if (gcdfound == -1) //things became dense
1459      {
1460        if (algExtension)
1461        {
1462          result= modGCDFq (F, G, a);
1463          if (extOfExt)
1464          {
1465            result= mapDown (result, primElem, imPrimElem, oldA, dest, source);
1466            prune1 (oldA);
1467          }
1468          return N (d*result);
1469        }
1470        if (CFFactory::gettype() == GaloisFieldDomain)
1471        {
1472          result= modGCDGF (F, G);
1473          if (passToGF)
1474          {
1475            CanonicalForm mipo= gf_mipo;
1476            setCharacteristic (p);
1477            Variable alpha= rootOf (mipo.mapinto());
1478            result= GF2FalphaRep (result, alpha);
1479            prune (alpha);
1480          }
1481          if (k > 1)
1482          {
1483            result= GFMapDown (result, k);
1484            setCharacteristic (p, k, gf_name);
1485          }
1486          return N (d*result);
1487        }
1488        else
1489        {
1490          if (p >= cf_getBigPrime(0))
1491            return N (d*sparseGCDFp (F,G));
1492          else
1493            return N (d*modGCDFp (F,G));
1494        }
1495      }
1496
1497      if (gcdfound == 1)
1498      {
1499        TIMING_START (termination_test);
1500        contcand= content (DD[2], Variable (1));
1501        cand = DD[2] / contcand;
1502        if (B_is_F)
1503          gcdfound = fdivides( cand, G ) && cand*(DD[1]/(lcD/contcand)) == F;
1504        else
1505          gcdfound = fdivides( cand, F ) && cand*(DD[1]/(lcD/contcand)) == G;
1506        TIMING_END_AND_PRINT (termination_test,
1507                              "time for termination test EZ_P: ");
1508
1509        if (passToGF && gcdfound)
1510        {
1511          CanonicalForm mipo= gf_mipo;
1512          setCharacteristic (p);
1513          Variable alpha= rootOf (mipo.mapinto());
1514          cand= GF2FalphaRep (cand, alpha);
1515          prune (alpha);
1516        }
1517        if (k > 1 && gcdfound)
1518        {
1519          cand= GFMapDown (cand, k);
1520          setCharacteristic (p, k, gf_name);
1521        }
1522        if (extOfExt && gcdfound)
1523        {
1524          cand= mapDown (cand, primElem, imPrimElem, oldA, dest, source);
1525          prune1 (oldA);
1526        }
1527      }
1528    }
1529    delta--;
1530    goodPointCount= 0;
1531  }
1532  return N (d*cand);
1533}
1534#endif
1535
Note: See TracBrowser for help on using the repository browser.