source: git/factory/cfEzgcd.cc @ 11a6ca

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