source: git/factory/fac_ezgcd.cc @ 517530

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