source: git/factory/cfEzgcd.cc @ afa00e

fieker-DuValspielwiese
Last change on this file since afa00e was 53e6c3, checked in by Martin Lee <martinlee84@…>, 10 years ago
deleted some redundant code
  • Property mode set to 100644
File size: 34.2 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        ((getCharacteristic() > degree (U,i)) || getCharacteristic() == 0))
311    {
312      termEstimate *= degree (U,i)*2;
313      termEstimate /= 3;
314    }
315    evaluation.append (A [i]);
316  }
317  if (termEstimate/getNumVars(U) > 500)
318    return -1;
319  CFList UEval;
320  CanonicalForm shiftedU= myShift2Zero (U, UEval, evaluation);
321
322  if (size (shiftedU)/getNumVars (U) > 500)
323    return -1;
324
325  CFArray shiftedLCs= CFArray (2);
326  CFList shiftedLCsEval1, shiftedLCsEval2;
327  shiftedLCs[0]= myShift2Zero (LCs[1], shiftedLCsEval1, evaluation);
328  shiftedLCs[1]= myShift2Zero (LCs[2], shiftedLCsEval2, evaluation);
329  factors.insert (1);
330  int liftBound= degree (UEval.getLast(), 2) + 1;
331  CFArray Pi;
332  CFMatrix M= CFMatrix (liftBound, factors.length() - 1);
333  CFList diophant;
334  CFArray lcs= CFArray (2);
335  lcs [0]= shiftedLCsEval1.getFirst();
336  lcs [1]= shiftedLCsEval2.getFirst();
337  nonMonicHenselLift12 (UEval.getFirst(), factors, liftBound, Pi, diophant, M,
338                        lcs, false);
339
340  for (CFListIterator i= factors; i.hasItem(); i++)
341  {
342    if (!fdivides (i.getItem(), UEval.getFirst()))
343      return 0;
344  }
345
346  int * liftBounds;
347  bool noOneToOne= false;
348  if (U.level() > 2)
349  {
350    liftBounds= new int [U.level() - 1]; /* index: 0.. U.level()-2 */
351    liftBounds[0]= liftBound;
352    for (int i= 1; i < U.level() - 1; i++)
353      liftBounds[i]= degree (shiftedU, Variable (i + 2)) + 1;
354    factors= nonMonicHenselLift2 (UEval, factors, liftBounds, U.level() - 1,
355                                  false, shiftedLCsEval1, shiftedLCsEval2, Pi,
356                                  diophant, noOneToOne);
357    delete [] liftBounds;
358    if (noOneToOne)
359      return 0;
360  }
361  G[1]= factors.getFirst();
362  G[2]= factors.getLast();
363  G[1]= myReverseShift (G[1], evaluation);
364  G[2]= myReverseShift (G[2], evaluation);
365  G[1]= NN (G[1]);
366  G[2]= NN (G[2]);
367  return 1;
368}
369
370static
371bool findeval (const CanonicalForm & F, const CanonicalForm & G,
372               CanonicalForm & Fb, CanonicalForm & Gb, CanonicalForm & Db,
373               REvaluation & b, int delta, int degF, int degG, int maxeval,
374               int & count, int& k, int bound, int& l)
375{
376  if( count == 0 && delta != 0)
377  {
378    if( count++ > maxeval )
379      return false;
380  }
381  if (count > 0)
382  {
383    b.nextpoint(k);
384    if (k == 0)
385      k++;
386    l++;
387    if (l > bound)
388    {
389      l= 1;
390      k++;
391      if (k > tmax (F.level(), G.level()) - 1)
392        return false;
393      b.nextpoint (k);
394    }
395    if (count++ > maxeval)
396      return false;
397  }
398  while( true )
399  {
400    Fb = b( F );
401    if( degree( Fb, 1 ) == degF )
402    {
403      Gb = b( G );
404      if( degree( Gb, 1 ) == degG )
405      {
406        Db = gcd( Fb, Gb );
407        if( delta > 0 )
408        {
409          if( degree( Db, 1 ) <= delta )
410            return true;
411        }
412        else
413        {
414          k++;
415          return true;
416        }
417      }
418    }
419    if (k == 0)
420      k++;
421    b.nextpoint(k);
422    l++;
423    if (l > bound)
424    {
425      l= 1;
426      k++;
427      if (k > tmax (F.level(), G.level()) - 1)
428        return false;
429      b.nextpoint (k);
430    }
431    if( count++ > maxeval )
432      return false;
433  }
434}
435
436/// real implementation of EZGCD over Z
437static CanonicalForm
438ezgcd ( const CanonicalForm & FF, const CanonicalForm & GG, REvaluation & b,
439        bool internal )
440{
441  bool isRat= isOn (SW_RATIONAL);
442
443  int maxNumVars= tmax (getNumVars (FF), getNumVars (GG));
444  int sizeF= size (FF);
445  int sizeG= size (GG);
446
447
448  if (!isRat)
449    On (SW_RATIONAL);
450  if (sizeF/maxNumVars > 500 && sizeG/maxNumVars > 500)
451  {
452    Off(SW_USE_EZGCD);
453    CanonicalForm result=gcd( FF, GG );
454    On(SW_USE_EZGCD);
455    if (!isRat)
456      Off (SW_RATIONAL);
457    result /= icontent (result);
458    DEBDECLEVEL( cerr, "ezgcd" );
459    return result;
460  }
461
462
463  CanonicalForm F, G, f, g, d, Fb, Gb, Db, Fbt, Gbt, Dbt, B0, B, D0, lcF, lcG,
464                lcD, cand, contcand, result;
465  CFArray DD( 1, 2 ), lcDD( 1, 2 );
466  int degF, degG, delta, t, count, maxeval;
467  REvaluation bt;
468  int gcdfound = 0;
469  Variable x = Variable(1);
470  count= 0;
471  maxeval= 200;
472  int o, l;
473  o= 0;
474  l= 1;
475
476  if (!isRat)
477    On (SW_RATIONAL);
478  F= FF*bCommonDen (FF);
479  G= GG*bCommonDen (GG);
480  if (!isRat)
481    Off (SW_RATIONAL);
482
483  TIMING_START (ez_compress)
484  CFMap M,N;
485  int smallestDegLev;
486  int best_level= compress4EZGCD (F, G, M, N, smallestDegLev);
487
488  if (best_level == 0)
489  {
490    DEBDECLEVEL( cerr, "ezgcd" );
491    return G.genOne();
492  }
493
494  F= M (F);
495  G= M (G);
496  TIMING_END_AND_PRINT (ez_compress, "time for compression in EZ: ")
497
498  DEBINCLEVEL( cerr, "ezgcd" );
499  DEBOUTLN( cerr, "FF = " << FF );
500  DEBOUTLN( cerr, "GG = " << GG );
501  TIMING_START (ez_content)
502  f = content( F, x ); g = content( G, x ); d = gcd( f, g );
503  DEBOUTLN( cerr, "f = " << f );
504  DEBOUTLN( cerr, "g = " << g );
505  F /= f; G /= g;
506  TIMING_END_AND_PRINT (ez_content, "time to extract content in EZ: ")
507  if ( F.isUnivariate() && G.isUnivariate() )
508  {
509    DEBDECLEVEL( cerr, "ezgcd" );
510    if(F.mvar()==G.mvar())
511      d*=gcd(F,G);
512    else
513      return N (d);
514    return N (d);
515  }
516  if ( F.isUnivariate())
517  {
518    g= content (G,G.mvar());
519    return N(d*gcd(F,g));
520  }
521  if ( G.isUnivariate())
522  {
523    f= content (F,F.mvar());
524    return N(d*gcd(G,f));
525  }
526
527  maxNumVars= tmax (getNumVars (F), getNumVars (G));
528  sizeF= size (F);
529  sizeG= size (G);
530
531  if (!isRat)
532    On (SW_RATIONAL);
533  if (sizeF/maxNumVars > 500 && sizeG/maxNumVars > 500)
534  {
535    Off(SW_USE_EZGCD);
536    result=gcd( F, G );
537    On(SW_USE_EZGCD);
538    if (!isRat)
539      Off (SW_RATIONAL);
540    result /= icontent (result);
541    DEBDECLEVEL( cerr, "ezgcd" );
542    return N (d*result);
543  }
544
545  int dummy= 0;
546  if ( gcd_test_one( F, G, false, dummy ) )
547  {
548    DEBDECLEVEL( cerr, "ezgcd" );
549    if (!isRat)
550      Off (SW_RATIONAL);
551    return N (d);
552  }
553  lcF = LC( F, x ); lcG = LC( G, x );
554  lcD = gcd( lcF, lcG );
555  delta = 0;
556  degF = degree( F, x ); degG = degree( G, x );
557  t = tmax( F.level(), G.level() );
558  if ( ! internal )
559    b = REvaluation( 2, t, IntRandom( 25 ) );
560  while ( ! gcdfound )
561  {
562    /// ---> A2
563    DEBOUTLN( cerr, "search for evaluation, delta = " << delta );
564    DEBOUTLN( cerr, "F = " << F );
565    DEBOUTLN( cerr, "G = " << G );
566    TIMING_START (ez_eval)
567    if (!findeval( F, G, Fb, Gb, Db, b, delta, degF, degG, maxeval, count,
568                   o, 25, l))
569    {
570      Off(SW_USE_EZGCD);
571      result=gcd( F, G );
572      On(SW_USE_EZGCD);
573      if (!isRat)
574        Off (SW_RATIONAL);
575      DEBDECLEVEL( cerr, "ezgcd" );
576      result /= icontent (result);
577      return N (d*result);
578    }
579    TIMING_END_AND_PRINT (ez_eval, "time to find eval point in EZ1: ")
580    DEBOUTLN( cerr, "found evaluation b = " << b );
581    DEBOUTLN( cerr, "F(b) = " << Fb );
582    DEBOUTLN( cerr, "G(b) = " << Gb );
583    DEBOUTLN( cerr, "D(b) = " << Db );
584    delta = degree( Db );
585    /// ---> A3
586    if (delta == degF)
587    {
588      if (degF <= degG  && fdivides (F, G))
589      {
590        DEBDECLEVEL( cerr, "ezgcd" );
591        if (!isRat)
592          Off (SW_RATIONAL);
593        return N (d*F);
594      }
595      else
596        delta--;
597    }
598    else if (delta == degG)
599    {
600      if (degG <= degF && fdivides( G, F ))
601      {
602        DEBDECLEVEL( cerr, "ezgcd" );
603        if (!isRat)
604          Off (SW_RATIONAL);
605        return N (d*G);
606      }
607      else
608        delta--;
609    }
610    if ( delta == 0 )
611    {
612      DEBDECLEVEL( cerr, "ezgcd" );
613      if (!isRat)
614        Off (SW_RATIONAL);
615      return N (d);
616    }
617    /// ---> A4
618    //deltaold = delta;
619    while ( 1 )
620    {
621      bt = b;
622      TIMING_START (ez_eval)
623      if (!findeval( F, G, Fbt, Gbt, Dbt, bt, delta, degF, degG, maxeval, count,
624                     o, 25,l ))
625      {
626        Off(SW_USE_EZGCD);
627        result=gcd( F, G );
628        On(SW_USE_EZGCD);
629        if (!isRat)
630          Off (SW_RATIONAL);
631        DEBDECLEVEL( cerr, "ezgcd" );
632        result /= icontent (result);
633        return N (d*result);
634      }
635      TIMING_END_AND_PRINT (ez_eval, "time to find eval point in EZ2: ")
636      int dd=degree( Dbt );
637      if ( dd /*degree( Dbt )*/ == 0 )
638      {
639        DEBDECLEVEL( cerr, "ezgcd" );
640        if (!isRat)
641          Off (SW_RATIONAL);
642        return N (d);
643      }
644      if ( dd /*degree( Dbt )*/ == delta )
645        break;
646      else  if ( dd /*degree( Dbt )*/ < delta )
647      {
648        delta = dd /*degree( Dbt )*/;
649        b = bt;
650        Db = Dbt; Fb = Fbt; Gb = Gbt;
651      }
652      DEBOUTLN( cerr, "now after A4, delta = " << delta );
653      /// ---> A5
654      if (delta == degF)
655      {
656        if (degF <= degG  && fdivides (F, G))
657        {
658          DEBDECLEVEL( cerr, "ezgcd" );
659          if (!isRat)
660            Off (SW_RATIONAL);
661          return N (d*F);
662        }
663        else
664          delta--;
665      }
666      else if (delta == degG)
667      {
668        if (degG <= degF && fdivides( G, F ))
669        {
670          DEBDECLEVEL( cerr, "ezgcd" );
671          if (!isRat)
672            Off (SW_RATIONAL);
673          return N (d*G);
674        }
675        else
676          delta--;
677      }
678      if ( delta == 0 )
679      {
680        DEBDECLEVEL( cerr, "ezgcd" );
681        if (!isRat)
682          Off (SW_RATIONAL);
683        return N (d);
684      }
685    }
686    if ( delta != degF && delta != degG )
687    {
688      /// ---> A6
689      bool B_is_F;
690      CanonicalForm xxx1, xxx2;
691      CanonicalForm buf;
692      DD[1] = Fb / Db;
693      buf= Gb/Db;
694      xxx1 = gcd( DD[1], Db );
695      xxx2 = gcd( buf, Db );
696      if (((xxx1.inCoeffDomain() && xxx2.inCoeffDomain()) &&
697          (size (F) <= size (G)))
698            || (xxx1.inCoeffDomain() && !xxx2.inCoeffDomain()))
699      {
700        B = F;
701        DD[2] = Db;
702        lcDD[1] = lcF;
703        lcDD[2] = lcD;
704        B_is_F = true;
705      }
706      else if (((xxx1.inCoeffDomain() && xxx2.inCoeffDomain()) &&
707                (size (G) < size (F)))
708                || (!xxx1.inCoeffDomain() && xxx2.inCoeffDomain()))
709      {
710        DD[1] = buf;
711        B = G;
712        DD[2] = Db;
713        lcDD[1] = lcG;
714        lcDD[2] = lcD;
715        B_is_F = false;
716      }
717      else
718      {
719        //special case
720        Off(SW_USE_EZGCD);
721        result=gcd( F, G );
722        On(SW_USE_EZGCD);
723        if (!isRat)
724          Off (SW_RATIONAL);
725        DEBDECLEVEL( cerr, "ezgcd" );
726        result /= icontent (result);
727        return N (d*result);
728      }
729      /// ---> A7
730      DD[2] = DD[2] * ( b( lcDD[2] ) / lc( DD[2] ) );
731      DD[1] = DD[1] * ( b( lcDD[1] ) / lc( DD[1] ) );
732      DEBOUTLN( cerr, "(hensel) B    = " << B );
733      DEBOUTLN( cerr, "(hensel) lcB  = " << LC( B, Variable(1) ) );
734      DEBOUTLN( cerr, "(hensel) b(B) = " << b(B) );
735      DEBOUTLN( cerr, "(hensel) DD   = " << DD );
736      DEBOUTLN( cerr, "(hensel) lcDD = " << lcDD );
737      TIMING_START (ez_hensel_lift)
738      gcdfound= Hensel (B*lcD, DD, b, lcDD);
739      TIMING_END_AND_PRINT (ez_hensel_lift, "time to hensel lift in EZ: ")
740      DEBOUTLN( cerr, "(hensel finished) DD   = " << DD );
741
742      if (gcdfound == -1)
743      {
744        Off (SW_USE_EZGCD);
745        result= gcd (F,G);
746        On (SW_USE_EZGCD);
747        if (!isRat)
748          Off (SW_RATIONAL);
749        DEBDECLEVEL( cerr, "ezgcd" );
750        result /= icontent (result);
751        return N (d*result);
752      }
753
754      if (gcdfound)
755      {
756        TIMING_START (ez_termination)
757        contcand= content (DD[2], Variable (1));
758        cand = DD[2] / contcand;
759        if (B_is_F)
760          gcdfound = fdivides( cand, G ) && cand*(DD[1]/(lcD/contcand)) == F;
761        else
762          gcdfound = fdivides( cand, F ) && cand*(DD[1]/(lcD/contcand)) == G;
763        TIMING_END_AND_PRINT (ez_termination,
764                              "time for termination test in EZ: ")
765      }
766      /// ---> A8 (gcdfound)
767    }
768    delta--;
769  }
770  /// ---> A9
771  DEBDECLEVEL( cerr, "ezgcd" );
772  cand *= bCommonDen (cand);
773  if (!isRat)
774    Off (SW_RATIONAL);
775  cand /= icontent (cand);
776  return N (d*cand);
777}
778#endif
779
780/// Extended Zassenhaus GCD over Z.
781/// In case things become too dense we switch to a modular algorithm.
782CanonicalForm
783ezgcd ( const CanonicalForm & FF, const CanonicalForm & GG )
784{
785#ifdef HAVE_NTL
786  REvaluation b;
787  return ezgcd( FF, GG, b, false );
788#else
789  Off (SW_USE_EZGCD);
790  return gcd (FF, GG);
791  On (SW_USE_EZGCD);
792#endif
793}
794
795#ifdef HAVE_NTL
796// parameters for heuristic
797static int maxNumEval= 200;
798static int sizePerVars1= 500; //try dense gcd if size/#variables is bigger
799
800/// Extended Zassenhaus GCD for finite fields.
801/// In case things become too dense we switch to a modular algorithm.
802CanonicalForm EZGCD_P( const CanonicalForm & FF, const CanonicalForm & GG )
803{
804  if (FF.isZero() && degree(GG) > 0) return GG/Lc(GG);
805  else if (GG.isZero() && degree (FF) > 0) return FF/Lc(FF);
806  else if (FF.isZero() && GG.isZero()) return FF.genOne();
807  if (FF.inBaseDomain() || GG.inBaseDomain()) return FF.genOne();
808  if (FF.isUnivariate() && fdivides(FF, GG)) return FF/Lc(FF);
809  if (GG.isUnivariate() && fdivides(GG, FF)) return GG/Lc(GG);
810  if (FF == GG) return FF/Lc(FF);
811
812  int maxNumVars= tmax (getNumVars (FF), getNumVars (GG));
813  Variable a, oldA;
814  int sizeF= size (FF);
815  int sizeG= size (GG);
816
817  if (sizeF/maxNumVars > sizePerVars1 && sizeG/maxNumVars > sizePerVars1)
818  {
819    if (hasFirstAlgVar (FF, a) || hasFirstAlgVar (GG, a))
820      return modGCDFq (FF, GG, a);
821    else if (CFFactory::gettype() == GaloisFieldDomain)
822      return modGCDGF (FF, GG);
823    else
824      return modGCDFp (FF, GG);
825  }
826
827  CanonicalForm F, G, f, g, d, Fb, Gb, Db, Fbt, Gbt, Dbt, B0, B, D0, lcF, lcG,
828                lcD;
829  CFArray DD( 1, 2 ), lcDD( 1, 2 );
830  int degF, degG, delta, count;
831  int maxeval;
832  maxeval= tmin((getCharacteristic()/
833                (int)(ilog2(getCharacteristic())*log2exp))*2, maxNumEval);
834  count= 0; // number of eval. used
835  REvaluation b, bt;
836  int gcdfound = 0;
837  Variable x = Variable(1);
838
839  F= FF;
840  G= GG;
841
842  CFMap M,N;
843  int smallestDegLev;
844  TIMING_START (ez_p_compress)
845  int best_level= compress4EZGCD (F, G, M, N, smallestDegLev);
846
847  if (best_level == 0) return G.genOne();
848
849  F= M (F);
850  G= M (G);
851  TIMING_END_AND_PRINT (ez_p_compress, "time for compression in EZ_P: ")
852
853  TIMING_START (ez_p_content)
854  f = content( F, x ); g = content( G, x );
855  d = gcd( f, g );
856  F /= f; G /= g;
857  TIMING_END_AND_PRINT (ez_p_content, "time to extract content in EZ_P: ")
858
859  if( F.isUnivariate() && G.isUnivariate() )
860  {
861    if( F.mvar() == G.mvar() )
862      d *= gcd( F, G );
863    else
864      return N (d);
865    return N (d);
866  }
867  if ( F.isUnivariate())
868  {
869    g= content (G,G.mvar());
870    return N(d*gcd(F,g));
871  }
872  if ( G.isUnivariate())
873  {
874    f= content (F,F.mvar());
875    return N(d*gcd(G,f));
876  }
877
878  maxNumVars= tmax (getNumVars (F), getNumVars (G));
879  sizeF= size (F);
880  sizeG= size (G);
881
882  if (sizeF/maxNumVars > sizePerVars1 && sizeG/maxNumVars > sizePerVars1)
883  {
884    if (hasFirstAlgVar (F, a) || hasFirstAlgVar (G, a))
885      return N (d*modGCDFq (F, G, a));
886    else if (CFFactory::gettype() == GaloisFieldDomain)
887      return N (d*modGCDGF (F, G));
888    else
889      return N (d*modGCDFp (F, G));
890  }
891
892  int dummy= 0;
893  if( gcd_test_one( F, G, false, dummy ) )
894  {
895    return N (d);
896  }
897
898  bool passToGF= false;
899  bool extOfExt= false;
900  int p= getCharacteristic();
901  bool algExtension= (hasFirstAlgVar(F,a) || hasFirstAlgVar(G,a));
902  int k= 1;
903  CanonicalForm primElem, imPrimElem;
904  CFList source, dest;
905  if (p < 50 && CFFactory::gettype() != GaloisFieldDomain && !algExtension)
906  {
907    if (p == 2)
908      setCharacteristic (2, 12, 'Z');
909    else if (p == 3)
910      setCharacteristic (3, 4, 'Z');
911    else if (p == 5 || p == 7)
912      setCharacteristic (p, 3, 'Z');
913    else
914      setCharacteristic (p, 2, 'Z');
915    passToGF= true;
916    F= F.mapinto();
917    G= G.mapinto();
918    maxeval= 2*ipower (p, getGFDegree());
919  }
920  else if (CFFactory::gettype() == GaloisFieldDomain &&
921           ipower (p , getGFDegree()) < 50)
922  {
923    k= getGFDegree();
924    if (ipower (p, 2*k) > 50)
925      setCharacteristic (p, 2*k, gf_name);
926    else
927      setCharacteristic (p, 3*k, gf_name);
928    F= GFMapUp (F, k);
929    G= GFMapUp (G, k);
930    maxeval= tmin (2*ipower (p, getGFDegree()), maxNumEval);
931  }
932  else if (p < 50 && algExtension && CFFactory::gettype() != GaloisFieldDomain)
933  {
934    int d= degree (getMipo (a));
935    oldA= a;
936    Variable v2;
937    if (p == 2 && d < 6)
938    {
939      if (fac_NTL_char != p)
940      {
941        fac_NTL_char= p;
942        zz_p::init (p);
943      }
944      bool primFail= false;
945      Variable vBuf;
946      primElem= primitiveElement (a, vBuf, primFail);
947      ASSERT (!primFail, "failure in integer factorizer");
948      if (d < 3)
949      {
950        zz_pX NTLIrredpoly;
951        BuildIrred (NTLIrredpoly, d*3);
952        CanonicalForm newMipo= convertNTLzzpX2CF (NTLIrredpoly, Variable (1));
953        v2= rootOf (newMipo);
954      }
955      else
956      {
957        zz_pX NTLIrredpoly;
958        BuildIrred (NTLIrredpoly, d*2);
959        CanonicalForm newMipo= convertNTLzzpX2CF (NTLIrredpoly, Variable (1));
960        v2= rootOf (newMipo);
961      }
962      imPrimElem= mapPrimElem (primElem, a, v2);
963      extOfExt= true;
964    }
965    else if ((p == 3 && d < 4) || ((p == 5 || p == 7) && d < 3))
966    {
967      if (fac_NTL_char != p)
968      {
969        fac_NTL_char= p;
970        zz_p::init (p);
971      }
972      bool primFail= false;
973      Variable vBuf;
974      primElem= primitiveElement (a, vBuf, primFail);
975      ASSERT (!primFail, "failure in integer factorizer");
976      zz_pX NTLIrredpoly;
977      BuildIrred (NTLIrredpoly, d*2);
978      CanonicalForm newMipo= convertNTLzzpX2CF (NTLIrredpoly, Variable (1));
979      v2= rootOf (newMipo);
980      imPrimElem= mapPrimElem (primElem, a, v2);
981      extOfExt= true;
982    }
983    if (extOfExt)
984    {
985      maxeval= tmin (2*ipower (p, degree (getMipo (v2))), maxNumEval);
986      F= mapUp (F, a, v2, primElem, imPrimElem, source, dest);
987      G= mapUp (G, a, v2, primElem, imPrimElem, source, dest);
988      a= v2;
989    }
990  }
991
992  lcF = LC( F, x ); lcG = LC( G, x );
993  lcD = gcd( lcF, lcG );
994
995  delta = 0;
996  degF = degree( F, x ); degG = degree( G, x );
997
998  if (algExtension)
999    b = REvaluation( 2, tmax(F.level(), G.level()), AlgExtRandomF( a ) );
1000  else
1001  { // both not in extension given by algebraic variable
1002    if (CFFactory::gettype() != GaloisFieldDomain)
1003      b = REvaluation( 2, tmax(F.level(), G.level()), FFRandom() );
1004    else
1005      b = REvaluation( 2, tmax(F.level(), G.level()), GFRandom() );
1006  }
1007
1008  CanonicalForm cand, contcand;
1009  CanonicalForm result;
1010  int o, t;
1011  o= 0;
1012  t= 1;
1013  int goodPointCount= 0;
1014  while( !gcdfound )
1015  {
1016    TIMING_START (ez_p_eval);
1017    if( !findeval( F, G, Fb, Gb, Db, b, delta, degF, degG, maxeval, count, o,
1018         maxeval/maxNumVars, t ))
1019    { // too many eval. used --> try another method
1020      Off (SW_USE_EZGCD_P);
1021      result= gcd (F,G);
1022      On (SW_USE_EZGCD_P);
1023      if (passToGF)
1024      {
1025        CanonicalForm mipo= gf_mipo;
1026        setCharacteristic (p);
1027        Variable alpha= rootOf (mipo.mapinto());
1028        result= GF2FalphaRep (result, alpha);
1029        prune (alpha);
1030      }
1031      if (k > 1)
1032      {
1033        result= GFMapDown (result, k);
1034        setCharacteristic (p, k, gf_name);
1035      }
1036      if (extOfExt)
1037      {
1038        result= mapDown (result, primElem, imPrimElem, oldA, dest, source);
1039        prune1 (oldA);
1040      }
1041      return N (d*result);
1042    }
1043    TIMING_END_AND_PRINT (ez_p_eval, "time for eval point search in EZ_P1: ");
1044    delta = degree( Db );
1045    if (delta == degF)
1046    {
1047      if (degF <= degG && fdivides (F, G))
1048      {
1049        if (passToGF)
1050        {
1051          CanonicalForm mipo= gf_mipo;
1052          setCharacteristic (p);
1053          Variable alpha= rootOf (mipo.mapinto());
1054          F= GF2FalphaRep (F, alpha);
1055          prune (alpha);
1056        }
1057        if (k > 1)
1058        {
1059          F= GFMapDown (F, k);
1060          setCharacteristic (p, k, gf_name);
1061        }
1062        if (extOfExt)
1063        {
1064          F= mapDown (F, primElem, imPrimElem, oldA, dest, source);
1065          prune1 (oldA);
1066        }
1067        return N (d*F);
1068      }
1069      else
1070        delta--;
1071    }
1072    else if (delta == degG)
1073    {
1074      if (degG <= degF && fdivides (G, F))
1075      {
1076        if (passToGF)
1077        {
1078          CanonicalForm mipo= gf_mipo;
1079          setCharacteristic (p);
1080          Variable alpha= rootOf (mipo.mapinto());
1081          G= GF2FalphaRep (G, alpha);
1082          prune (alpha);
1083        }
1084        if (k > 1)
1085        {
1086          G= GFMapDown (G, k);
1087          setCharacteristic (p, k, gf_name);
1088        }
1089        if (extOfExt)
1090        {
1091          G= mapDown (G, primElem, imPrimElem, oldA, dest, source);
1092          prune1 (oldA);
1093        }
1094        return N (d*G);
1095      }
1096      else
1097        delta--;
1098    }
1099    if( delta == 0 )
1100    {
1101      if (passToGF)
1102        setCharacteristic (p);
1103      if (k > 1)
1104        setCharacteristic (p, k, gf_name);
1105      return N (d);
1106    }
1107    while( true )
1108    {
1109      bt = b;
1110      TIMING_START (ez_p_eval);
1111      if( !findeval(F,G,Fbt,Gbt,Dbt, bt, delta, degF, degG, maxeval, count, o,
1112           maxeval/maxNumVars, t ))
1113      { // too many eval. used --> try another method
1114        Off (SW_USE_EZGCD_P);
1115        result= gcd (F,G);
1116        On (SW_USE_EZGCD_P);
1117        if (passToGF)
1118        {
1119          CanonicalForm mipo= gf_mipo;
1120          setCharacteristic (p);
1121          Variable alpha= rootOf (mipo.mapinto());
1122          result= GF2FalphaRep (result, alpha);
1123          prune (alpha);
1124        }
1125        if (k > 1)
1126        {
1127          result= GFMapDown (result, k);
1128          setCharacteristic (p, k, gf_name);
1129        }
1130        if (extOfExt)
1131        {
1132          result= mapDown (result, primElem, imPrimElem, oldA, dest, source);
1133          prune1 (oldA);
1134        }
1135        return N (d*result);
1136      }
1137      TIMING_END_AND_PRINT (ez_p_eval, "time for eval point search in EZ_P2: ");
1138      int dd = degree( Dbt );
1139      if( dd == 0 )
1140      {
1141        if (passToGF)
1142          setCharacteristic (p);
1143        if (k > 1)
1144          setCharacteristic (p, k, gf_name);
1145        return N (d);
1146      }
1147      if( dd == delta )
1148      {
1149        goodPointCount++;
1150        if (goodPointCount == 5)
1151          break;
1152      }
1153      if( dd < delta )
1154      {
1155        goodPointCount= 0;
1156        delta = dd;
1157        b = bt;
1158        Db = Dbt; Fb = Fbt; Gb = Gbt;
1159      }
1160      if (delta == degF)
1161      {
1162        if (degF <= degG && fdivides (F, G))
1163        {
1164          if (passToGF)
1165          {
1166            CanonicalForm mipo= gf_mipo;
1167            setCharacteristic (p);
1168            Variable alpha= rootOf (mipo.mapinto());
1169            F= GF2FalphaRep (F, alpha);
1170            prune (alpha);
1171          }
1172          if (k > 1)
1173          {
1174            F= GFMapDown (F, k);
1175            setCharacteristic (p, k, gf_name);
1176          }
1177          if (extOfExt)
1178          {
1179            F= mapDown (F, primElem, imPrimElem, oldA, dest, source);
1180            prune1 (oldA);
1181          }
1182          return N (d*F);
1183        }
1184        else
1185          delta--;
1186      }
1187      else if (delta == degG)
1188      {
1189        if (degG <= degF && fdivides (G, F))
1190        {
1191          if (passToGF)
1192          {
1193            CanonicalForm mipo= gf_mipo;
1194            setCharacteristic (p);
1195            Variable alpha= rootOf (mipo.mapinto());
1196            G= GF2FalphaRep (G, alpha);
1197            prune (alpha);
1198          }
1199          if (k > 1)
1200          {
1201            G= GFMapDown (G, k);
1202            setCharacteristic (p, k, gf_name);
1203          }
1204          if (extOfExt)
1205          {
1206            G= mapDown (G, primElem, imPrimElem, oldA, dest, source);
1207            prune1 (oldA);
1208          }
1209          return N (d*G);
1210        }
1211        else
1212          delta--;
1213      }
1214      if( delta == 0 )
1215      {
1216        if (passToGF)
1217          setCharacteristic (p);
1218        if (k > 1)
1219          setCharacteristic (p, k, gf_name);
1220        return N (d);
1221      }
1222    }
1223    if( delta != degF && delta != degG )
1224    {
1225      bool B_is_F;
1226      CanonicalForm xxx1, xxx2;
1227      CanonicalForm buf;
1228      DD[1] = Fb / Db;
1229      buf= Gb/Db;
1230      xxx1 = gcd( DD[1], Db );
1231      xxx2 = gcd( buf, Db );
1232      if (((xxx1.inCoeffDomain() && xxx2.inCoeffDomain()) &&
1233          (size (F) <= size (G)))
1234          || (xxx1.inCoeffDomain() && !xxx2.inCoeffDomain()))
1235      {
1236        B = F;
1237        DD[2] = Db;
1238        lcDD[1] = lcF;
1239        lcDD[2] = lcD;
1240        B_is_F = true;
1241      }
1242      else if (((xxx1.inCoeffDomain() && xxx2.inCoeffDomain()) &&
1243               (size (G) < size (F)))
1244               || (!xxx1.inCoeffDomain() && xxx2.inCoeffDomain()))
1245      {
1246        DD[1] = buf;
1247        B = G;
1248        DD[2] = Db;
1249        lcDD[1] = lcG;
1250        lcDD[2] = lcD;
1251        B_is_F = false;
1252      }
1253      else // special case handling
1254      {
1255        Off (SW_USE_EZGCD_P);
1256        result= gcd (F,G);
1257        On (SW_USE_EZGCD_P);
1258        if (passToGF)
1259        {
1260          CanonicalForm mipo= gf_mipo;
1261          setCharacteristic (p);
1262          Variable alpha= rootOf (mipo.mapinto());
1263          result= GF2FalphaRep (result, alpha);
1264          prune (alpha);
1265        }
1266        if (k > 1)
1267        {
1268          result= GFMapDown (result, k);
1269          setCharacteristic (p, k, gf_name);
1270        }
1271        if (extOfExt)
1272        {
1273          result= mapDown (result, primElem, imPrimElem, oldA, dest, source);
1274          prune1 (oldA);
1275        }
1276        return N (d*result);
1277      }
1278      DD[2] = DD[2] * ( b( lcDD[2] ) / lc( DD[2] ) );
1279      DD[1] = DD[1] * ( b( lcDD[1] ) / lc( DD[1] ) );
1280
1281      if (size (B*lcDD[2])/maxNumVars > sizePerVars1)
1282      {
1283        if (algExtension)
1284        {
1285          result= modGCDFq (F, G, a);
1286          if (extOfExt)
1287          {
1288            result= mapDown (result, primElem, imPrimElem, oldA, dest, source);
1289            prune1 (oldA);
1290          }
1291          return N (d*result);
1292        }
1293        if (CFFactory::gettype() == GaloisFieldDomain)
1294        {
1295          result= modGCDGF (F, G);
1296          if (passToGF)
1297          {
1298            CanonicalForm mipo= gf_mipo;
1299            setCharacteristic (p);
1300            Variable alpha= rootOf (mipo.mapinto());
1301            result= GF2FalphaRep (result, alpha);
1302            prune (alpha);
1303          }
1304          if (k > 1)
1305          {
1306            result= GFMapDown (result, k);
1307            setCharacteristic (p, k, gf_name);
1308          }
1309          return N (d*result);
1310        }
1311        else
1312          return N (d*modGCDFp (F,G));
1313      }
1314
1315      TIMING_START (ez_p_hensel_lift);
1316      gcdfound= Hensel (B*lcD, DD, b, lcDD);
1317      TIMING_END_AND_PRINT (ez_p_hensel_lift, "time for Hensel lift in EZ_P: ");
1318
1319      if (gcdfound == -1) //things became dense
1320      {
1321        if (algExtension)
1322        {
1323          result= modGCDFq (F, G, a);
1324          if (extOfExt)
1325          {
1326            result= mapDown (result, primElem, imPrimElem, oldA, dest, source);
1327            prune1 (oldA);
1328          }
1329          return N (d*result);
1330        }
1331        if (CFFactory::gettype() == GaloisFieldDomain)
1332        {
1333          result= modGCDGF (F, G);
1334          if (passToGF)
1335          {
1336            CanonicalForm mipo= gf_mipo;
1337            setCharacteristic (p);
1338            Variable alpha= rootOf (mipo.mapinto());
1339            result= GF2FalphaRep (result, alpha);
1340            prune (alpha);
1341          }
1342          if (k > 1)
1343          {
1344            result= GFMapDown (result, k);
1345            setCharacteristic (p, k, gf_name);
1346          }
1347          return N (d*result);
1348        }
1349        else
1350        {
1351          if (p >= cf_getBigPrime(0))
1352            return N (d*sparseGCDFp (F,G));
1353          else
1354            return N (d*modGCDFp (F,G));
1355        }
1356      }
1357
1358      if (gcdfound == 1)
1359      {
1360        TIMING_START (termination_test);
1361        contcand= content (DD[2], Variable (1));
1362        cand = DD[2] / contcand;
1363        if (B_is_F)
1364          gcdfound = fdivides( cand, G ) && cand*(DD[1]/(lcD/contcand)) == F;
1365        else
1366          gcdfound = fdivides( cand, F ) && cand*(DD[1]/(lcD/contcand)) == G;
1367        TIMING_END_AND_PRINT (termination_test,
1368                              "time for termination test EZ_P: ");
1369
1370        if (passToGF && gcdfound)
1371        {
1372          CanonicalForm mipo= gf_mipo;
1373          setCharacteristic (p);
1374          Variable alpha= rootOf (mipo.mapinto());
1375          cand= GF2FalphaRep (cand, alpha);
1376          prune (alpha);
1377        }
1378        if (k > 1 && gcdfound)
1379        {
1380          cand= GFMapDown (cand, k);
1381          setCharacteristic (p, k, gf_name);
1382        }
1383        if (extOfExt && gcdfound)
1384        {
1385          cand= mapDown (cand, primElem, imPrimElem, oldA, dest, source);
1386          prune1 (oldA);
1387        }
1388      }
1389    }
1390    delta--;
1391    goodPointCount= 0;
1392  }
1393  return N (d*cand);
1394}
1395#endif
1396
Note: See TracBrowser for help on using the repository browser.