source: git/factory/fac_ezgcd.cc @ 9f84ad

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