source: git/factory/fac_ezgcd.cc @ ce41efa

spielwiese
Last change on this file since ce41efa was ce41efa, checked in by Martin Lee <martinlee84@…>, 11 years ago
chg: estimate term growth through shifting evaluation point in ezgcd fix: stupid bug in ezgcd
  • Property mode set to 100644
File size: 16.7 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    return N (d);
483  }
484
485  int maxNumVars= tmax (getNumVars (F), getNumVars (G));
486  int sizeF= size (F);
487  int sizeG= size (G);
488
489
490  if (!isRat)
491    On (SW_RATIONAL);
492  if (sizeF/maxNumVars > 500 && sizeG/maxNumVars > 500)
493  {
494    Off(SW_USE_EZGCD);
495    result=gcd( F, G );
496    On(SW_USE_EZGCD);
497    if (!isRat)
498      Off (SW_RATIONAL);
499    result /= icontent (result);
500    DEBDECLEVEL( cerr, "ezgcd" );
501    return N (d*result);
502  }
503
504  int dummy= 0;
505  if ( gcd_test_one( F, G, false, dummy ) )
506  {
507    DEBDECLEVEL( cerr, "ezgcd" );
508    if (!isRat)
509      Off (SW_RATIONAL);
510    return N (d);
511  }
512  lcF = LC( F, x ); lcG = LC( G, x );
513  lcD = gcd( lcF, lcG );
514  delta = 0;
515  degF = degree( F, x ); degG = degree( G, x );
516  t = tmax( F.level(), G.level() );
517  if ( ! internal )
518    b = REvaluation( 2, t, IntRandom( 25 ) );
519  while ( ! gcdfound )
520  {
521    /// ---> A2
522    DEBOUTLN( cerr, "search for evaluation, delta = " << delta );
523    DEBOUTLN( cerr, "F = " << F );
524    DEBOUTLN( cerr, "G = " << G );
525    TIMING_START (ez_eval)
526    if (!findeval( F, G, Fb, Gb, Db, b, delta, degF, degG, maxeval, count,
527                   o, 25, l))
528    {
529      Off(SW_USE_EZGCD);
530      result=gcd( F, G );
531      On(SW_USE_EZGCD);
532      if (!isRat)
533        Off (SW_RATIONAL);
534      DEBDECLEVEL( cerr, "ezgcd" );
535      result /= icontent (result);
536      return N (d*result);
537    }
538    TIMING_END_AND_PRINT (ez_eval, "time to find eval point in EZ1: ")
539    DEBOUTLN( cerr, "found evaluation b = " << b );
540    DEBOUTLN( cerr, "F(b) = " << Fb );
541    DEBOUTLN( cerr, "G(b) = " << Gb );
542    DEBOUTLN( cerr, "D(b) = " << Db );
543    delta = degree( Db );
544    /// ---> A3
545    if ( delta == 0 )
546    {
547      DEBDECLEVEL( cerr, "ezgcd" );
548      if (!isRat)
549        Off (SW_RATIONAL);
550      return N (d);
551    }
552    /// ---> A4
553    //deltaold = delta;
554    while ( 1 ) 
555    {
556      bt = b;
557      TIMING_START (ez_eval)
558      if (!findeval( F, G, Fbt, Gbt, Dbt, bt, delta, degF, degG, maxeval, count,
559                     o, 25,l ))
560      {
561        Off(SW_USE_EZGCD);
562        result=gcd( F, G );
563        On(SW_USE_EZGCD);
564        if (!isRat)
565          Off (SW_RATIONAL);
566        DEBDECLEVEL( cerr, "ezgcd" );
567        result /= icontent (result);
568        return N (d*result);
569      }
570      TIMING_END_AND_PRINT (ez_eval, "time to find eval point in EZ2: ")
571      int dd=degree( Dbt );
572      if ( dd /*degree( Dbt )*/ == 0 )
573      {
574        DEBDECLEVEL( cerr, "ezgcd" );
575        if (!isRat)
576          Off (SW_RATIONAL);
577        return N (d);
578      }
579      if ( dd /*degree( Dbt )*/ == delta )
580        break;
581      else  if ( dd /*degree( Dbt )*/ < delta )
582      {
583        delta = dd /*degree( Dbt )*/;
584        b = bt;
585        Db = Dbt; Fb = Fbt; Gb = Gbt;
586      }
587      DEBOUTLN( cerr, "now after A4, delta = " << delta );
588      /// ---> A5
589      if (delta == degF)
590      {
591        if (degF <= degG  && fdivides (F, G))
592        {
593          DEBDECLEVEL( cerr, "ezgcd" );
594          if (!isRat)
595            Off (SW_RATIONAL);
596          return N (d*F);
597        }
598        else
599          delta--;
600      }
601      else if (delta == degG)
602      {
603        if (degG <= degF && fdivides( G, F ))
604        {
605          DEBDECLEVEL( cerr, "ezgcd" );
606          if (!isRat)
607            Off (SW_RATIONAL);
608          return N (d*G);
609        }
610        else
611          delta--;
612      }
613    }
614    if ( delta != degF && delta != degG )
615    {
616      /// ---> A6
617      bool B_is_F;
618      CanonicalForm xxx1, xxx2;
619      CanonicalForm buf;
620      DD[1] = Fb / Db;
621      buf= Gb/Db;
622      xxx1 = gcd( DD[1], Db );
623      xxx2 = gcd( buf, Db );
624      if (((xxx1.inCoeffDomain() && xxx2.inCoeffDomain()) &&
625          (size (F) <= size (G)))
626            || (xxx1.inCoeffDomain() && !xxx2.inCoeffDomain()))
627      {
628        B = F;
629        DD[2] = Db;
630        lcDD[1] = lcF;
631        lcDD[2] = lcD;
632        B_is_F = true;
633      }
634      else if (((xxx1.inCoeffDomain() && xxx2.inCoeffDomain()) &&
635                (size (G) < size (F)))
636                || (!xxx1.inCoeffDomain() && xxx2.inCoeffDomain()))
637      {
638        DD[1] = buf;
639        B = G;
640        DD[2] = Db;
641        lcDD[1] = lcG;
642        lcDD[2] = lcD;
643        B_is_F = false;
644      }
645      else
646      {
647        //special case
648        Off(SW_USE_EZGCD);
649        result=gcd( F, G );
650        On(SW_USE_EZGCD);
651        if (!isRat)
652          Off (SW_RATIONAL);
653        DEBDECLEVEL( cerr, "ezgcd" );
654        result /= icontent (result);
655        return N (d*result);
656      }
657      /// ---> A7
658      DD[2] = DD[2] * ( b( lcDD[2] ) / lc( DD[2] ) );
659      DD[1] = DD[1] * ( b( lcDD[1] ) / lc( DD[1] ) );
660      DEBOUTLN( cerr, "(hensel) B    = " << B );
661      DEBOUTLN( cerr, "(hensel) lcB  = " << LC( B, Variable(1) ) );
662      DEBOUTLN( cerr, "(hensel) b(B) = " << b(B) );
663      DEBOUTLN( cerr, "(hensel) DD   = " << DD );
664      DEBOUTLN( cerr, "(hensel) lcDD = " << lcDD );
665      TIMING_START (ez_hensel_lift)
666      gcdfound= Hensel (B*lcD, DD, b, lcDD);
667      TIMING_END_AND_PRINT (ez_hensel_lift, "time to hensel lift in EZ: ")
668      DEBOUTLN( cerr, "(hensel finished) DD   = " << DD );
669
670      if (gcdfound == -1)
671      {
672        Off (SW_USE_EZGCD);
673        result= gcd (F,G);
674        On (SW_USE_EZGCD);
675        if (!isRat)
676          Off (SW_RATIONAL);
677        DEBDECLEVEL( cerr, "ezgcd" );
678        result /= icontent (result);
679        return N (d*result);
680      }
681
682      if (gcdfound)
683      {
684        TIMING_START (ez_termination)
685        contcand= content (DD[2], Variable (1));
686        cand = DD[2] / contcand;
687        if (B_is_F)
688          gcdfound = fdivides( cand, G ) && cand*(DD[1]/(lcD/contcand)) == F;
689        else
690          gcdfound = fdivides( cand, F ) && cand*(DD[1]/(lcD/contcand)) == G;
691        TIMING_END_AND_PRINT (ez_termination,
692                              "time for termination test in EZ: ")
693      }
694      /// ---> A8 (gcdfound)
695    }
696    delta--;
697  }
698  /// ---> A9
699  DEBDECLEVEL( cerr, "ezgcd" );
700  cand *= bCommonDen (cand);
701  if (!isRat)
702    Off (SW_RATIONAL);
703  cand /= icontent (cand);
704  return N (d*cand);
705}
706#endif
707
708CanonicalForm
709ezgcd ( const CanonicalForm & FF, const CanonicalForm & GG )
710{
711#ifdef HAVE_NTL
712  REvaluation b;
713  return ezgcd( FF, GG, b, false );
714#else
715  Off (SW_USE_EZGCD);
716  return gcd (FF, GG);
717  On (SW_USE_EZGCD);
718#endif
719}
720
Note: See TracBrowser for help on using the repository browser.