source: git/factory/fac_ezgcd.cc @ 0b447e

spielwiese
Last change on this file since 0b447e was b770bf, checked in by Martin Lee <martinlee84@…>, 11 years ago
chg: faster gcd computation in EZGCD in corner cases
  • Property mode set to 100644
File size: 16.9 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  long termEstimate= size (U);
299  for (int i= A.min(); i <= A.max(); i++)
300  {
301    if (!A[i].isZero())
302    {
303      termEstimate *= degree (U,i)*2;
304      termEstimate /= 3;
305    }
306    evaluation.append (A [i]);
307  }
308  if (termEstimate/getNumVars(U) > 500)
309    return -1;
310  CFList UEval;
311  CanonicalForm shiftedU= myShift2Zero (U, UEval, evaluation);
312
313  if (size (shiftedU)/getNumVars (U) > 500)
314    return -1;
315
316  CFArray shiftedLCs= CFArray (2);
317  CFList shiftedLCsEval1, shiftedLCsEval2;
318  shiftedLCs[0]= myShift2Zero (LCs[1], shiftedLCsEval1, evaluation);
319  shiftedLCs[1]= myShift2Zero (LCs[2], shiftedLCsEval2, evaluation);
320  factors.insert (1);
321  int liftBound= degree (UEval.getLast(), 2) + 1;
322  CFArray Pi;
323  CFMatrix M= CFMatrix (liftBound, factors.length() - 1);
324  CFList diophant;
325  CFArray lcs= CFArray (2);
326  lcs [0]= shiftedLCsEval1.getFirst();
327  lcs [1]= shiftedLCsEval2.getFirst();
328  nonMonicHenselLift12 (UEval.getFirst(), factors, liftBound, Pi, diophant, M,
329                        lcs, false);
330
331  for (CFListIterator i= factors; i.hasItem(); i++)
332  {
333    if (!fdivides (i.getItem(), UEval.getFirst()))
334      return 0;
335  }
336
337  int * liftBounds;
338  bool noOneToOne= false;
339  if (U.level() > 2)
340  {
341    liftBounds= new int [U.level() - 1]; /* index: 0.. U.level()-2 */
342    liftBounds[0]= liftBound;
343    for (int i= 1; i < U.level() - 1; i++)
344      liftBounds[i]= degree (shiftedU, Variable (i + 2)) + 1;
345    factors= nonMonicHenselLift2 (UEval, factors, liftBounds, U.level() - 1,
346                                  false, shiftedLCsEval1, shiftedLCsEval2, Pi,
347                                  diophant, noOneToOne);
348    delete [] liftBounds;
349    if (noOneToOne)
350      return 0;
351  }
352  G[1]= factors.getFirst();
353  G[2]= factors.getLast();
354  G[1]= myReverseShift (G[1], evaluation);
355  G[2]= myReverseShift (G[2], evaluation);
356  G[1]= NN (G[1]);
357  G[2]= NN (G[2]);
358  return 1;
359}
360
361static
362bool findeval (const CanonicalForm & F, const CanonicalForm & G,
363               CanonicalForm & Fb, CanonicalForm & Gb, CanonicalForm & Db,
364               REvaluation & b, int delta, int degF, int degG, int maxeval,
365               int & count, int& k, int bound, int& l)
366{
367  if( count == 0 && delta != 0)
368  {
369    if( count++ > maxeval )
370      return false;
371  }
372  if (count > 0)
373  {
374    b.nextpoint(k);
375    if (k == 0)
376      k++;
377    l++;
378    if (l > bound)
379    {
380      l= 1;
381      k++;
382      if (k > tmax (F.level(), G.level()) - 1)
383        return false;
384      b.nextpoint (k);
385    }
386    if (count++ > maxeval)
387      return false;
388  }
389  while( true )
390  {
391    Fb = b( F );
392    if( degree( Fb, 1 ) == degF )
393    {
394      Gb = b( G );
395      if( degree( Gb, 1 ) == degG )
396      {
397        Db = gcd( Fb, Gb );
398        if( delta > 0 )
399        {
400          if( degree( Db, 1 ) <= delta )
401            return true;
402        }
403        else
404        {
405          k++;
406          return true;
407        }
408      }
409    }
410    if (k == 0)
411      k++;
412    b.nextpoint(k);
413    l++;
414    if (l > bound)
415    {
416      l= 1;
417      k++;
418      if (k > tmax (F.level(), G.level()) - 1)
419        return false;
420      b.nextpoint (k);
421    }
422    if( count++ > maxeval )
423      return false;
424  }
425}
426
427static CanonicalForm
428ezgcd ( const CanonicalForm & FF, const CanonicalForm & GG, REvaluation & b,
429        bool internal )
430{
431  bool isRat= isOn (SW_RATIONAL);
432  CanonicalForm F, G, f, g, d, Fb, Gb, Db, Fbt, Gbt, Dbt, B0, B, D0, lcF, lcG,
433                lcD, cand, contcand, result;
434  CFArray DD( 1, 2 ), lcDD( 1, 2 );
435  int degF, degG, delta, t, count, maxeval;
436  REvaluation bt;
437  int gcdfound = 0;
438  Variable x = Variable(1);
439  count= 0;
440  maxeval= 200;
441  int o, l;
442  o= 0;
443  l= 1;
444
445  if (!isRat)
446    On (SW_RATIONAL);
447  F= FF*bCommonDen (FF);
448  G= GG*bCommonDen (GG);
449  if (!isRat)
450    Off (SW_RATIONAL);
451
452  TIMING_START (ez_compress)
453  CFMap M,N;
454  int smallestDegLev;
455  int best_level= compress4EZGCD (F, G, M, N, smallestDegLev);
456
457  if (best_level == 0)
458  {
459    DEBDECLEVEL( cerr, "ezgcd" );
460    return G.genOne();
461  }
462
463  F= M (F);
464  G= M (G);
465  TIMING_END_AND_PRINT (ez_compress, "time for compression in EZ: ")
466
467  DEBINCLEVEL( cerr, "ezgcd" );
468  DEBOUTLN( cerr, "FF = " << FF );
469  DEBOUTLN( cerr, "GG = " << GG );
470  TIMING_START (ez_content)
471  f = content( F, x ); g = content( G, x ); d = gcd( f, g );
472  d /= icontent (d);
473  DEBOUTLN( cerr, "f = " << f );
474  DEBOUTLN( cerr, "g = " << g );
475  F /= f; G /= g;
476  TIMING_END_AND_PRINT (ez_content, "time to extract content in EZ: ")
477  if ( F.isUnivariate() && G.isUnivariate() )
478  {
479    DEBDECLEVEL( cerr, "ezgcd" );
480    if(F.mvar()==G.mvar())
481      d*=gcd(F,G);
482    else
483      return N (d);
484    return N (d);
485  }
486  if ( F.isUnivariate())
487  {
488    g= content (G,G.mvar());
489    return N(d*gcd(F,g));
490  }
491  if ( G.isUnivariate())
492  {
493    f= content (F,F.mvar());
494    return N(d*gcd(G,f));
495  }
496
497  int maxNumVars= tmax (getNumVars (F), getNumVars (G));
498  int sizeF= size (F);
499  int sizeG= size (G);
500
501
502  if (!isRat)
503    On (SW_RATIONAL);
504  if (sizeF/maxNumVars > 500 && sizeG/maxNumVars > 500)
505  {
506    Off(SW_USE_EZGCD);
507    result=gcd( F, G );
508    On(SW_USE_EZGCD);
509    if (!isRat)
510      Off (SW_RATIONAL);
511    result /= icontent (result);
512    DEBDECLEVEL( cerr, "ezgcd" );
513    return N (d*result);
514  }
515
516  int dummy= 0;
517  if ( gcd_test_one( F, G, false, dummy ) )
518  {
519    DEBDECLEVEL( cerr, "ezgcd" );
520    if (!isRat)
521      Off (SW_RATIONAL);
522    return N (d);
523  }
524  lcF = LC( F, x ); lcG = LC( G, x );
525  lcD = gcd( lcF, lcG );
526  delta = 0;
527  degF = degree( F, x ); degG = degree( G, x );
528  t = tmax( F.level(), G.level() );
529  if ( ! internal )
530    b = REvaluation( 2, t, IntRandom( 25 ) );
531  while ( ! gcdfound )
532  {
533    /// ---> A2
534    DEBOUTLN( cerr, "search for evaluation, delta = " << delta );
535    DEBOUTLN( cerr, "F = " << F );
536    DEBOUTLN( cerr, "G = " << G );
537    TIMING_START (ez_eval)
538    if (!findeval( F, G, Fb, Gb, Db, b, delta, degF, degG, maxeval, count,
539                   o, 25, l))
540    {
541      Off(SW_USE_EZGCD);
542      result=gcd( F, G );
543      On(SW_USE_EZGCD);
544      if (!isRat)
545        Off (SW_RATIONAL);
546      DEBDECLEVEL( cerr, "ezgcd" );
547      result /= icontent (result);
548      return N (d*result);
549    }
550    TIMING_END_AND_PRINT (ez_eval, "time to find eval point in EZ1: ")
551    DEBOUTLN( cerr, "found evaluation b = " << b );
552    DEBOUTLN( cerr, "F(b) = " << Fb );
553    DEBOUTLN( cerr, "G(b) = " << Gb );
554    DEBOUTLN( cerr, "D(b) = " << Db );
555    delta = degree( Db );
556    /// ---> A3
557    if ( delta == 0 )
558    {
559      DEBDECLEVEL( cerr, "ezgcd" );
560      if (!isRat)
561        Off (SW_RATIONAL);
562      return N (d);
563    }
564    /// ---> A4
565    //deltaold = delta;
566    while ( 1 ) 
567    {
568      bt = b;
569      TIMING_START (ez_eval)
570      if (!findeval( F, G, Fbt, Gbt, Dbt, bt, delta, degF, degG, maxeval, count,
571                     o, 25,l ))
572      {
573        Off(SW_USE_EZGCD);
574        result=gcd( F, G );
575        On(SW_USE_EZGCD);
576        if (!isRat)
577          Off (SW_RATIONAL);
578        DEBDECLEVEL( cerr, "ezgcd" );
579        result /= icontent (result);
580        return N (d*result);
581      }
582      TIMING_END_AND_PRINT (ez_eval, "time to find eval point in EZ2: ")
583      int dd=degree( Dbt );
584      if ( dd /*degree( Dbt )*/ == 0 )
585      {
586        DEBDECLEVEL( cerr, "ezgcd" );
587        if (!isRat)
588          Off (SW_RATIONAL);
589        return N (d);
590      }
591      if ( dd /*degree( Dbt )*/ == delta )
592        break;
593      else  if ( dd /*degree( Dbt )*/ < delta )
594      {
595        delta = dd /*degree( Dbt )*/;
596        b = bt;
597        Db = Dbt; Fb = Fbt; Gb = Gbt;
598      }
599      DEBOUTLN( cerr, "now after A4, delta = " << delta );
600      /// ---> A5
601      if (delta == degF)
602      {
603        if (degF <= degG  && fdivides (F, G))
604        {
605          DEBDECLEVEL( cerr, "ezgcd" );
606          if (!isRat)
607            Off (SW_RATIONAL);
608          return N (d*F);
609        }
610        else
611          delta--;
612      }
613      else if (delta == degG)
614      {
615        if (degG <= degF && fdivides( G, F ))
616        {
617          DEBDECLEVEL( cerr, "ezgcd" );
618          if (!isRat)
619            Off (SW_RATIONAL);
620          return N (d*G);
621        }
622        else
623          delta--;
624      }
625    }
626    if ( delta != degF && delta != degG )
627    {
628      /// ---> A6
629      bool B_is_F;
630      CanonicalForm xxx1, xxx2;
631      CanonicalForm buf;
632      DD[1] = Fb / Db;
633      buf= Gb/Db;
634      xxx1 = gcd( DD[1], Db );
635      xxx2 = gcd( buf, Db );
636      if (((xxx1.inCoeffDomain() && xxx2.inCoeffDomain()) &&
637          (size (F) <= size (G)))
638            || (xxx1.inCoeffDomain() && !xxx2.inCoeffDomain()))
639      {
640        B = F;
641        DD[2] = Db;
642        lcDD[1] = lcF;
643        lcDD[2] = lcD;
644        B_is_F = true;
645      }
646      else if (((xxx1.inCoeffDomain() && xxx2.inCoeffDomain()) &&
647                (size (G) < size (F)))
648                || (!xxx1.inCoeffDomain() && xxx2.inCoeffDomain()))
649      {
650        DD[1] = buf;
651        B = G;
652        DD[2] = Db;
653        lcDD[1] = lcG;
654        lcDD[2] = lcD;
655        B_is_F = false;
656      }
657      else
658      {
659        //special case
660        Off(SW_USE_EZGCD);
661        result=gcd( F, G );
662        On(SW_USE_EZGCD);
663        if (!isRat)
664          Off (SW_RATIONAL);
665        DEBDECLEVEL( cerr, "ezgcd" );
666        result /= icontent (result);
667        return N (d*result);
668      }
669      /// ---> A7
670      DD[2] = DD[2] * ( b( lcDD[2] ) / lc( DD[2] ) );
671      DD[1] = DD[1] * ( b( lcDD[1] ) / lc( DD[1] ) );
672      DEBOUTLN( cerr, "(hensel) B    = " << B );
673      DEBOUTLN( cerr, "(hensel) lcB  = " << LC( B, Variable(1) ) );
674      DEBOUTLN( cerr, "(hensel) b(B) = " << b(B) );
675      DEBOUTLN( cerr, "(hensel) DD   = " << DD );
676      DEBOUTLN( cerr, "(hensel) lcDD = " << lcDD );
677      TIMING_START (ez_hensel_lift)
678      gcdfound= Hensel (B*lcD, DD, b, lcDD);
679      TIMING_END_AND_PRINT (ez_hensel_lift, "time to hensel lift in EZ: ")
680      DEBOUTLN( cerr, "(hensel finished) DD   = " << DD );
681
682      if (gcdfound == -1)
683      {
684        Off (SW_USE_EZGCD);
685        result= gcd (F,G);
686        On (SW_USE_EZGCD);
687        if (!isRat)
688          Off (SW_RATIONAL);
689        DEBDECLEVEL( cerr, "ezgcd" );
690        result /= icontent (result);
691        return N (d*result);
692      }
693
694      if (gcdfound)
695      {
696        TIMING_START (ez_termination)
697        contcand= content (DD[2], Variable (1));
698        cand = DD[2] / contcand;
699        if (B_is_F)
700          gcdfound = fdivides( cand, G ) && cand*(DD[1]/(lcD/contcand)) == F;
701        else
702          gcdfound = fdivides( cand, F ) && cand*(DD[1]/(lcD/contcand)) == G;
703        TIMING_END_AND_PRINT (ez_termination,
704                              "time for termination test in EZ: ")
705      }
706      /// ---> A8 (gcdfound)
707    }
708    delta--;
709  }
710  /// ---> A9
711  DEBDECLEVEL( cerr, "ezgcd" );
712  cand *= bCommonDen (cand);
713  if (!isRat)
714    Off (SW_RATIONAL);
715  cand /= icontent (cand);
716  return N (d*cand);
717}
718#endif
719
720CanonicalForm
721ezgcd ( const CanonicalForm & FF, const CanonicalForm & GG )
722{
723#ifdef HAVE_NTL
724  REvaluation b;
725  return ezgcd( FF, GG, b, false );
726#else
727  Off (SW_USE_EZGCD);
728  return gcd (FF, GG);
729  On (SW_USE_EZGCD);
730#endif
731}
732
Note: See TracBrowser for help on using the repository browser.