source: git/factory/cfEzgcd.cc @ 8861479

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