source: git/factory/cfGcdAlgExt.cc @ eebdf2

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