source: git/factory/cfEzgcd.cc @ 5d8c8b

spielwiese
Last change on this file since 5d8c8b was 5d8c8b, checked in by Hans Schoenemann <hannes@…>, 5 years ago
opt: handle gcd of monomials in ezgcd/EZGCD_P
  • Property mode set to 100644
File size: 35.6 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    return gcd_mon( FF, GG );
494  }
495  else if (sizeG==1)
496  {
497    return gcd_mon( GG, FF );
498  }
499  if (!isRat)
500    On (SW_RATIONAL);
501  if (sizeF/maxNumVars > 500 && sizeG/maxNumVars > 500)
502  {
503    Off(SW_USE_EZGCD);
504    CanonicalForm result=gcd( FF, GG );
505    On(SW_USE_EZGCD);
506    if (!isRat)
507      Off (SW_RATIONAL);
508    result /= icontent (result);
509    DEBDECLEVEL( cerr, "ezgcd" );
510    return result;
511  }
512
513
514  CanonicalForm F, G, f, g, d, Fb, Gb, Db, Fbt, Gbt, Dbt, B0, B, D0, lcF, lcG,
515                lcD, cand, contcand, result;
516  CFArray DD( 1, 2 ), lcDD( 1, 2 );
517  int degF, degG, delta, t, count, maxeval;
518  REvaluation bt;
519  int gcdfound = 0;
520  Variable x = Variable(1);
521  count= 0;
522  maxeval= 200;
523  int o, l;
524  o= 0;
525  l= 1;
526
527  if (!isRat)
528    On (SW_RATIONAL);
529  F= FF*bCommonDen (FF);
530  G= GG*bCommonDen (GG);
531  if (!isRat)
532    Off (SW_RATIONAL);
533
534  TIMING_START (ez_compress)
535  CFMap M,N;
536  int smallestDegLev;
537  int best_level= compress4EZGCD (F, G, M, N, smallestDegLev);
538
539  if (best_level == 0)
540  {
541    DEBDECLEVEL( cerr, "ezgcd" );
542    return G.genOne();
543  }
544
545  F= M (F);
546  G= M (G);
547  TIMING_END_AND_PRINT (ez_compress, "time for compression in EZ: ")
548
549  DEBINCLEVEL( cerr, "ezgcd" );
550  DEBOUTLN( cerr, "FF = " << FF );
551  DEBOUTLN( cerr, "GG = " << GG );
552  TIMING_START (ez_content)
553  f = content( F, x ); g = content( G, x ); d = gcd( f, g );
554  DEBOUTLN( cerr, "f = " << f );
555  DEBOUTLN( cerr, "g = " << g );
556  F /= f; G /= g;
557  TIMING_END_AND_PRINT (ez_content, "time to extract content in EZ: ")
558  if ( F.isUnivariate() )
559  {
560    if ( G.isUnivariate() )
561    {
562      DEBDECLEVEL( cerr, "ezgcd" );
563      if(F.mvar()==G.mvar())
564        d*=gcd(F,G);
565      else
566        return N (d);
567      return N (d);
568    }
569    else
570    {
571      g= content (G,G.mvar());
572      return N(d*gcd(F,g));
573    }
574  }
575  if ( G.isUnivariate())
576  {
577    f= content (F,F.mvar());
578    return N(d*gcd(G,f));
579  }
580
581  maxNumVars= tmax (getNumVars (F), getNumVars (G));
582  sizeF= size (F);
583  sizeG= size (G);
584
585  if (!isRat)
586    On (SW_RATIONAL);
587  if (sizeF/maxNumVars > 500 && sizeG/maxNumVars > 500)
588  {
589    Off(SW_USE_EZGCD);
590    result=gcd( F, G );
591    On(SW_USE_EZGCD);
592    if (!isRat)
593      Off (SW_RATIONAL);
594    result /= icontent (result);
595    DEBDECLEVEL( cerr, "ezgcd" );
596    return N (d*result);
597  }
598
599  int dummy= 0;
600  if ( gcd_test_one( F, G, false, dummy ) )
601  {
602    DEBDECLEVEL( cerr, "ezgcd" );
603    if (!isRat)
604      Off (SW_RATIONAL);
605    return N (d);
606  }
607  lcF = LC( F, x ); lcG = LC( G, x );
608  lcD = gcd( lcF, lcG );
609  delta = 0;
610  degF = degree( F, x ); degG = degree( G, x );
611  t = tmax( F.level(), G.level() );
612  if ( ! internal )
613    b = REvaluation( 2, t, IntRandom( 25 ) );
614  while ( ! gcdfound )
615  {
616    /// ---> A2
617    DEBOUTLN( cerr, "search for evaluation, delta = " << delta );
618    DEBOUTLN( cerr, "F = " << F );
619    DEBOUTLN( cerr, "G = " << G );
620    TIMING_START (ez_eval)
621    if (!findeval( F, G, Fb, Gb, Db, b, 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 EZ1: ")
634    DEBOUTLN( cerr, "found evaluation b = " << b );
635    DEBOUTLN( cerr, "F(b) = " << Fb );
636    DEBOUTLN( cerr, "G(b) = " << Gb );
637    DEBOUTLN( cerr, "D(b) = " << Db );
638    delta = degree( Db );
639    /// ---> A3
640    if (delta == degF)
641    {
642      if (degF <= degG  && fdivides (F, G))
643      {
644        DEBDECLEVEL( cerr, "ezgcd" );
645        if (!isRat)
646          Off (SW_RATIONAL);
647        return N (d*F);
648      }
649      else
650        delta--;
651    }
652    else if (delta == degG)
653    {
654      if (degG <= degF && fdivides( G, F ))
655      {
656        DEBDECLEVEL( cerr, "ezgcd" );
657        if (!isRat)
658          Off (SW_RATIONAL);
659        return N (d*G);
660      }
661      else
662        delta--;
663    }
664    if ( delta == 0 )
665    {
666      DEBDECLEVEL( cerr, "ezgcd" );
667      if (!isRat)
668        Off (SW_RATIONAL);
669      return N (d);
670    }
671    /// ---> A4
672    //deltaold = delta;
673    while ( 1 )
674    {
675      bt = b;
676      TIMING_START (ez_eval)
677      if (!findeval( F, G, Fbt, Gbt, Dbt, bt, delta, degF, degG, maxeval, count,
678                     o, 25,l ))
679      {
680        Off(SW_USE_EZGCD);
681        result=gcd( F, G );
682        On(SW_USE_EZGCD);
683        if (!isRat)
684          Off (SW_RATIONAL);
685        DEBDECLEVEL( cerr, "ezgcd" );
686        result /= icontent (result);
687        return N (d*result);
688      }
689      TIMING_END_AND_PRINT (ez_eval, "time to find eval point in EZ2: ")
690      int dd=degree( Dbt );
691      if ( dd /*degree( Dbt )*/ == 0 )
692      {
693        DEBDECLEVEL( cerr, "ezgcd" );
694        if (!isRat)
695          Off (SW_RATIONAL);
696        return N (d);
697      }
698      if ( dd /*degree( Dbt )*/ == delta )
699        break;
700      else  if ( dd /*degree( Dbt )*/ < delta )
701      {
702        delta = dd /*degree( Dbt )*/;
703        b = bt;
704        Db = Dbt; Fb = Fbt; Gb = Gbt;
705      }
706      DEBOUTLN( cerr, "now after A4, delta = " << delta );
707      /// ---> A5
708      if (delta == degF)
709      {
710        if (degF <= degG  && fdivides (F, G))
711        {
712          DEBDECLEVEL( cerr, "ezgcd" );
713          if (!isRat)
714            Off (SW_RATIONAL);
715          return N (d*F);
716        }
717        else
718          delta--;
719      }
720      else if (delta == degG)
721      {
722        if (degG <= degF && fdivides( G, F ))
723        {
724          DEBDECLEVEL( cerr, "ezgcd" );
725          if (!isRat)
726            Off (SW_RATIONAL);
727          return N (d*G);
728        }
729        else
730          delta--;
731      }
732      if ( delta == 0 )
733      {
734        DEBDECLEVEL( cerr, "ezgcd" );
735        if (!isRat)
736          Off (SW_RATIONAL);
737        return N (d);
738      }
739    }
740    if ( delta != degF && delta != degG )
741    {
742      /// ---> A6
743      bool B_is_F;
744      CanonicalForm xxx1, xxx2;
745      CanonicalForm buf;
746      DD[1] = Fb / Db;
747      buf= Gb/Db;
748      xxx1 = gcd( DD[1], Db );
749      xxx2 = gcd( buf, Db );
750      if (((xxx1.inCoeffDomain() && xxx2.inCoeffDomain()) &&
751          (size (F) <= size (G)))
752            || (xxx1.inCoeffDomain() && !xxx2.inCoeffDomain()))
753      {
754        B = F;
755        DD[2] = Db;
756        lcDD[1] = lcF;
757        lcDD[2] = lcD;
758        B_is_F = true;
759      }
760      else if (((xxx1.inCoeffDomain() && xxx2.inCoeffDomain()) &&
761                (size (G) < size (F)))
762                || (!xxx1.inCoeffDomain() && xxx2.inCoeffDomain()))
763      {
764        DD[1] = buf;
765        B = G;
766        DD[2] = Db;
767        lcDD[1] = lcG;
768        lcDD[2] = lcD;
769        B_is_F = false;
770      }
771      else
772      {
773        //special case
774        Off(SW_USE_EZGCD);
775        result=gcd( F, G );
776        On(SW_USE_EZGCD);
777        if (!isRat)
778          Off (SW_RATIONAL);
779        DEBDECLEVEL( cerr, "ezgcd" );
780        result /= icontent (result);
781        return N (d*result);
782      }
783      /// ---> A7
784      DD[2] = DD[2] * ( b( lcDD[2] ) / lc( DD[2] ) );
785      DD[1] = DD[1] * ( b( lcDD[1] ) / lc( DD[1] ) );
786      DEBOUTLN( cerr, "(hensel) B    = " << B );
787      DEBOUTLN( cerr, "(hensel) lcB  = " << LC( B, Variable(1) ) );
788      DEBOUTLN( cerr, "(hensel) b(B) = " << b(B) );
789      DEBOUTLN( cerr, "(hensel) DD   = " << DD );
790      DEBOUTLN( cerr, "(hensel) lcDD = " << lcDD );
791      TIMING_START (ez_hensel_lift)
792      gcdfound= Hensel (B*lcD, DD, b, lcDD);
793      TIMING_END_AND_PRINT (ez_hensel_lift, "time to hensel lift in EZ: ")
794      DEBOUTLN( cerr, "(hensel finished) DD   = " << DD );
795
796      if (gcdfound == -1)
797      {
798        Off (SW_USE_EZGCD);
799        result= gcd (F,G);
800        On (SW_USE_EZGCD);
801        if (!isRat)
802          Off (SW_RATIONAL);
803        DEBDECLEVEL( cerr, "ezgcd" );
804        result /= icontent (result);
805        return N (d*result);
806      }
807
808      if (gcdfound)
809      {
810        TIMING_START (ez_termination)
811        contcand= content (DD[2], Variable (1));
812        cand = DD[2] / contcand;
813        if (B_is_F)
814          gcdfound = fdivides( cand, G ) && cand*(DD[1]/(lcD/contcand)) == F;
815        else
816          gcdfound = fdivides( cand, F ) && cand*(DD[1]/(lcD/contcand)) == G;
817        TIMING_END_AND_PRINT (ez_termination,
818                              "time for termination test in EZ: ")
819      }
820      /// ---> A8 (gcdfound)
821    }
822    delta--;
823  }
824  /// ---> A9
825  DEBDECLEVEL( cerr, "ezgcd" );
826  cand *= bCommonDen (cand);
827  if (!isRat)
828    Off (SW_RATIONAL);
829  cand /= icontent (cand);
830  return N (d*cand);
831}
832#endif
833
834/// Extended Zassenhaus GCD over Z.
835/// In case things become too dense we switch to a modular algorithm.
836CanonicalForm
837ezgcd ( const CanonicalForm & FF, const CanonicalForm & GG )
838{
839#ifdef HAVE_NTL
840  REvaluation b;
841  return ezgcd( FF, GG, b, false );
842#else
843  Off (SW_USE_EZGCD);
844  return gcd (FF, GG);
845  On (SW_USE_EZGCD);
846#endif
847}
848
849#ifdef HAVE_NTL
850// parameters for heuristic
851static int maxNumEval= 200;
852static int sizePerVars1= 500; //try dense gcd if size/#variables is bigger
853
854/// Extended Zassenhaus GCD for finite fields.
855/// In case things become too dense we switch to a modular algorithm.
856CanonicalForm EZGCD_P( const CanonicalForm & FF, const CanonicalForm & GG )
857{
858  if (FF.isZero() && degree(GG) > 0) return GG/Lc(GG);
859  else if (GG.isZero() && degree (FF) > 0) return FF/Lc(FF);
860  else if (FF.isZero() && GG.isZero()) return FF.genOne();
861  if (FF.inBaseDomain() || GG.inBaseDomain()) return FF.genOne();
862  if (FF.isUnivariate() && fdivides(FF, GG)) return FF/Lc(FF);
863  if (GG.isUnivariate() && fdivides(GG, FF)) return GG/Lc(GG);
864  if (FF == GG) return FF/Lc(FF);
865
866  int maxNumVars= tmax (getNumVars (FF), getNumVars (GG));
867  Variable a, oldA;
868  int sizeF= size (FF);
869  int sizeG= size (GG);
870
871  if (sizeF==1)
872  {
873    return gcd_mon( FF, GG );
874  }
875  else if (sizeG==1)
876  {
877    return gcd_mon( GG, FF );
878  }
879
880  if (sizeF/maxNumVars > sizePerVars1 && sizeG/maxNumVars > sizePerVars1)
881  {
882    if (hasFirstAlgVar (FF, a) || hasFirstAlgVar (GG, a))
883      return modGCDFq (FF, GG, a);
884    else if (CFFactory::gettype() == GaloisFieldDomain)
885      return modGCDGF (FF, GG);
886    else
887      return modGCDFp (FF, GG);
888  }
889
890  CanonicalForm F, G, f, g, d, Fb, Gb, Db, Fbt, Gbt, Dbt, B0, B, D0, lcF, lcG,
891                lcD;
892  CFArray DD( 1, 2 ), lcDD( 1, 2 );
893  int degF, degG, delta, count;
894  int maxeval;
895  maxeval= tmin((getCharacteristic()/
896                (int)(ilog2(getCharacteristic())*log2exp))*2, maxNumEval);
897  count= 0; // number of eval. used
898  REvaluation b, bt;
899  int gcdfound = 0;
900  Variable x = Variable(1);
901
902  F= FF;
903  G= GG;
904
905  CFMap M,N;
906  int smallestDegLev;
907  TIMING_DEFINE(ez_p_compress);
908  TIMING_START (ez_p_compress);
909  int best_level= compress4EZGCD (F, G, M, N, smallestDegLev);
910
911  if (best_level == 0) return G.genOne();
912
913  F= M (F);
914  G= M (G);
915  TIMING_END_AND_PRINT (ez_p_compress, "time for compression in EZ_P: ")
916
917  TIMING_DEFINE (ez_p_content)
918  TIMING_START (ez_p_content)
919  f = content( F, x ); g = content( G, x );
920  d = gcd( f, g );
921  F /= f; G /= g;
922  TIMING_END_AND_PRINT (ez_p_content, "time to extract content in EZ_P: ")
923
924  if( F.isUnivariate() && G.isUnivariate() )
925  {
926    if( F.mvar() == G.mvar() )
927      d *= gcd( F, G );
928    else
929      return N (d);
930    return N (d);
931  }
932  if ( F.isUnivariate())
933  {
934    g= content (G,G.mvar());
935    return N(d*gcd(F,g));
936  }
937  if ( G.isUnivariate())
938  {
939    f= content (F,F.mvar());
940    return N(d*gcd(G,f));
941  }
942
943  maxNumVars= tmax (getNumVars (F), getNumVars (G));
944  sizeF= size (F);
945  sizeG= size (G);
946
947  if (sizeF/maxNumVars > sizePerVars1 && sizeG/maxNumVars > sizePerVars1)
948  {
949    if (hasFirstAlgVar (F, a) || hasFirstAlgVar (G, a))
950      return N (d*modGCDFq (F, G, a));
951    else if (CFFactory::gettype() == GaloisFieldDomain)
952      return N (d*modGCDGF (F, G));
953    else
954      return N (d*modGCDFp (F, G));
955  }
956
957  int dummy= 0;
958  if( gcd_test_one( F, G, false, dummy ) )
959  {
960    return N (d);
961  }
962
963  bool passToGF= false;
964  bool extOfExt= false;
965  int p= getCharacteristic();
966  bool algExtension= (hasFirstAlgVar(F,a) || hasFirstAlgVar(G,a));
967  int k= 1;
968  CanonicalForm primElem, imPrimElem;
969  CFList source, dest;
970  if (p < 50 && CFFactory::gettype() != GaloisFieldDomain && !algExtension)
971  {
972    if (p == 2)
973      setCharacteristic (2, 12, 'Z');
974    else if (p == 3)
975      setCharacteristic (3, 4, 'Z');
976    else if (p == 5 || p == 7)
977      setCharacteristic (p, 3, 'Z');
978    else
979      setCharacteristic (p, 2, 'Z');
980    passToGF= true;
981    F= F.mapinto();
982    G= G.mapinto();
983    maxeval= 2*ipower (p, getGFDegree());
984  }
985  else if (CFFactory::gettype() == GaloisFieldDomain &&
986           ipower (p , getGFDegree()) < 50)
987  {
988    k= getGFDegree();
989    if (ipower (p, 2*k) > 50)
990      setCharacteristic (p, 2*k, gf_name);
991    else
992      setCharacteristic (p, 3*k, gf_name);
993    F= GFMapUp (F, k);
994    G= GFMapUp (G, k);
995    maxeval= tmin (2*ipower (p, getGFDegree()), maxNumEval);
996  }
997  else if (p < 50 && algExtension && CFFactory::gettype() != GaloisFieldDomain)
998  {
999    int d= degree (getMipo (a));
1000    oldA= a;
1001    Variable v2;
1002    if (p == 2 && d < 6)
1003    {
1004      if (fac_NTL_char != p)
1005      {
1006        fac_NTL_char= p;
1007        zz_p::init (p);
1008      }
1009      bool primFail= false;
1010      Variable vBuf;
1011      primElem= primitiveElement (a, vBuf, primFail);
1012      ASSERT (!primFail, "failure in integer factorizer");
1013      if (d < 3)
1014      {
1015        zz_pX NTLIrredpoly;
1016        BuildIrred (NTLIrredpoly, d*3);
1017        CanonicalForm newMipo= convertNTLzzpX2CF (NTLIrredpoly, Variable (1));
1018        v2= rootOf (newMipo);
1019      }
1020      else
1021      {
1022        zz_pX NTLIrredpoly;
1023        BuildIrred (NTLIrredpoly, d*2);
1024        CanonicalForm newMipo= convertNTLzzpX2CF (NTLIrredpoly, Variable (1));
1025        v2= rootOf (newMipo);
1026      }
1027      imPrimElem= mapPrimElem (primElem, a, v2);
1028      extOfExt= true;
1029    }
1030    else if ((p == 3 && d < 4) || ((p == 5 || p == 7) && d < 3))
1031    {
1032      if (fac_NTL_char != p)
1033      {
1034        fac_NTL_char= p;
1035        zz_p::init (p);
1036      }
1037      bool primFail= false;
1038      Variable vBuf;
1039      primElem= primitiveElement (a, vBuf, primFail);
1040      ASSERT (!primFail, "failure in integer factorizer");
1041      zz_pX NTLIrredpoly;
1042      BuildIrred (NTLIrredpoly, d*2);
1043      CanonicalForm newMipo= convertNTLzzpX2CF (NTLIrredpoly, Variable (1));
1044      v2= rootOf (newMipo);
1045      imPrimElem= mapPrimElem (primElem, a, v2);
1046      extOfExt= true;
1047    }
1048    if (extOfExt)
1049    {
1050      maxeval= tmin (2*ipower (p, degree (getMipo (v2))), maxNumEval);
1051      F= mapUp (F, a, v2, primElem, imPrimElem, source, dest);
1052      G= mapUp (G, a, v2, primElem, imPrimElem, source, dest);
1053      a= v2;
1054    }
1055  }
1056
1057  lcF = LC( F, x ); lcG = LC( G, x );
1058  lcD = gcd( lcF, lcG );
1059
1060  delta = 0;
1061  degF = degree( F, x ); degG = degree( G, x );
1062
1063  if (algExtension)
1064    b = REvaluation( 2, tmax(F.level(), G.level()), AlgExtRandomF( a ) );
1065  else
1066  { // both not in extension given by algebraic variable
1067    if (CFFactory::gettype() != GaloisFieldDomain)
1068      b = REvaluation( 2, tmax(F.level(), G.level()), FFRandom() );
1069    else
1070      b = REvaluation( 2, tmax(F.level(), G.level()), GFRandom() );
1071  }
1072
1073  CanonicalForm cand, contcand;
1074  CanonicalForm result;
1075  int o, t;
1076  o= 0;
1077  t= 1;
1078  int goodPointCount= 0;
1079  TIMING_DEFINE(ez_p_eval);
1080  while( !gcdfound )
1081  {
1082    TIMING_START (ez_p_eval);
1083    if( !findeval( F, G, Fb, Gb, Db, b, delta, degF, degG, maxeval, count, o,
1084         maxeval/maxNumVars, t ))
1085    { // too many eval. used --> try another method
1086      Off (SW_USE_EZGCD_P);
1087      result= gcd (F,G);
1088      On (SW_USE_EZGCD_P);
1089      if (passToGF)
1090      {
1091        CanonicalForm mipo= gf_mipo;
1092        setCharacteristic (p);
1093        Variable alpha= rootOf (mipo.mapinto());
1094        result= GF2FalphaRep (result, alpha);
1095        prune (alpha);
1096      }
1097      if (k > 1)
1098      {
1099        result= GFMapDown (result, k);
1100        setCharacteristic (p, k, gf_name);
1101      }
1102      if (extOfExt)
1103      {
1104        result= mapDown (result, primElem, imPrimElem, oldA, dest, source);
1105        prune1 (oldA);
1106      }
1107      return N (d*result);
1108    }
1109    TIMING_END_AND_PRINT (ez_p_eval, "time for eval point search in EZ_P1: ");
1110    delta = degree( Db );
1111    if (delta == degF)
1112    {
1113      if (degF <= degG && fdivides (F, G))
1114      {
1115        if (passToGF)
1116        {
1117          CanonicalForm mipo= gf_mipo;
1118          setCharacteristic (p);
1119          Variable alpha= rootOf (mipo.mapinto());
1120          F= GF2FalphaRep (F, alpha);
1121          prune (alpha);
1122        }
1123        if (k > 1)
1124        {
1125          F= GFMapDown (F, k);
1126          setCharacteristic (p, k, gf_name);
1127        }
1128        if (extOfExt)
1129        {
1130          F= mapDown (F, primElem, imPrimElem, oldA, dest, source);
1131          prune1 (oldA);
1132        }
1133        return N (d*F);
1134      }
1135      else
1136        delta--;
1137    }
1138    else if (delta == degG)
1139    {
1140      if (degG <= degF && fdivides (G, F))
1141      {
1142        if (passToGF)
1143        {
1144          CanonicalForm mipo= gf_mipo;
1145          setCharacteristic (p);
1146          Variable alpha= rootOf (mipo.mapinto());
1147          G= GF2FalphaRep (G, alpha);
1148          prune (alpha);
1149        }
1150        if (k > 1)
1151        {
1152          G= GFMapDown (G, k);
1153          setCharacteristic (p, k, gf_name);
1154        }
1155        if (extOfExt)
1156        {
1157          G= mapDown (G, primElem, imPrimElem, oldA, dest, source);
1158          prune1 (oldA);
1159        }
1160        return N (d*G);
1161      }
1162      else
1163        delta--;
1164    }
1165    if( delta == 0 )
1166    {
1167      if (passToGF)
1168        setCharacteristic (p);
1169      if (k > 1)
1170        setCharacteristic (p, k, gf_name);
1171      return N (d);
1172    }
1173    while( true )
1174    {
1175      bt = b;
1176      TIMING_START (ez_p_eval);
1177      if( !findeval(F,G,Fbt,Gbt,Dbt, bt, delta, degF, degG, maxeval, count, o,
1178           maxeval/maxNumVars, t ))
1179      { // too many eval. used --> try another method
1180        Off (SW_USE_EZGCD_P);
1181        result= gcd (F,G);
1182        On (SW_USE_EZGCD_P);
1183        if (passToGF)
1184        {
1185          CanonicalForm mipo= gf_mipo;
1186          setCharacteristic (p);
1187          Variable alpha= rootOf (mipo.mapinto());
1188          result= GF2FalphaRep (result, alpha);
1189          prune (alpha);
1190        }
1191        if (k > 1)
1192        {
1193          result= GFMapDown (result, k);
1194          setCharacteristic (p, k, gf_name);
1195        }
1196        if (extOfExt)
1197        {
1198          result= mapDown (result, primElem, imPrimElem, oldA, dest, source);
1199          prune1 (oldA);
1200        }
1201        return N (d*result);
1202      }
1203      TIMING_END_AND_PRINT (ez_p_eval, "time for eval point search in EZ_P2: ");
1204      int dd = degree( Dbt );
1205      if( dd == 0 )
1206      {
1207        if (passToGF)
1208          setCharacteristic (p);
1209        if (k > 1)
1210          setCharacteristic (p, k, gf_name);
1211        return N (d);
1212      }
1213      if( dd == delta )
1214      {
1215        goodPointCount++;
1216        if (goodPointCount == 5)
1217          break;
1218      }
1219      if( dd < delta )
1220      {
1221        goodPointCount= 0;
1222        delta = dd;
1223        b = bt;
1224        Db = Dbt; Fb = Fbt; Gb = Gbt;
1225      }
1226      if (delta == degF)
1227      {
1228        if (degF <= degG && fdivides (F, G))
1229        {
1230          if (passToGF)
1231          {
1232            CanonicalForm mipo= gf_mipo;
1233            setCharacteristic (p);
1234            Variable alpha= rootOf (mipo.mapinto());
1235            F= GF2FalphaRep (F, alpha);
1236            prune (alpha);
1237          }
1238          if (k > 1)
1239          {
1240            F= GFMapDown (F, k);
1241            setCharacteristic (p, k, gf_name);
1242          }
1243          if (extOfExt)
1244          {
1245            F= mapDown (F, primElem, imPrimElem, oldA, dest, source);
1246            prune1 (oldA);
1247          }
1248          return N (d*F);
1249        }
1250        else
1251          delta--;
1252      }
1253      else if (delta == degG)
1254      {
1255        if (degG <= degF && fdivides (G, F))
1256        {
1257          if (passToGF)
1258          {
1259            CanonicalForm mipo= gf_mipo;
1260            setCharacteristic (p);
1261            Variable alpha= rootOf (mipo.mapinto());
1262            G= GF2FalphaRep (G, alpha);
1263            prune (alpha);
1264          }
1265          if (k > 1)
1266          {
1267            G= GFMapDown (G, k);
1268            setCharacteristic (p, k, gf_name);
1269          }
1270          if (extOfExt)
1271          {
1272            G= mapDown (G, primElem, imPrimElem, oldA, dest, source);
1273            prune1 (oldA);
1274          }
1275          return N (d*G);
1276        }
1277        else
1278          delta--;
1279      }
1280      if( delta == 0 )
1281      {
1282        if (passToGF)
1283          setCharacteristic (p);
1284        if (k > 1)
1285          setCharacteristic (p, k, gf_name);
1286        return N (d);
1287      }
1288    }
1289    if( delta != degF && delta != degG )
1290    {
1291      bool B_is_F;
1292      CanonicalForm xxx1, xxx2;
1293      CanonicalForm buf;
1294      DD[1] = Fb / Db;
1295      buf= Gb/Db;
1296      xxx1 = gcd( DD[1], Db );
1297      xxx2 = gcd( buf, Db );
1298      if (((xxx1.inCoeffDomain() && xxx2.inCoeffDomain()) &&
1299          (size (F) <= size (G)))
1300          || (xxx1.inCoeffDomain() && !xxx2.inCoeffDomain()))
1301      {
1302        B = F;
1303        DD[2] = Db;
1304        lcDD[1] = lcF;
1305        lcDD[2] = lcD;
1306        B_is_F = true;
1307      }
1308      else if (((xxx1.inCoeffDomain() && xxx2.inCoeffDomain()) &&
1309               (size (G) < size (F)))
1310               || (!xxx1.inCoeffDomain() && xxx2.inCoeffDomain()))
1311      {
1312        DD[1] = buf;
1313        B = G;
1314        DD[2] = Db;
1315        lcDD[1] = lcG;
1316        lcDD[2] = lcD;
1317        B_is_F = false;
1318      }
1319      else // special case handling
1320      {
1321        Off (SW_USE_EZGCD_P);
1322        result= gcd (F,G);
1323        On (SW_USE_EZGCD_P);
1324        if (passToGF)
1325        {
1326          CanonicalForm mipo= gf_mipo;
1327          setCharacteristic (p);
1328          Variable alpha= rootOf (mipo.mapinto());
1329          result= GF2FalphaRep (result, alpha);
1330          prune (alpha);
1331        }
1332        if (k > 1)
1333        {
1334          result= GFMapDown (result, k);
1335          setCharacteristic (p, k, gf_name);
1336        }
1337        if (extOfExt)
1338        {
1339          result= mapDown (result, primElem, imPrimElem, oldA, dest, source);
1340          prune1 (oldA);
1341        }
1342        return N (d*result);
1343      }
1344      DD[2] = DD[2] * ( b( lcDD[2] ) / lc( DD[2] ) );
1345      DD[1] = DD[1] * ( b( lcDD[1] ) / lc( DD[1] ) );
1346
1347      if (size (B*lcDD[2])/maxNumVars > sizePerVars1)
1348      {
1349        if (algExtension)
1350        {
1351          result= modGCDFq (F, G, a);
1352          if (extOfExt)
1353          {
1354            result= mapDown (result, primElem, imPrimElem, oldA, dest, source);
1355            prune1 (oldA);
1356          }
1357          return N (d*result);
1358        }
1359        if (CFFactory::gettype() == GaloisFieldDomain)
1360        {
1361          result= modGCDGF (F, G);
1362          if (passToGF)
1363          {
1364            CanonicalForm mipo= gf_mipo;
1365            setCharacteristic (p);
1366            Variable alpha= rootOf (mipo.mapinto());
1367            result= GF2FalphaRep (result, alpha);
1368            prune (alpha);
1369          }
1370          if (k > 1)
1371          {
1372            result= GFMapDown (result, k);
1373            setCharacteristic (p, k, gf_name);
1374          }
1375          return N (d*result);
1376        }
1377        else
1378          return N (d*modGCDFp (F,G));
1379      }
1380
1381      TIMING_DEFINE(ez_p_hensel_lift);
1382      TIMING_START (ez_p_hensel_lift);
1383      gcdfound= Hensel (B*lcD, DD, b, lcDD);
1384      TIMING_END_AND_PRINT (ez_p_hensel_lift, "time for Hensel lift in EZ_P: ");
1385
1386      if (gcdfound == -1) //things became dense
1387      {
1388        if (algExtension)
1389        {
1390          result= modGCDFq (F, G, a);
1391          if (extOfExt)
1392          {
1393            result= mapDown (result, primElem, imPrimElem, oldA, dest, source);
1394            prune1 (oldA);
1395          }
1396          return N (d*result);
1397        }
1398        if (CFFactory::gettype() == GaloisFieldDomain)
1399        {
1400          result= modGCDGF (F, G);
1401          if (passToGF)
1402          {
1403            CanonicalForm mipo= gf_mipo;
1404            setCharacteristic (p);
1405            Variable alpha= rootOf (mipo.mapinto());
1406            result= GF2FalphaRep (result, alpha);
1407            prune (alpha);
1408          }
1409          if (k > 1)
1410          {
1411            result= GFMapDown (result, k);
1412            setCharacteristic (p, k, gf_name);
1413          }
1414          return N (d*result);
1415        }
1416        else
1417        {
1418          if (p >= cf_getBigPrime(0))
1419            return N (d*sparseGCDFp (F,G));
1420          else
1421            return N (d*modGCDFp (F,G));
1422        }
1423      }
1424
1425      if (gcdfound == 1)
1426      {
1427        TIMING_DEFINE(termination_test);
1428        TIMING_START (termination_test);
1429        contcand= content (DD[2], Variable (1));
1430        cand = DD[2] / contcand;
1431        if (B_is_F)
1432          gcdfound = fdivides( cand, G ) && cand*(DD[1]/(lcD/contcand)) == F;
1433        else
1434          gcdfound = fdivides( cand, F ) && cand*(DD[1]/(lcD/contcand)) == G;
1435        TIMING_END_AND_PRINT (termination_test,
1436                              "time for termination test EZ_P: ");
1437
1438        if (passToGF && gcdfound)
1439        {
1440          CanonicalForm mipo= gf_mipo;
1441          setCharacteristic (p);
1442          Variable alpha= rootOf (mipo.mapinto());
1443          cand= GF2FalphaRep (cand, alpha);
1444          prune (alpha);
1445        }
1446        if (k > 1 && gcdfound)
1447        {
1448          cand= GFMapDown (cand, k);
1449          setCharacteristic (p, k, gf_name);
1450        }
1451        if (extOfExt && gcdfound)
1452        {
1453          cand= mapDown (cand, primElem, imPrimElem, oldA, dest, source);
1454          prune1 (oldA);
1455        }
1456      }
1457    }
1458    delta--;
1459    goodPointCount= 0;
1460  }
1461  return N (d*cand);
1462}
1463#endif
1464
Note: See TracBrowser for help on using the repository browser.