source: git/factory/cfEzgcd.cc @ 928be8

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