source: git/factory/fac_ezgcd.cc @ 2024b69

spielwiese
Last change on this file since 2024b69 was 2a95b2, checked in by Martin Lee <martinlee84@…>, 12 years ago
chg: added more timing infos to gcd functions
  • Property mode set to 100644
File size: 16.5 KB
Line 
1/*****************************************************************************\
2 * Computer Algebra System SINGULAR
3\*****************************************************************************/
4/** @file fac_ezgcd.cc
5 *
6 * This file implements the GCD of two multivariate polynomials over Q using
7 * EZ-GCD as described in "Algorithms for Computer Algebra" by Geddes, Czapor,
8 * Labahnn
9 *
10 * @author Martin Lee
11 *
12 **/
13/*****************************************************************************/
14
15#include "config.h"
16
17#include "timing.h"
18#include "cf_assert.h"
19#include "debug.h"
20
21#include "cf_defs.h"
22#include "canonicalform.h"
23#include "fac_util.h" // HEADER for this file
24#include "cf_algorithm.h"
25#include "cf_reval.h"
26#include "cf_random.h"
27#include "cf_primes.h"
28#include "templates/ftmpl_functions.h"
29#include "cf_map.h"
30#include "facHensel.h"
31
32#ifdef HAVE_NTL
33
34TIMING_DEFINE_PRINT(ez_eval)
35TIMING_DEFINE_PRINT(ez_compress)
36TIMING_DEFINE_PRINT(ez_hensel_lift)
37TIMING_DEFINE_PRINT(ez_content)
38TIMING_DEFINE_PRINT(ez_termination)
39
40static
41int compress4EZGCD (const CanonicalForm& F, const CanonicalForm& G, CFMap & M,
42                    CFMap & N, int& both_non_zero)
43{
44  int n= tmax (F.level(), G.level());
45  int * degsf= new int [n + 1];
46  int * degsg= new int [n + 1];
47
48  for (int i = 0; i <= n; i++)
49    degsf[i]= degsg[i]= 0;
50
51  degsf= degrees (F, degsf);
52  degsg= degrees (G, degsg);
53
54  both_non_zero= 0;
55  int f_zero= 0;
56  int g_zero= 0;
57
58  for (int i= 1; i <= n; i++)
59  {
60    if (degsf[i] != 0 && degsg[i] != 0)
61    {
62      both_non_zero++;
63      continue;
64    }
65    if (degsf[i] == 0 && degsg[i] != 0 && i <= G.level())
66    {
67      f_zero++;
68      continue;
69    }
70    if (degsg[i] == 0 && degsf[i] && i <= F.level())
71    {
72      g_zero++;
73      continue;
74    }
75  }
76
77  if (both_non_zero == 0)
78  {
79    delete [] degsf;
80    delete [] degsg;
81    return 0;
82  }
83
84  // map Variables which do not occur in both polynomials to higher levels
85  int k= 1;
86  int l= 1;
87  for (int i= 1; i <= n; i++)
88  {
89    if (degsf[i] != 0 && degsg[i] == 0 && i <= F.level())
90    {
91      if (k + both_non_zero != i)
92      {
93        M.newpair (Variable (i), Variable (k + both_non_zero));
94        N.newpair (Variable (k + both_non_zero), Variable (i));
95      }
96      k++;
97    }
98    if (degsf[i] == 0 && degsg[i] != 0 && i <= G.level())
99    {
100      if (l + g_zero + both_non_zero != i)
101      {
102        M.newpair (Variable (i), Variable (l + g_zero + both_non_zero));
103        N.newpair (Variable (l + g_zero + both_non_zero), Variable (i));
104      }
105      l++;
106    }
107  }
108
109  // sort Variables x_{i} in decreasing order of
110  // min(deg_{x_{i}}(f),deg_{x_{i}}(g))
111  int m= tmin (F.level(), G.level());
112  int max_min_deg;
113  k= both_non_zero;
114  l= 0;
115  int i= 1;
116  while (k > 0)
117  {
118    max_min_deg= tmin (degsf[i], degsg[i]);
119    while (max_min_deg == 0)
120    {
121      i++;
122      max_min_deg= tmin (degsf[i], degsg[i]);
123    }
124    for (int j= i + 1; j <= m; j++)
125    {
126      if ((tmin (degsf[j],degsg[j]) < max_min_deg) &&
127          (tmin (degsf[j], degsg[j]) != 0))
128      {
129        max_min_deg= tmin (degsf[j], degsg[j]);
130        l= j;
131      }
132    }
133
134    if (l != 0)
135    {
136      if (l != k)
137      {
138        M.newpair (Variable (l), Variable(k));
139        N.newpair (Variable (k), Variable(l));
140        degsf[l]= 0;
141        degsg[l]= 0;
142        l= 0;
143      }
144      else
145      {
146        degsf[l]= 0;
147        degsg[l]= 0;
148        l= 0;
149      }
150    }
151    else if (l == 0)
152    {
153      if (i != k)
154      {
155        M.newpair (Variable (i), Variable (k));
156        N.newpair (Variable (k), Variable (i));
157        degsf[i]= 0;
158        degsg[i]= 0;
159      }
160      else
161      {
162        degsf[i]= 0;
163        degsg[i]= 0;
164      }
165      i++;
166    }
167    k--;
168  }
169
170  delete [] degsf;
171  delete [] degsg;
172
173  return both_non_zero;
174}
175
176static inline
177CanonicalForm myShift2Zero (const CanonicalForm& F, CFList& Feval,
178                            const CFList& evaluation)
179{
180  CanonicalForm A= F;
181  int k= 2;
182  for (CFListIterator i= evaluation; i.hasItem(); i++, k++)
183    A= A (Variable (k) + i.getItem(), k);
184
185  CanonicalForm buf= A;
186  Feval= CFList();
187  Feval.append (buf);
188  for (k= evaluation.length() + 1; k > 2; k--)
189  {
190    buf= mod (buf, Variable (k));
191    Feval.insert (buf);
192  }
193  return A;
194}
195
196static inline
197CanonicalForm myReverseShift (const CanonicalForm& F, const CFList& evaluation)
198{
199  int l= evaluation.length() + 1;
200  CanonicalForm result= F;
201  CFListIterator j= evaluation;
202  for (int i= 2; i < l + 1; i++, j++)
203  {
204    if (F.level() < i)
205      continue;
206    result= result (Variable (i) - j.getItem(), i);
207  }
208  return result;
209}
210
211static inline
212Evaluation optimize4Lift (const CanonicalForm& F, CFMap & M,
213                    CFMap & N, const Evaluation& A)
214{
215  int n= F.level();
216  int * degsf= new int [n + 1];
217
218  for (int i = 0; i <= n; i++)
219    degsf[i]= 0;
220
221  degsf= degrees (F, degsf);
222
223  Evaluation result= Evaluation (A.min(), A.max());
224  ASSERT (A.min() == 2, "expected A.min() == 2");
225  int max_deg;
226  int k= n;
227  int l= 1;
228  int i= 2;
229  int pos= 2;
230  while (k > 1)
231  {
232    max_deg= degsf [i];
233    while (max_deg == 0)
234    {
235      i++;
236      max_deg= degsf [i];
237    }
238    l= i;
239    for (int j= i + 1; j <=  n; j++)
240    {
241      if (degsf[j] > max_deg)
242      {
243        max_deg= degsf[j];
244        l= j;
245      }
246    }
247
248    if (l <= n)
249    {
250      if (l != pos)
251      {
252        result.setValue (pos, A [l]);
253        M.newpair (Variable (l), Variable (pos));
254        N.newpair (Variable (pos), Variable (l));
255        degsf[l]= 0;
256        l= 2;
257        if (k == 2 && n == 3)
258        {
259          result.setValue (l, A [pos]);
260          M.newpair (Variable (pos), Variable (l));
261          N.newpair (Variable (l), Variable (pos));
262          degsf[pos]= 0;
263        }
264      }
265      else
266      {
267        result.setValue (l, A [l]);
268        degsf [l]= 0;
269      }
270    }
271    pos++;
272    k--;
273    l= 2;
274  }
275
276  delete [] degsf;
277
278  return result;
279}
280
281static inline
282int Hensel (const CanonicalForm & UU, CFArray & G, const Evaluation & AA,
283            const CFArray& LeadCoeffs )
284{
285  CFList factors;
286  factors.append (G[1]);
287  factors.append (G[2]);
288
289  CFMap NN, MM;
290  Evaluation A= optimize4Lift (UU, MM, NN, AA);
291
292  CanonicalForm U= MM (UU);
293  CFArray LCs= CFArray (1,2);
294  LCs [1]= MM (LeadCoeffs [1]);
295  LCs [2]= MM (LeadCoeffs [2]);
296
297  CFList evaluation;
298  for (int i= A.min(); i <= A.max(); i++)
299    evaluation.append (A [i]);
300  CFList UEval;
301  CanonicalForm shiftedU= myShift2Zero (U, UEval, evaluation);
302
303  if (size (shiftedU)/getNumVars (U) > 500)
304    return -1;
305
306  CFArray shiftedLCs= CFArray (2);
307  CFList shiftedLCsEval1, shiftedLCsEval2;
308  shiftedLCs[0]= myShift2Zero (LCs[1], shiftedLCsEval1, evaluation);
309  shiftedLCs[1]= myShift2Zero (LCs[2], shiftedLCsEval2, evaluation);
310  factors.insert (1);
311  int liftBound= degree (UEval.getLast(), 2) + 1;
312  CFArray Pi;
313  CFMatrix M= CFMatrix (liftBound, factors.length() - 1);
314  CFList diophant;
315  CFArray lcs= CFArray (2);
316  lcs [0]= shiftedLCsEval1.getFirst();
317  lcs [1]= shiftedLCsEval2.getFirst();
318  nonMonicHenselLift12 (UEval.getFirst(), factors, liftBound, Pi, diophant, M,
319                        lcs, false);
320
321  for (CFListIterator i= factors; i.hasItem(); i++)
322  {
323    if (!fdivides (i.getItem(), UEval.getFirst()))
324      return 0;
325  }
326
327  int * liftBounds;
328  bool noOneToOne= false;
329  if (U.level() > 2)
330  {
331    liftBounds= new int [U.level() - 1]; /* index: 0.. U.level()-2 */
332    liftBounds[0]= liftBound;
333    for (int i= 1; i < U.level() - 1; i++)
334      liftBounds[i]= degree (shiftedU, Variable (i + 2)) + 1;
335    factors= nonMonicHenselLift2 (UEval, factors, liftBounds, U.level() - 1,
336                                  false, shiftedLCsEval1, shiftedLCsEval2, Pi,
337                                  diophant, noOneToOne);
338    delete [] liftBounds;
339    if (noOneToOne)
340      return 0;
341  }
342  G[1]= factors.getFirst();
343  G[2]= factors.getLast();
344  G[1]= myReverseShift (G[1], evaluation);
345  G[2]= myReverseShift (G[2], evaluation);
346  G[1]= NN (G[1]);
347  G[2]= NN (G[2]);
348  return 1;
349}
350
351static
352bool findeval (const CanonicalForm & F, const CanonicalForm & G,
353               CanonicalForm & Fb, CanonicalForm & Gb, CanonicalForm & Db,
354               REvaluation & b, int delta, int degF, int degG, int maxeval,
355               int & count, int& k, int bound, int& l)
356{
357  if( count == 0 && delta != 0)
358  {
359    if( count++ > maxeval )
360      return false;
361  }
362  if (count > 0)
363  {
364    b.nextpoint(k);
365    if (k == 0)
366      k++;
367    l++;
368    if (l > bound)
369    {
370      l= 1;
371      k++;
372      if (k > tmax (F.level(), G.level()) - 1)
373        return false;
374      b.nextpoint (k);
375    }
376    if (count++ > maxeval)
377      return false;
378  }
379  while( true )
380  {
381    Fb = b( F );
382    if( degree( Fb, 1 ) == degF )
383    {
384      Gb = b( G );
385      if( degree( Gb, 1 ) == degG )
386      {
387        Db = gcd( Fb, Gb );
388        if( delta > 0 )
389        {
390          if( degree( Db, 1 ) <= delta )
391            return true;
392        }
393        else
394        {
395          k++;
396          return true;
397        }
398      }
399    }
400    if (k == 0)
401      k++;
402    b.nextpoint(k);
403    l++;
404    if (l > bound)
405    {
406      l= 1;
407      k++;
408      if (k > tmax (F.level(), G.level()) - 1)
409        return false;
410      b.nextpoint (k);
411    }
412    if( count++ > maxeval )
413      return false;
414  }
415}
416
417static CanonicalForm
418ezgcd ( const CanonicalForm & FF, const CanonicalForm & GG, REvaluation & b,
419        bool internal )
420{
421  bool isRat= isOn (SW_RATIONAL);
422  CanonicalForm F, G, f, g, d, Fb, Gb, Db, Fbt, Gbt, Dbt, B0, B, D0, lcF, lcG,
423                lcD, cand, contcand, result;
424  CFArray DD( 1, 2 ), lcDD( 1, 2 );
425  int degF, degG, delta, t, count, maxeval;
426  REvaluation bt;
427  bool gcdfound = false;
428  Variable x = Variable(1);
429  count= 0;
430  maxeval= 200;
431  int o, l;
432  o= 0;
433  l= 1;
434
435  if (!isRat)
436    On (SW_RATIONAL);
437  F= FF*bCommonDen (FF);
438  G= GG*bCommonDen (GG);
439  if (!isRat)
440    Off (SW_RATIONAL);
441
442  TIMING_START (ez_compress)
443  CFMap M,N;
444  int smallestDegLev;
445  int best_level= compress4EZGCD (F, G, M, N, smallestDegLev);
446
447  if (best_level == 0)
448  {
449    DEBDECLEVEL( cerr, "ezgcd" );
450    return G.genOne();
451  }
452
453  F= M (F);
454  G= M (G);
455  TIMING_END_AND_PRINT (ez_compress, "time for compression in EZ: ")
456
457  DEBINCLEVEL( cerr, "ezgcd" );
458  DEBOUTLN( cerr, "FF = " << FF );
459  DEBOUTLN( cerr, "GG = " << GG );
460  TIMING_START (ez_content)
461  f = content( F, x ); g = content( G, x ); d = gcd( f, g );
462  d /= icontent (d);
463  DEBOUTLN( cerr, "f = " << f );
464  DEBOUTLN( cerr, "g = " << g );
465  F /= f; G /= g;
466  TIMING_END_AND_PRINT (ez_content, "time to extract content in EZ: ")
467  if ( F.isUnivariate() && G.isUnivariate() )
468  {
469    DEBDECLEVEL( cerr, "ezgcd" );
470    if(F.mvar()==G.mvar())
471      d*=gcd(F,G);
472    return N (d);
473  }
474
475  int maxNumVars= tmax (getNumVars (F), getNumVars (G));
476  int sizeF= size (F);
477  int sizeG= size (G);
478
479
480  if (!isRat)
481    On (SW_RATIONAL);
482  if (sizeF/maxNumVars > 500 && sizeG/maxNumVars > 500)
483  {
484    Off(SW_USE_EZGCD);
485    result=gcd( F, G );
486    On(SW_USE_EZGCD);
487    if (!isRat)
488      Off (SW_RATIONAL);
489    result /= icontent (result);
490    DEBDECLEVEL( cerr, "ezgcd" );
491    return N (d*result);
492  }
493
494  if ( gcd_test_one( F, G, false ) )
495  {
496    DEBDECLEVEL( cerr, "ezgcd" );
497    if (!isRat)
498      Off (SW_RATIONAL);
499    return N (d);
500  }
501  lcF = LC( F, x ); lcG = LC( G, x );
502  lcD = gcd( lcF, lcG );
503  delta = 0;
504  degF = degree( F, x ); degG = degree( G, x );
505  t = tmax( F.level(), G.level() );
506  if ( ! internal )
507    b = REvaluation( 2, t, IntRandom( 25 ) );
508  while ( ! gcdfound )
509  {
510    /// ---> A2
511    DEBOUTLN( cerr, "search for evaluation, delta = " << delta );
512    DEBOUTLN( cerr, "F = " << F );
513    DEBOUTLN( cerr, "G = " << G );
514    TIMING_START (ez_eval)
515    if (!findeval( F, G, Fb, Gb, Db, b, delta, degF, degG, maxeval, count,
516                   o, 25, l))
517    {
518      Off(SW_USE_EZGCD);
519      result=gcd( F, G );
520      On(SW_USE_EZGCD);
521      if (!isRat)
522        Off (SW_RATIONAL);
523      DEBDECLEVEL( cerr, "ezgcd" );
524      result /= icontent (result);
525      return N (d*result);
526    }
527    TIMING_END_AND_PRINT (ez_eval, "time to find eval point in EZ1: ")
528    DEBOUTLN( cerr, "found evaluation b = " << b );
529    DEBOUTLN( cerr, "F(b) = " << Fb );
530    DEBOUTLN( cerr, "G(b) = " << Gb );
531    DEBOUTLN( cerr, "D(b) = " << Db );
532    delta = degree( Db );
533    /// ---> A3
534    if ( delta == 0 )
535    {
536      DEBDECLEVEL( cerr, "ezgcd" );
537      if (!isRat)
538        Off (SW_RATIONAL);
539      return N (d);
540    }
541    /// ---> A4
542    //deltaold = delta;
543    while ( 1 ) 
544    {
545      bt = b;
546      TIMING_START (ez_eval)
547      if (!findeval( F, G, Fbt, Gbt, Dbt, bt, delta, degF, degG, maxeval, count,
548                     o, 25,l ))
549      {
550        Off(SW_USE_EZGCD);
551        result=gcd( F, G );
552        On(SW_USE_EZGCD);
553        if (!isRat)
554          Off (SW_RATIONAL);
555        DEBDECLEVEL( cerr, "ezgcd" );
556        result /= icontent (result);
557        return N (d*result);
558      }
559      TIMING_END_AND_PRINT (ez_eval, "time to find eval point in EZ2: ")
560      int dd=degree( Dbt );
561      if ( dd /*degree( Dbt )*/ == 0 )
562      {
563        DEBDECLEVEL( cerr, "ezgcd" );
564        if (!isRat)
565          Off (SW_RATIONAL);
566        return N (d);
567      }
568      if ( dd /*degree( Dbt )*/ == delta )
569        break;
570      else  if ( dd /*degree( Dbt )*/ < delta )
571      {
572        delta = dd /*degree( Dbt )*/;
573        b = bt;
574        Db = Dbt; Fb = Fbt; Gb = Gbt;
575      }
576      DEBOUTLN( cerr, "now after A4, delta = " << delta );
577      /// ---> A5
578      if (delta == degF)
579      {
580        if (degF <= degG  && fdivides (F, G))
581        {
582          DEBDECLEVEL( cerr, "ezgcd" );
583          if (!isRat)
584            Off (SW_RATIONAL);
585          return N (d*F);
586        }
587        else
588          delta--;
589      }
590      else if (delta == degG)
591      {
592        if (degG <= degF && fdivides( G, F ))
593        {
594          DEBDECLEVEL( cerr, "ezgcd" );
595          if (!isRat)
596            Off (SW_RATIONAL);
597          return N (d*G);
598        }
599        else
600          delta--;
601      }
602    }
603    if ( delta != degF && delta != degG )
604    {
605      /// ---> A6
606      bool B_is_F;
607      CanonicalForm xxx1, xxx2;
608      CanonicalForm buf;
609      DD[1] = Fb / Db;
610      buf= Gb/Db;
611      xxx1 = gcd( DD[1], Db );
612      xxx2 = gcd( buf, Db );
613      if (((xxx1.inCoeffDomain() && xxx2.inCoeffDomain()) &&
614          (size (F) <= size (G)))
615            || (xxx1.inCoeffDomain() && !xxx2.inCoeffDomain()))
616      {
617        B = F;
618        DD[2] = Db;
619        lcDD[1] = lcF;
620        lcDD[2] = lcD;
621        B_is_F = true;
622      }
623      else if (((xxx1.inCoeffDomain() && xxx2.inCoeffDomain()) &&
624                (size (G) < size (F)))
625                || (!xxx1.inCoeffDomain() && xxx2.inCoeffDomain()))
626      {
627        DD[1] = buf;
628        B = G;
629        DD[2] = Db;
630        lcDD[1] = lcG;
631        lcDD[2] = lcD;
632        B_is_F = false;
633      }
634      else
635      {
636        //special case
637        Off(SW_USE_EZGCD);
638        result=gcd( F, G );
639        On(SW_USE_EZGCD);
640        if (!isRat)
641          Off (SW_RATIONAL);
642        DEBDECLEVEL( cerr, "ezgcd" );
643        result /= icontent (result);
644        return N (d*result);
645      }
646      /// ---> A7
647      DD[2] = DD[2] * ( b( lcDD[2] ) / lc( DD[2] ) );
648      DD[1] = DD[1] * ( b( lcDD[1] ) / lc( DD[1] ) );
649      DEBOUTLN( cerr, "(hensel) B    = " << B );
650      DEBOUTLN( cerr, "(hensel) lcB  = " << LC( B, Variable(1) ) );
651      DEBOUTLN( cerr, "(hensel) b(B) = " << b(B) );
652      DEBOUTLN( cerr, "(hensel) DD   = " << DD );
653      DEBOUTLN( cerr, "(hensel) lcDD = " << lcDD );
654      TIMING_START (ez_hensel_lift)
655      gcdfound= Hensel (B*lcD, DD, b, lcDD);
656      TIMING_END_AND_PRINT (ez_hensel_lift, "time to hensel lift in EZ: ")
657      DEBOUTLN( cerr, "(hensel finished) DD   = " << DD );
658
659      if (gcdfound == -1)
660      {
661        Off (SW_USE_EZGCD);
662        result= gcd (F,G);
663        On (SW_USE_EZGCD);
664        if (!isRat)
665          Off (SW_RATIONAL);
666        DEBDECLEVEL( cerr, "ezgcd" );
667        result /= icontent (result);
668        return N (d*result);
669      }
670
671      if (gcdfound)
672      {
673        TIMING_START (ez_termination)
674        contcand= content (DD[2], Variable (1));
675        cand = DD[2] / contcand;
676        if (B_is_F)
677          gcdfound = fdivides( cand, G ) && cand*(DD[1]/(lcD/contcand)) == F;
678        else
679          gcdfound = fdivides( cand, F ) && cand*(DD[1]/(lcD/contcand)) == G;
680        TIMING_END_AND_PRINT (ez_termination,
681                              "time for termination test in EZ: ")
682      }
683      /// ---> A8 (gcdfound)
684    }
685    delta--;
686  }
687  /// ---> A9
688  DEBDECLEVEL( cerr, "ezgcd" );
689  cand *= bCommonDen (cand);
690  if (!isRat)
691    Off (SW_RATIONAL);
692  cand /= icontent (cand);
693  return N (d*cand);
694}
695#endif
696
697CanonicalForm
698ezgcd ( const CanonicalForm & FF, const CanonicalForm & GG )
699{
700#ifdef HAVE_NTL
701  REvaluation b;
702  return ezgcd( FF, GG, b, false );
703#else
704  Off (SW_USE_EZGCD);
705  return gcd (FF, GG);
706  On (SW_USE_EZGCD);
707#endif
708}
709
Note: See TracBrowser for help on using the repository browser.