source: git/factory/cfGcdAlgExt.cc @ d1a302

fieker-DuValspielwiese
Last change on this file since d1a302 was 88355eb, checked in by Martin Lee <martinlee84@…>, 10 years ago
chg: got rid of unnecessary swapping variables
  • Property mode set to 100644
File size: 27.3 KB
Line 
1
2#include "config.h"
3
4
5#ifndef NOSTREAMIO
6#ifdef HAVE_CSTDIO
7#include <cstdio>
8#else
9#include <stdio.h>
10#endif
11#ifdef HAVE_IOSTREAM_H
12#include <iostream.h>
13#elif defined(HAVE_IOSTREAM)
14#include <iostream>
15#endif
16#endif
17
18#include "cf_assert.h"
19#include "timing.h"
20
21#include "templates/ftmpl_functions.h"
22#include "cf_defs.h"
23#include "canonicalform.h"
24#include "cf_iter.h"
25#include "cf_primes.h"
26#include "cf_algorithm.h"
27#include "cfGcdAlgExt.h"
28#include "cfUnivarGcd.h"
29#include "cf_map.h"
30#include "cf_generator.h"
31#include "facMul.h"
32#include "cfNTLzzpEXGCD.h"
33
34#ifdef HAVE_NTL
35#include "NTLconvert.h"
36#endif
37
38#ifdef HAVE_FLINT
39#include "FLINTconvert.h"
40#endif
41
42TIMING_DEFINE_PRINT(alg_content_p)
43TIMING_DEFINE_PRINT(alg_content)
44TIMING_DEFINE_PRINT(alg_compress)
45TIMING_DEFINE_PRINT(alg_termination)
46TIMING_DEFINE_PRINT(alg_termination_p)
47TIMING_DEFINE_PRINT(alg_reconstruction)
48TIMING_DEFINE_PRINT(alg_newton_p)
49TIMING_DEFINE_PRINT(alg_recursion_p)
50TIMING_DEFINE_PRINT(alg_gcd_p)
51TIMING_DEFINE_PRINT(alg_euclid_p)
52
53/// compressing two polynomials F and G, M is used for compressing,
54/// N to reverse the compression
55static
56int myCompress (const CanonicalForm& F, const CanonicalForm& G, CFMap & M,
57                CFMap & N, bool topLevel)
58{
59  int n= tmax (F.level(), G.level());
60  int * degsf= new int [n + 1];
61  int * degsg= new int [n + 1];
62
63  for (int i = 0; i <= n; i++)
64    degsf[i]= degsg[i]= 0;
65
66  degsf= degrees (F, degsf);
67  degsg= degrees (G, degsg);
68
69  int both_non_zero= 0;
70  int f_zero= 0;
71  int g_zero= 0;
72  int both_zero= 0;
73
74  if (topLevel)
75  {
76    for (int i= 1; i <= n; i++)
77    {
78      if (degsf[i] != 0 && degsg[i] != 0)
79      {
80        both_non_zero++;
81        continue;
82      }
83      if (degsf[i] == 0 && degsg[i] != 0 && i <= G.level())
84      {
85        f_zero++;
86        continue;
87      }
88      if (degsg[i] == 0 && degsf[i] && i <= F.level())
89      {
90        g_zero++;
91        continue;
92      }
93    }
94
95    if (both_non_zero == 0)
96    {
97      delete [] degsf;
98      delete [] degsg;
99      return 0;
100    }
101
102    // map Variables which do not occur in both polynomials to higher levels
103    int k= 1;
104    int l= 1;
105    for (int i= 1; i <= n; i++)
106    {
107      if (degsf[i] != 0 && degsg[i] == 0 && i <= F.level())
108      {
109        if (k + both_non_zero != i)
110        {
111          M.newpair (Variable (i), Variable (k + both_non_zero));
112          N.newpair (Variable (k + both_non_zero), Variable (i));
113        }
114        k++;
115      }
116      if (degsf[i] == 0 && degsg[i] != 0 && i <= G.level())
117      {
118        if (l + g_zero + both_non_zero != i)
119        {
120          M.newpair (Variable (i), Variable (l + g_zero + both_non_zero));
121          N.newpair (Variable (l + g_zero + both_non_zero), Variable (i));
122        }
123        l++;
124      }
125    }
126
127    // sort Variables x_{i} in increasing order of
128    // min(deg_{x_{i}}(f),deg_{x_{i}}(g))
129    int m= tmax (F.level(), G.level());
130    int min_max_deg;
131    k= both_non_zero;
132    l= 0;
133    int i= 1;
134    while (k > 0)
135    {
136      if (degsf [i] != 0 && degsg [i] != 0)
137        min_max_deg= tmax (degsf[i], degsg[i]);
138      else
139        min_max_deg= 0;
140      while (min_max_deg == 0)
141      {
142        i++;
143        min_max_deg= tmax (degsf[i], degsg[i]);
144        if (degsf [i] != 0 && degsg [i] != 0)
145          min_max_deg= tmax (degsf[i], degsg[i]);
146        else
147          min_max_deg= 0;
148      }
149      for (int j= i + 1; j <=  m; j++)
150      {
151        if (tmax (degsf[j],degsg[j]) <= min_max_deg && degsf[j] != 0 && degsg [j] != 0)
152        {
153          min_max_deg= tmax (degsf[j], degsg[j]);
154          l= j;
155        }
156      }
157      if (l != 0)
158      {
159        if (l != k)
160        {
161          M.newpair (Variable (l), Variable(k));
162          N.newpair (Variable (k), Variable(l));
163          degsf[l]= 0;
164          degsg[l]= 0;
165          l= 0;
166        }
167        else
168        {
169          degsf[l]= 0;
170          degsg[l]= 0;
171          l= 0;
172        }
173      }
174      else if (l == 0)
175      {
176        if (i != k)
177        {
178          M.newpair (Variable (i), Variable (k));
179          N.newpair (Variable (k), Variable (i));
180          degsf[i]= 0;
181          degsg[i]= 0;
182        }
183        else
184        {
185          degsf[i]= 0;
186          degsg[i]= 0;
187        }
188        i++;
189      }
190      k--;
191    }
192  }
193  else
194  {
195    //arrange Variables such that no gaps occur
196    for (int i= 1; i <= n; i++)
197    {
198      if (degsf[i] == 0 && degsg[i] == 0)
199      {
200        both_zero++;
201        continue;
202      }
203      else
204      {
205        if (both_zero != 0)
206        {
207          M.newpair (Variable (i), Variable (i - both_zero));
208          N.newpair (Variable (i - both_zero), Variable (i));
209        }
210      }
211    }
212  }
213
214  delete [] degsf;
215  delete [] degsg;
216
217  return 1;
218}
219
220void tryInvert( const CanonicalForm & F, const CanonicalForm & M, CanonicalForm & inv, bool & fail )
221{ // F, M are required to be "univariate" polynomials in an algebraic variable
222  // we try to invert F modulo M
223  if(F.inBaseDomain())
224  {
225    if(F.isZero())
226    {
227      fail = true;
228      return;
229    }
230    inv = 1/F;
231    return;
232  }
233  CanonicalForm b;
234  Variable a = M.mvar();
235  Variable x = Variable(1);
236  if(!extgcd( replacevar( F, a, x ), replacevar( M, a, x ), inv, b ).isOne())
237    fail = true;
238  else
239    inv = replacevar( inv, x, a ); // change back to alg var
240}
241
242#ifndef HAVE_NTL
243void tryDivrem (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
244                CanonicalForm& R, CanonicalForm& inv, const CanonicalForm& mipo,
245                bool& fail)
246{
247  if (F.inCoeffDomain())
248  {
249    Q= 0;
250    R= F;
251    return;
252  }
253
254  CanonicalForm A, B;
255  Variable x= F.mvar();
256  A= F;
257  B= G;
258  int degA= degree (A, x);
259  int degB= degree (B, x);
260
261  if (degA < degB)
262  {
263    R= A;
264    Q= 0;
265    return;
266  }
267
268  tryInvert (Lc (B), mipo, inv, fail);
269  if (fail)
270    return;
271
272  R= A;
273  Q= 0;
274  CanonicalForm Qi;
275  for (int i= degA -degB; i >= 0; i--)
276  {
277    if (degree (R, x) == i + degB)
278    {
279      Qi= Lc (R)*inv*power (x, i);
280      Qi= reduce (Qi, mipo);
281      R -= Qi*B;
282      R= reduce (R, mipo);
283      Q += Qi;
284    }
285  }
286}
287
288void tryEuclid( const CanonicalForm & A, const CanonicalForm & B, const CanonicalForm & M, CanonicalForm & result, bool & fail )
289{
290  CanonicalForm P;
291  if(A.inCoeffDomain())
292  {
293    tryInvert( A, M, P, fail );
294    if(fail)
295      return;
296    result = 1;
297    return;
298  }
299  if(B.inCoeffDomain())
300  {
301    tryInvert( B, M, P, fail );
302    if(fail)
303      return;
304    result = 1;
305    return;
306  }
307  // here: both not inCoeffDomain
308  if( A.degree() > B.degree() )
309  {
310    P = A; result = B;
311  }
312  else
313  {
314    P = B; result = A;
315  }
316  CanonicalForm inv;
317  if( result.isZero() )
318  {
319    tryInvert( Lc(P), M, inv, fail );
320    if(fail)
321      return;
322    result = inv*P; // monify result (not reduced, yet)
323    result= reduce (result, M);
324    return;
325  }
326  Variable x = P.mvar();
327  CanonicalForm rem, Q;
328  // here: degree(P) >= degree(result)
329  while(true)
330  {
331    tryDivrem (P, result, Q, rem, inv, M, fail);
332    if (fail)
333      return;
334    if( rem.isZero() )
335    {
336      result *= inv;
337      result= reduce (result, M);
338      return;
339    }
340    if(result.degree(x) >= rem.degree(x))
341    {
342      P = result;
343      result = rem;
344    }
345    else
346      P = rem;
347  }
348}
349#endif
350
351CanonicalForm QGCD( const CanonicalForm & F, const CanonicalForm & G );
352int * leadDeg(const CanonicalForm & f, int *degs);
353bool isLess(int *a, int *b, int lower, int upper);
354bool isEqual(int *a, int *b, int lower, int upper);
355CanonicalForm firstLC(const CanonicalForm & f);
356static CanonicalForm trycontent ( const CanonicalForm & f, const Variable & x, const CanonicalForm & M, bool & fail );
357static CanonicalForm tryvcontent ( const CanonicalForm & f, const Variable & x, const CanonicalForm & M, bool & fail );
358static CanonicalForm trycf_content ( const CanonicalForm & f, const CanonicalForm & g, const CanonicalForm & M, bool & fail );
359
360static inline CanonicalForm
361tryNewtonInterp (const CanonicalForm & alpha, const CanonicalForm & u,
362              const CanonicalForm & newtonPoly, const CanonicalForm & oldInterPoly,
363              const Variable & x, const CanonicalForm& M, bool& fail)
364{
365  CanonicalForm interPoly;
366
367  CanonicalForm inv;
368  tryInvert (newtonPoly (alpha, x), M, inv, fail);
369  if (fail)
370    return 0;
371
372  interPoly= oldInterPoly+reduce ((u - oldInterPoly (alpha, x))*inv*newtonPoly, M);
373  return interPoly;
374}
375
376void tryBrownGCD( const CanonicalForm & F, const CanonicalForm & G, const CanonicalForm & M, CanonicalForm & result, bool & fail, bool topLevel )
377{ // assume F,G are multivariate polys over Z/p(a) for big prime p, M "univariate" polynomial in an algebraic variable
378  // M is assumed to be monic
379  if(F.isZero())
380  {
381    if(G.isZero())
382    {
383      result = G; // G is zero
384      return;
385    }
386    if(G.inCoeffDomain())
387    {
388      tryInvert(G,M,result,fail);
389      if(fail)
390        return;
391      result = 1;
392      return;
393    }
394    // try to make G monic modulo M
395    CanonicalForm inv;
396    tryInvert(Lc(G),M,inv,fail);
397    if(fail)
398      return;
399    result = inv*G;
400    result= reduce (result, M);
401    return;
402  }
403  if(G.isZero()) // F is non-zero
404  {
405    if(F.inCoeffDomain())
406    {
407      tryInvert(F,M,result,fail);
408      if(fail)
409        return;
410      result = 1;
411      return;
412    }
413    // try to make F monic modulo M
414    CanonicalForm inv;
415    tryInvert(Lc(F),M,inv,fail);
416    if(fail)
417      return;
418    result = inv*F;
419    result= reduce (result, M);
420    return;
421  }
422  // here: F,G both nonzero
423  if(F.inCoeffDomain())
424  {
425    tryInvert(F,M,result,fail);
426    if(fail)
427      return;
428    result = 1;
429    return;
430  }
431  if(G.inCoeffDomain())
432  {
433    tryInvert(G,M,result,fail);
434    if(fail)
435      return;
436    result = 1;
437    return;
438  }
439  TIMING_START (alg_compress)
440  CFMap MM,NN;
441  int lev= myCompress (F, G, MM, NN, topLevel);
442  if (lev == 0)
443  {
444    result= 1;
445    return;
446  }
447  CanonicalForm f=MM(F);
448  CanonicalForm g=MM(G);
449  TIMING_END_AND_PRINT (alg_compress, "time to compress in alg gcd: ")
450  // here: f,g are compressed
451  // compute largest variable in f or g (least one is Variable(1))
452  int mv = f.level();
453  if(g.level() > mv)
454    mv = g.level();
455  // here: mv is level of the largest variable in f, g
456  Variable v1= Variable (1);
457#ifdef HAVE_NTL
458  Variable v= M.mvar();
459  if (fac_NTL_char != getCharacteristic())
460  {
461    fac_NTL_char= getCharacteristic();
462    zz_p::init (getCharacteristic());
463  }
464  zz_pX NTLMipo= convertFacCF2NTLzzpX (M);
465  zz_pE::init (NTLMipo);
466  zz_pEX NTLResult;
467  zz_pEX NTLF;
468  zz_pEX NTLG;
469#endif
470  if(mv == 1) // f,g univariate
471  {
472    TIMING_START (alg_euclid_p)
473#ifdef HAVE_NTL
474    NTLF= convertFacCF2NTLzz_pEX (f, NTLMipo);
475    NTLG= convertFacCF2NTLzz_pEX (g, NTLMipo);
476    tryNTLGCD (NTLResult, NTLF, NTLG, fail);
477    if (fail)
478      return;
479    result= convertNTLzz_pEX2CF (NTLResult, f.mvar(), v);
480#else
481    tryEuclid(f,g,M,result,fail);
482    if(fail)
483      return;
484#endif
485    TIMING_END_AND_PRINT (alg_euclid_p, "time for euclidean alg mod p: ")
486    result= NN (reduce (result, M)); // do not forget to map back
487    return;
488  }
489  TIMING_START (alg_content_p)
490  // here: mv > 1
491  CanonicalForm cf = tryvcontent(f, Variable(2), M, fail); // cf is univariate poly in var(1)
492  if(fail)
493    return;
494  CanonicalForm cg = tryvcontent(g, Variable(2), M, fail);
495  if(fail)
496    return;
497  CanonicalForm c;
498#ifdef HAVE_NTL
499  NTLF= convertFacCF2NTLzz_pEX (cf, NTLMipo);
500  NTLG= convertFacCF2NTLzz_pEX (cg, NTLMipo);
501  tryNTLGCD (NTLResult, NTLF, NTLG, fail);
502  if (fail)
503    return;
504  c= convertNTLzz_pEX2CF (NTLResult, v1, v);
505#else
506  tryEuclid(cf,cg,M,c,fail);
507  if(fail)
508    return;
509#endif
510  // f /= cf
511  f.tryDiv (cf, M, fail);
512  if(fail)
513    return;
514  // g /= cg
515  g.tryDiv (cg, M, fail);
516  if(fail)
517    return;
518  TIMING_END_AND_PRINT (alg_content_p, "time for content in alg gcd mod p: ")
519  if(f.inCoeffDomain())
520  {
521    tryInvert(f,M,result,fail);
522    if(fail)
523      return;
524    result = NN(c);
525    return;
526  }
527  if(g.inCoeffDomain())
528  {
529    tryInvert(g,M,result,fail);
530    if(fail)
531      return;
532    result = NN(c);
533    return;
534  }
535  int *L = new int[mv+1]; // L is addressed by i from 2 to mv
536  int *N = new int[mv+1];
537  for(int i=2; i<=mv; i++)
538    L[i] = N[i] = 0;
539  L = leadDeg(f, L);
540  N = leadDeg(g, N);
541  CanonicalForm gamma;
542  TIMING_START (alg_euclid_p)
543#ifdef HAVE_NTL
544  NTLF= convertFacCF2NTLzz_pEX (firstLC (f), NTLMipo);
545  NTLG= convertFacCF2NTLzz_pEX (firstLC (g), NTLMipo);
546  tryNTLGCD (NTLResult, NTLF, NTLG, fail);
547  if (fail)
548    return;
549  gamma= convertNTLzz_pEX2CF (NTLResult, v1, v);
550#else
551  tryEuclid( firstLC(f), firstLC(g), M, gamma, fail );
552  if(fail)
553    return;
554#endif
555  TIMING_END_AND_PRINT (alg_euclid_p, "time for gcd of lcs in alg mod p: ")
556  for(int i=2; i<=mv; i++) // entries at i=0,1 not visited
557    if(N[i] < L[i])
558      L[i] = N[i];
559  // L is now upper bound for degrees of gcd
560  int *dg_im = new int[mv+1]; // for the degree vector of the image we don't need any entry at i=1
561  for(int i=2; i<=mv; i++)
562    dg_im[i] = 0; // initialize
563  CanonicalForm gamma_image, m=1;
564  CanonicalForm gm=0;
565  CanonicalForm g_image, alpha, gnew;
566  FFGenerator gen = FFGenerator();
567  Variable x= Variable (1);
568  bool divides= true;
569  for(FFGenerator gen = FFGenerator(); gen.hasItems(); gen.next())
570  {
571    alpha = gen.item();
572    gamma_image = reduce(gamma(alpha, x),M); // plug in alpha for var(1)
573    if(gamma_image.isZero()) // skip lc-bad points var(1)-alpha
574      continue;
575    TIMING_START (alg_recursion_p)
576    tryBrownGCD( f(alpha, x), g(alpha, x), M, g_image, fail, false ); // recursive call with one var less
577    TIMING_END_AND_PRINT (alg_recursion_p,
578                         "time for recursive calls in alg gcd mod p: ")
579    if(fail)
580      return;
581    g_image = reduce(g_image, M);
582    if(g_image.inCoeffDomain()) // early termination
583    {
584      tryInvert(g_image,M,result,fail);
585      if(fail)
586        return;
587      result = NN(c);
588      return;
589    }
590    for(int i=2; i<=mv; i++)
591      dg_im[i] = 0; // reset (this is necessary, because some entries may not be updated by call to leadDeg)
592    dg_im = leadDeg(g_image, dg_im); // dg_im cannot be NIL-pointer
593    if(isEqual(dg_im, L, 2, mv))
594    {
595      CanonicalForm inv;
596      tryInvert (firstLC (g_image), M, inv, fail);
597      if (fail)
598        return;
599      g_image *= inv;
600      g_image *= gamma_image; // multiply by multiple of image lc(gcd)
601      g_image= reduce (g_image, M);
602      TIMING_START (alg_newton_p)
603      gnew= tryNewtonInterp (alpha, g_image, m, gm, x, M, fail);
604      TIMING_END_AND_PRINT (alg_newton_p,
605                            "time for Newton interpolation in alg gcd mod p: ")
606      // gnew = gm mod m
607      // gnew = g_image mod var(1)-alpha
608      // mnew = m * (var(1)-alpha)
609      if(fail)
610        return;
611      m *= (x - alpha);
612      if((firstLC(gnew) == gamma) || (gnew == gm)) // gnew did not change
613      {
614        TIMING_START (alg_termination_p)
615        cf = tryvcontent(gnew, Variable(2), M, fail);
616        if(fail)
617          return;
618        divides = true;
619        g_image= gnew;
620        g_image.tryDiv (cf, M, fail);
621        if(fail)
622          return;
623        divides= tryFdivides (g_image,f, M, fail); // trial division (f)
624        if(fail)
625          return;
626        if(divides)
627        {
628          bool divides2= tryFdivides (g_image,g, M, fail); // trial division (g)
629          if(fail)
630            return;
631          if(divides2)
632          {
633            result = NN(reduce (c*g_image, M));
634            TIMING_END_AND_PRINT (alg_termination_p,
635                      "time for successful termination test in alg gcd mod p: ")
636            return;
637          }
638        }
639        TIMING_END_AND_PRINT (alg_termination_p,
640                    "time for unsuccessful termination test in alg gcd mod p: ")
641      }
642      gm = gnew;
643      continue;
644    }
645
646    if(isLess(L, dg_im, 2, mv)) // dg_im > L --> current point unlucky
647      continue;
648
649    // here: isLess(dg_im, L, 2, mv) --> all previous points were unlucky
650    m = CanonicalForm(1); // reset
651    gm = 0; // reset
652    for(int i=2; i<=mv; i++) // tighten bound
653      L[i] = dg_im[i];
654  }
655  // we are out of evaluation points
656  fail = true;
657}
658
659static CanonicalForm
660myicontent ( const CanonicalForm & f, const CanonicalForm & c )
661{
662#ifdef HAVE_NTL
663    if (f.isOne() || c.isOne())
664      return 1;
665    if ( f.inBaseDomain() && c.inBaseDomain())
666    {
667      if (c.isZero()) return abs(f);
668      return bgcd( f, c );
669    }
670    else if ( (f.inCoeffDomain() && c.inCoeffDomain()) ||
671              (f.inCoeffDomain() && c.inBaseDomain()) ||
672              (f.inBaseDomain() && c.inCoeffDomain()))
673    {
674      if (c.isZero()) return abs (f);
675#ifdef HAVE_FLINT
676      fmpz_poly_t FLINTf, FLINTc;
677      convertFacCF2Fmpz_poly_t (FLINTf, f);
678      convertFacCF2Fmpz_poly_t (FLINTc, c);
679      fmpz_poly_gcd (FLINTc, FLINTc, FLINTf);
680      CanonicalForm result;
681      if (f.inCoeffDomain())
682        result= convertFmpz_poly_t2FacCF (FLINTc, f.mvar());
683      else
684        result= convertFmpz_poly_t2FacCF (FLINTc, c.mvar());
685      fmpz_poly_clear (FLINTc);
686      fmpz_poly_clear (FLINTf);
687      return result;
688#else
689      ZZX NTLf= convertFacCF2NTLZZX (f);
690      ZZX NTLc= convertFacCF2NTLZZX (c);
691      NTLc= GCD (NTLc, NTLf);
692      if (f.inCoeffDomain())
693        return convertNTLZZX2CF(NTLc,f.mvar());
694      else
695        return convertNTLZZX2CF(NTLc,c.mvar());
696#endif
697    }
698    else
699    {
700        CanonicalForm g = c;
701        for ( CFIterator i = f; i.hasTerms() && ! g.isOne(); i++ )
702            g = myicontent( i.coeff(), g );
703        return g;
704    }
705#else
706    return 1;
707#endif
708}
709
710CanonicalForm
711myicontent ( const CanonicalForm & f )
712{
713#ifdef HAVE_NTL
714    return myicontent( f, 0 );
715#else
716    return 1;
717#endif
718}
719
720CanonicalForm QGCD( const CanonicalForm & F, const CanonicalForm & G )
721{ // f,g in Q(a)[x1,...,xn]
722  if(F.isZero())
723  {
724    if(G.isZero())
725      return G; // G is zero
726    if(G.inCoeffDomain())
727      return CanonicalForm(1);
728    CanonicalForm lcinv= 1/Lc (G);
729    return G*lcinv; // return monic G
730  }
731  if(G.isZero()) // F is non-zero
732  {
733    if(F.inCoeffDomain())
734      return CanonicalForm(1);
735    CanonicalForm lcinv= 1/Lc (F);
736    return F*lcinv; // return monic F
737  }
738  if(F.inCoeffDomain() || G.inCoeffDomain())
739    return CanonicalForm(1);
740  // here: both NOT inCoeffDomain
741  CanonicalForm f, g, tmp, M, q, D, Dp, cl, newq, mipo;
742  int p, i;
743  int *bound, *other; // degree vectors
744  bool fail;
745  bool off_rational=!isOn(SW_RATIONAL);
746  On( SW_RATIONAL ); // needed by bCommonDen
747  f = F * bCommonDen(F);
748  g = G * bCommonDen(G);
749  TIMING_START (alg_content)
750  CanonicalForm contf= myicontent (f);
751  CanonicalForm contg= myicontent (g);
752  f /= contf;
753  g /= contg;
754  CanonicalForm gcdcfcg= myicontent (contf, contg);
755  TIMING_END_AND_PRINT (alg_content, "time for content in alg gcd: ")
756  Variable a, b;
757  if(hasFirstAlgVar(f,a))
758  {
759    if(hasFirstAlgVar(g,b))
760    {
761      if(b.level() > a.level())
762        a = b;
763    }
764  }
765  else
766  {
767    if(!hasFirstAlgVar(g,a))// both not in extension
768    {
769      Off( SW_RATIONAL );
770      Off( SW_USE_QGCD );
771      tmp = gcdcfcg*gcd( f, g );
772      On( SW_USE_QGCD );
773      if (off_rational) Off(SW_RATIONAL);
774      return tmp;
775    }
776  }
777  // here: a is the biggest alg. var in f and g AND some of f,g is in extension
778  setReduce(a,false); // do not reduce expressions modulo mipo
779  tmp = getMipo(a);
780  M = tmp * bCommonDen(tmp);
781  // here: f, g in Z[a][x1,...,xn], M in Z[a] not necessarily monic
782  Off( SW_RATIONAL ); // needed by mod
783  // calculate upper bound for degree vector of gcd
784  int mv = f.level(); i = g.level();
785  if(i > mv)
786    mv = i;
787  // here: mv is level of the largest variable in f, g
788  bound = new int[mv+1]; // 'bound' could be indexed from 0 to mv, but we will only use from 1 to mv
789  other = new int[mv+1];
790  for(int i=1; i<=mv; i++) // initialize 'bound', 'other' with zeros
791    bound[i] = other[i] = 0;
792  bound = leadDeg(f,bound); // 'bound' is set the leading degree vector of f
793  other = leadDeg(g,other);
794  for(int i=1; i<=mv; i++) // entry at i=0 not visited
795    if(other[i] < bound[i])
796      bound[i] = other[i];
797  // now 'bound' is the smaller vector
798  cl = lc(M) * lc(f) * lc(g);
799  q = 1;
800  D = 0;
801  CanonicalForm test= 0;
802  bool equal= false;
803  for( i=cf_getNumBigPrimes()-1; i>-1; i-- )
804  {
805    p = cf_getBigPrime(i);
806    if( mod( cl, p ).isZero() ) // skip lc-bad primes
807      continue;
808    fail = false;
809    setCharacteristic(p);
810    mipo = mapinto(M);
811    mipo /= mipo.lc();
812    // here: mipo is monic
813    TIMING_START (alg_gcd_p)
814    tryBrownGCD( mapinto(f), mapinto(g), mipo, Dp, fail );
815    TIMING_END_AND_PRINT (alg_gcd_p, "time for alg gcd mod p: ")
816    if( fail ) // mipo splits in char p
817      continue;
818    if( Dp.inCoeffDomain() ) // early termination
819    {
820      tryInvert(Dp,mipo,tmp,fail); // check if zero divisor
821      if(fail)
822        continue;
823      setReduce(a,true);
824      if (off_rational) Off(SW_RATIONAL); else On(SW_RATIONAL);
825      setCharacteristic(0);
826      return gcdcfcg;
827    }
828    setCharacteristic(0);
829    // here: Dp NOT inCoeffDomain
830    for(int i=1; i<=mv; i++)
831      other[i] = 0; // reset (this is necessary, because some entries may not be updated by call to leadDeg)
832    other = leadDeg(Dp,other);
833
834    if(isEqual(bound, other, 1, mv)) // equal
835    {
836      chineseRemainder( D, q, mapinto(Dp), p, tmp, newq );
837      // tmp = Dp mod p
838      // tmp = D mod q
839      // newq = p*q
840      q = newq;
841      if( D != tmp )
842        D = tmp;
843      On( SW_RATIONAL );
844      TIMING_START (alg_reconstruction)
845      tmp = Farey( D, q ); // Farey
846      tmp *= bCommonDen (tmp);
847      TIMING_END_AND_PRINT (alg_reconstruction,
848                            "time for rational reconstruction in alg gcd: ")
849      setReduce(a,true); // reduce expressions modulo mipo
850      On( SW_RATIONAL ); // needed by fdivides
851      if (test != tmp)
852        test= tmp;
853      else
854        equal= true; // modular image did not add any new information
855      TIMING_START (alg_termination)
856#ifdef HAVE_NTL
857#ifdef HAVE_FLINT
858      if (equal && tmp.isUnivariate() && f.isUnivariate() && g.isUnivariate()
859          && f.level() == tmp.level() && tmp.level() == g.level())
860      {
861        CanonicalForm Q, R;
862        newtonDivrem (f, tmp, Q, R);
863        if (R.isZero())
864        {
865          newtonDivrem (g, tmp, Q, R);
866          if (R.isZero())
867          {
868            Off (SW_RATIONAL);
869            setReduce (a,true);
870            if (off_rational) Off(SW_RATIONAL); else On(SW_RATIONAL);
871            TIMING_END_AND_PRINT (alg_termination,
872                                 "time for successful termination test in alg gcd: ")
873            return tmp*gcdcfcg;
874          }
875        }
876      }
877      else
878#endif
879#endif
880      if(equal && fdivides( tmp, f ) && fdivides( tmp, g )) // trial division
881      {
882        Off( SW_RATIONAL );
883        setReduce(a,true);
884        if (off_rational) Off(SW_RATIONAL); else On(SW_RATIONAL);
885        TIMING_END_AND_PRINT (alg_termination,
886                            "time for successful termination test in alg gcd: ")
887        return tmp*gcdcfcg;
888      }
889      TIMING_END_AND_PRINT (alg_termination,
890                          "time for unsuccessful termination test in alg gcd: ")
891      Off( SW_RATIONAL );
892      setReduce(a,false); // do not reduce expressions modulo mipo
893      continue;
894    }
895    if( isLess(bound, other, 1, mv) ) // current prime unlucky
896      continue;
897    // here: isLess(other, bound, 1, mv) ) ==> all previous primes unlucky
898    q = p;
899    D = mapinto(Dp); // shortcut CRA // shortcut CRA
900    for(int i=1; i<=mv; i++) // tighten bound
901      bound[i] = other[i];
902  }
903  // hopefully, we never reach this point
904  setReduce(a,true);
905  Off( SW_USE_QGCD );
906  D = gcdcfcg*gcd( f, g );
907  On( SW_USE_QGCD );
908  if (off_rational) Off(SW_RATIONAL); else On(SW_RATIONAL);
909  return D;
910}
911
912
913int * leadDeg(const CanonicalForm & f, int *degs)
914{ // leading degree vector w.r.t. lex. monomial order x(i+1) > x(i)
915  // if f is in a coeff domain, the zero pointer is returned
916  // 'a' should point to an array of sufficient size level(f)+1
917  if(f.inCoeffDomain())
918    return 0;
919  CanonicalForm tmp = f;
920  do
921  {
922    degs[tmp.level()] = tmp.degree();
923    tmp = LC(tmp);
924  }
925  while(!tmp.inCoeffDomain());
926  return degs;
927}
928
929
930bool isLess(int *a, int *b, int lower, int upper)
931{ // compares the degree vectors a,b on the specified part. Note: x(i+1) > x(i)
932  for(int i=upper; i>=lower; i--)
933    if(a[i] == b[i])
934      continue;
935    else
936      return a[i] < b[i];
937  return true;
938}
939
940
941bool isEqual(int *a, int *b, int lower, int upper)
942{ // compares the degree vectors a,b on the specified part. Note: x(i+1) > x(i)
943  for(int i=lower; i<=upper; i++)
944    if(a[i] != b[i])
945      return false;
946  return true;
947}
948
949
950CanonicalForm firstLC(const CanonicalForm & f)
951{ // returns the leading coefficient (LC) of level <= 1
952  CanonicalForm ret = f;
953  while(ret.level() > 1)
954    ret = LC(ret);
955  return ret;
956}
957
958#ifndef HAVE_NTL
959void tryExtgcd( const CanonicalForm & F, const CanonicalForm & G, const CanonicalForm & M, CanonicalForm & result, CanonicalForm & s, CanonicalForm & t, bool & fail )
960{ // F, G are univariate polynomials (i.e. they have exactly one polynomial variable)
961  // F and G must have the same level AND level > 0
962  // we try to calculate gcd(F,G) = s*F + t*G
963  // if a zero divisor is encountered, 'fail' is set to one
964  // M is assumed to be monic
965  CanonicalForm P;
966  if(F.inCoeffDomain())
967  {
968    tryInvert( F, M, P, fail );
969    if(fail)
970      return;
971    result = 1;
972    s = P; t = 0;
973    return;
974  }
975  if(G.inCoeffDomain())
976  {
977    tryInvert( G, M, P, fail );
978    if(fail)
979      return;
980    result = 1;
981    s = 0; t = P;
982    return;
983  }
984  // here: both not inCoeffDomain
985  CanonicalForm inv, rem, tmp, u, v, q, sum=0;
986  if( F.degree() > G.degree() )
987  {
988    P = F; result = G;  s=v=0; t=u=1;
989  }
990  else
991  {
992    P = G; result = F; s=v=1; t=u=0;
993  }
994  Variable x = P.mvar();
995  // here: degree(P) >= degree(result)
996  while(true)
997  {
998    tryDivrem (P, result, q, rem, inv, M, fail);
999    if(fail)
1000      return;
1001    if( rem.isZero() )
1002    {
1003      s*=inv;
1004      s= reduce (s, M);
1005      t*=inv;
1006      t= reduce (t, M);
1007      result *= inv; // monify result
1008      result= reduce (result, M);
1009      return;
1010    }
1011    sum += q;
1012    if(result.degree(x) >= rem.degree(x))
1013    {
1014      P=result;
1015      result=rem;
1016      tmp=u-sum*s;
1017      u=s;
1018      s=tmp;
1019      tmp=v-sum*t;
1020      v=t;
1021      t=tmp;
1022      sum = 0; // reset
1023    }
1024    else
1025      P = rem;
1026  }
1027}
1028#endif
1029
1030static CanonicalForm trycontent ( const CanonicalForm & f, const Variable & x, const CanonicalForm & M, bool & fail )
1031{ // as 'content', but takes care of zero divisors
1032  ASSERT( x.level() > 0, "cannot calculate content with respect to algebraic variable" );
1033  Variable y = f.mvar();
1034  if ( y == x )
1035    return trycf_content( f, 0, M, fail );
1036  if ( y < x )
1037     return f;
1038  return swapvar( trycontent( swapvar( f, y, x ), y, M, fail ), y, x );
1039}
1040
1041
1042static CanonicalForm tryvcontent ( const CanonicalForm & f, const Variable & x, const CanonicalForm & M, bool & fail )
1043{ // as vcontent, but takes care of zero divisors
1044  ASSERT( x.level() > 0, "cannot calculate vcontent with respect to algebraic variable" );
1045  if ( f.mvar() <= x )
1046    return trycontent( f, x, M, fail );
1047  CFIterator i;
1048  CanonicalForm d = 0, e, ret;
1049  for ( i = f; i.hasTerms() && ! d.isOne() && ! fail; i++ )
1050  {
1051    e = tryvcontent( i.coeff(), x, M, fail );
1052    if(fail)
1053      break;
1054    tryBrownGCD( d, e, M, ret, fail );
1055    d = ret;
1056  }
1057  return d;
1058}
1059
1060
1061static CanonicalForm trycf_content ( const CanonicalForm & f, const CanonicalForm & g, const CanonicalForm & M, bool & fail )
1062{ // as cf_content, but takes care of zero divisors
1063  if ( f.inPolyDomain() || ( f.inExtension() && ! getReduce( f.mvar() ) ) )
1064  {
1065    CFIterator i = f;
1066    CanonicalForm tmp = g, result;
1067    while ( i.hasTerms() && ! tmp.isOne() && ! fail )
1068    {
1069      tryBrownGCD( i.coeff(), tmp, M, result, fail );
1070      tmp = result;
1071      i++;
1072    }
1073    return result;
1074  }
1075  return abs( f );
1076}
1077
Note: See TracBrowser for help on using the repository browser.