source: git/factory/algext.cc @ 9fd8d2

spielwiese
Last change on this file since 9fd8d2 was 9fd8d2, checked in by Martin Lee <martinlee84@…>, 10 years ago
chg: deleted unused tryCRA and tryExtgcd
  • Property mode set to 100644
File size: 27.9 KB
Line 
1#ifdef HAVE_CONFIG_H
2#include "config.h"
3#endif /* HAVE_CONFIG_H */
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 "algext.h"
28#include "cf_map.h"
29#include "cf_generator.h"
30#include "facMul.h"
31#include "facNTLzzpEXGCD.h"
32
33#ifdef HAVE_NTL
34#include "NTLconvert.h"
35#endif
36
37#ifdef HAVE_FLINT
38#include "FLINTconvert.h"
39#endif
40
41TIMING_DEFINE_PRINT(alg_content_p)
42TIMING_DEFINE_PRINT(alg_content)
43TIMING_DEFINE_PRINT(alg_compress)
44TIMING_DEFINE_PRINT(alg_termination)
45TIMING_DEFINE_PRINT(alg_termination_p)
46TIMING_DEFINE_PRINT(alg_reconstruction)
47TIMING_DEFINE_PRINT(alg_newton_p)
48TIMING_DEFINE_PRINT(alg_recursion_p)
49TIMING_DEFINE_PRINT(alg_gcd_p)
50TIMING_DEFINE_PRINT(alg_euclid_p)
51
52/// compressing two polynomials F and G, M is used for compressing,
53/// N to reverse the compression
54static
55int 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
350bool hasFirstAlgVar( const CanonicalForm & f, Variable & a )
351{
352  if( f.inBaseDomain() ) // f has NO alg. variable
353    return false;
354  if( f.level()<0 ) // f has only alg. vars, so take the first one
355  {
356    a = f.mvar();
357    return true;
358  }
359  for(CFIterator i=f; i.hasTerms(); i++)
360    if( hasFirstAlgVar( i.coeff(), a ))
361      return true; // 'a' is already set
362  return false;
363}
364
365CanonicalForm QGCD( const CanonicalForm & F, const CanonicalForm & G );
366int * leadDeg(const CanonicalForm & f, int *degs);
367bool isLess(int *a, int *b, int lower, int upper);
368bool isEqual(int *a, int *b, int lower, int upper);
369CanonicalForm firstLC(const CanonicalForm & f);
370static CanonicalForm trycontent ( const CanonicalForm & f, const Variable & x, const CanonicalForm & M, bool & fail );
371static CanonicalForm tryvcontent ( const CanonicalForm & f, const Variable & x, const CanonicalForm & M, bool & fail );
372static CanonicalForm trycf_content ( const CanonicalForm & f, const CanonicalForm & g, const CanonicalForm & M, bool & fail );
373
374static inline CanonicalForm
375tryNewtonInterp (const CanonicalForm alpha, const CanonicalForm u,
376              const CanonicalForm newtonPoly, const CanonicalForm oldInterPoly,
377              const Variable & x, const CanonicalForm& M, bool& fail)
378{
379  CanonicalForm interPoly;
380
381  CanonicalForm inv;
382  tryInvert (newtonPoly (alpha, x), M, inv, fail);
383  if (fail)
384    return 0;
385
386  interPoly= oldInterPoly+reduce ((u - oldInterPoly (alpha, x))*inv*newtonPoly, M);
387  return interPoly;
388}
389
390void tryBrownGCD( const CanonicalForm & F, const CanonicalForm & G, const CanonicalForm & M, CanonicalForm & result, bool & fail, bool topLevel )
391{ // assume F,G are multivariate polys over Z/p(a) for big prime p, M "univariate" polynomial in an algebraic variable
392  // M is assumed to be monic
393  if(F.isZero())
394  {
395    if(G.isZero())
396    {
397      result = G; // G is zero
398      return;
399    }
400    if(G.inCoeffDomain())
401    {
402      tryInvert(G,M,result,fail);
403      if(fail)
404        return;
405      result = 1;
406      return;
407    }
408    // try to make G monic modulo M
409    CanonicalForm inv;
410    tryInvert(Lc(G),M,inv,fail);
411    if(fail)
412      return;
413    result = inv*G;
414    result= reduce (result, M);
415    return;
416  }
417  if(G.isZero()) // F is non-zero
418  {
419    if(F.inCoeffDomain())
420    {
421      tryInvert(F,M,result,fail);
422      if(fail)
423        return;
424      result = 1;
425      return;
426    }
427    // try to make F monic modulo M
428    CanonicalForm inv;
429    tryInvert(Lc(F),M,inv,fail);
430    if(fail)
431      return;
432    result = inv*F;
433    result= reduce (result, M);
434    return;
435  }
436  // here: F,G both nonzero
437  if(F.inCoeffDomain())
438  {
439    tryInvert(F,M,result,fail);
440    if(fail)
441      return;
442    result = 1;
443    return;
444  }
445  if(G.inCoeffDomain())
446  {
447    tryInvert(G,M,result,fail);
448    if(fail)
449      return;
450    result = 1;
451    return;
452  }
453  TIMING_START (alg_compress)
454  CFMap MM,NN;
455  int lev= myCompress (F, G, MM, NN, topLevel);
456  if (lev == 0)
457  {
458    result= 1;
459    return;
460  }
461  CanonicalForm f=MM(F);
462  CanonicalForm g=MM(G);
463  TIMING_END_AND_PRINT (alg_compress, "time to compress in alg gcd: ")
464  // here: f,g are compressed
465  // compute largest variable in f or g (least one is Variable(1))
466  int mv = f.level();
467  if(g.level() > mv)
468    mv = g.level();
469  // here: mv is level of the largest variable in f, g
470  Variable v1= Variable (1);
471#ifdef HAVE_NTL
472  Variable v= M.mvar();
473  if (fac_NTL_char != getCharacteristic())
474  {
475    fac_NTL_char= getCharacteristic();
476    zz_p::init (getCharacteristic());
477  }
478  zz_pX NTLMipo= convertFacCF2NTLzzpX (M);
479  zz_pE::init (NTLMipo);
480  zz_pEX NTLResult;
481  zz_pEX NTLF;
482  zz_pEX NTLG;
483#endif
484  if(mv == 1) // f,g univariate
485  {
486    TIMING_START (alg_euclid_p)
487#ifdef HAVE_NTL
488    NTLF= convertFacCF2NTLzz_pEX (f, NTLMipo);
489    NTLG= convertFacCF2NTLzz_pEX (g, NTLMipo);
490    tryNTLGCD (NTLResult, NTLF, NTLG, fail);
491    if (fail)
492      return;
493    result= convertNTLzz_pEX2CF (NTLResult, f.mvar(), v);
494#else
495    tryEuclid(f,g,M,result,fail);
496    if(fail)
497      return;
498#endif
499    TIMING_END_AND_PRINT (alg_euclid_p, "time for euclidean alg mod p: ")
500    result= NN (reduce (result, M)); // do not forget to map back
501    return;
502  }
503  TIMING_START (alg_content_p)
504  // here: mv > 1
505  CanonicalForm cf = tryvcontent(f, Variable(2), M, fail); // cf is univariate poly in var(1)
506  if(fail)
507    return;
508  CanonicalForm cg = tryvcontent(g, Variable(2), M, fail);
509  if(fail)
510    return;
511  CanonicalForm c;
512#ifdef HAVE_NTL
513  NTLF= convertFacCF2NTLzz_pEX (cf, NTLMipo);
514  NTLG= convertFacCF2NTLzz_pEX (cg, NTLMipo);
515  tryNTLGCD (NTLResult, NTLF, NTLG, fail);
516  if (fail)
517    return;
518  c= convertNTLzz_pEX2CF (NTLResult, v1, v);
519#else
520  tryEuclid(cf,cg,M,c,fail);
521  if(fail)
522    return;
523#endif
524  // f /= cf
525  f.tryDiv (cf, M, fail);
526  if(fail)
527    return;
528  // g /= cg
529  g.tryDiv (cg, M, fail);
530  if(fail)
531    return;
532  TIMING_END_AND_PRINT (alg_content_p, "time for content in alg gcd mod p: ")
533  if(f.inCoeffDomain())
534  {
535    tryInvert(f,M,result,fail);
536    if(fail)
537      return;
538    result = NN(c);
539    return;
540  }
541  if(g.inCoeffDomain())
542  {
543    tryInvert(g,M,result,fail);
544    if(fail)
545      return;
546    result = NN(c);
547    return;
548  }
549  int *L = new int[mv+1]; // L is addressed by i from 2 to mv
550  int *N = new int[mv+1];
551  for(int i=2; i<=mv; i++)
552    L[i] = N[i] = 0;
553  L = leadDeg(f, L);
554  N = leadDeg(g, N);
555  CanonicalForm gamma;
556  TIMING_START (alg_euclid_p)
557#ifdef HAVE_NTL
558  NTLF= convertFacCF2NTLzz_pEX (firstLC (f), NTLMipo);
559  NTLG= convertFacCF2NTLzz_pEX (firstLC (g), NTLMipo);
560  tryNTLGCD (NTLResult, NTLF, NTLG, fail);
561  if (fail)
562    return;
563  gamma= convertNTLzz_pEX2CF (NTLResult, v1, v);
564#else
565  tryEuclid( firstLC(f), firstLC(g), M, gamma, fail );
566  if(fail)
567    return;
568#endif
569  TIMING_END_AND_PRINT (alg_euclid_p, "time for gcd of lcs in alg mod p: ")
570  for(int i=2; i<=mv; i++) // entries at i=0,1 not visited
571    if(N[i] < L[i])
572      L[i] = N[i];
573  // L is now upper bound for degrees of gcd
574  int *dg_im = new int[mv+1]; // for the degree vector of the image we don't need any entry at i=1
575  for(int i=2; i<=mv; i++)
576    dg_im[i] = 0; // initialize
577  CanonicalForm gamma_image, m=1;
578  CanonicalForm gm=0;
579  CanonicalForm g_image, alpha, gnew;
580  FFGenerator gen = FFGenerator();
581  Variable x= Variable (1);
582  bool divides= true;
583  for(FFGenerator gen = FFGenerator(); gen.hasItems(); gen.next())
584  {
585    alpha = gen.item();
586    gamma_image = reduce(gamma(alpha, x),M); // plug in alpha for var(1)
587    if(gamma_image.isZero()) // skip lc-bad points var(1)-alpha
588      continue;
589    TIMING_START (alg_recursion_p)
590    tryBrownGCD( f(alpha, x), g(alpha, x), M, g_image, fail, false ); // recursive call with one var less
591    TIMING_END_AND_PRINT (alg_recursion_p,
592                         "time for recursive calls in alg gcd mod p: ")
593    if(fail)
594      return;
595    g_image = reduce(g_image, M);
596    if(g_image.inCoeffDomain()) // early termination
597    {
598      tryInvert(g_image,M,result,fail);
599      if(fail)
600        return;
601      result = NN(c);
602      return;
603    }
604    for(int i=2; i<=mv; i++)
605      dg_im[i] = 0; // reset (this is necessary, because some entries may not be updated by call to leadDeg)
606    dg_im = leadDeg(g_image, dg_im); // dg_im cannot be NIL-pointer
607    if(isEqual(dg_im, L, 2, mv))
608    {
609      CanonicalForm inv;
610      tryInvert (firstLC (g_image), M, inv, fail);
611      if (fail)
612        return;
613      g_image *= inv;
614      g_image *= gamma_image; // multiply by multiple of image lc(gcd)
615      g_image= reduce (g_image, M);
616      TIMING_START (alg_newton_p)
617      gnew= tryNewtonInterp (alpha, g_image, m, gm, x, M, fail);
618      TIMING_END_AND_PRINT (alg_newton_p,
619                            "time for Newton interpolation in alg gcd mod p: ")
620      // gnew = gm mod m
621      // gnew = g_image mod var(1)-alpha
622      // mnew = m * (var(1)-alpha)
623      if(fail)
624        return;
625      m *= (x - alpha);
626      if((firstLC(gnew) == gamma) || (gnew == gm)) // gnew did not change
627      {
628        TIMING_START (alg_termination_p)
629        cf = tryvcontent(gnew, Variable(2), M, fail);
630        if(fail)
631          return;
632        divides = true;
633        g_image= gnew;
634        g_image.tryDiv (cf, M, fail);
635        if(fail)
636          return;
637        divides= tryFdivides (g_image,f, M, fail); // trial division (f)
638        if(fail)
639          return;
640        if(divides)
641        {
642          bool divides2= tryFdivides (g_image,g, M, fail); // trial division (g)
643          if(fail)
644            return;
645          if(divides2)
646          {
647            result = NN(reduce (c*g_image, M));
648            TIMING_END_AND_PRINT (alg_termination_p,
649                      "time for successful termination test in alg gcd mod p: ")
650            return;
651          }
652        }
653        TIMING_END_AND_PRINT (alg_termination_p,
654                    "time for unsuccessful termination test in alg gcd mod p: ")
655      }
656      gm = gnew;
657      continue;
658    }
659
660    if(isLess(L, dg_im, 2, mv)) // dg_im > L --> current point unlucky
661      continue;
662
663    // here: isLess(dg_im, L, 2, mv) --> all previous points were unlucky
664    m = CanonicalForm(1); // reset
665    gm = 0; // reset
666    for(int i=2; i<=mv; i++) // tighten bound
667      L[i] = dg_im[i];
668  }
669  // we are out of evaluation points
670  fail = true;
671}
672
673static CanonicalForm
674myicontent ( const CanonicalForm & f, const CanonicalForm & c )
675{
676#ifdef HAVE_NTL
677    if (f.isOne() || c.isOne())
678      return 1;
679    if ( f.inBaseDomain() && c.inBaseDomain())
680    {
681      if (c.isZero()) return abs(f);
682      return bgcd( f, c );
683    }
684    else if ( (f.inCoeffDomain() && c.inCoeffDomain()) ||
685              (f.inCoeffDomain() && c.inBaseDomain()) ||
686              (f.inBaseDomain() && c.inCoeffDomain()))
687    {
688      if (c.isZero()) return abs (f);
689#ifdef HAVE_FLINT
690      fmpz_poly_t FLINTf, FLINTc;
691      convertFacCF2Fmpz_poly_t (FLINTf, f);
692      convertFacCF2Fmpz_poly_t (FLINTc, c);
693      fmpz_poly_gcd (FLINTc, FLINTc, FLINTf);
694      CanonicalForm result;
695      if (f.inCoeffDomain())
696        result= convertFmpz_poly_t2FacCF (FLINTc, f.mvar());
697      else
698        result= convertFmpz_poly_t2FacCF (FLINTc, c.mvar());
699      fmpz_poly_clear (FLINTc);
700      fmpz_poly_clear (FLINTf);
701      return result;
702#else
703      ZZX NTLf= convertFacCF2NTLZZX (f);
704      ZZX NTLc= convertFacCF2NTLZZX (c);
705      NTLc= GCD (NTLc, NTLf);
706      if (f.inCoeffDomain())
707        return convertNTLZZX2CF(NTLc,f.mvar());
708      else
709        return convertNTLZZX2CF(NTLc,c.mvar());
710#endif
711    }
712    else
713    {
714        CanonicalForm g = c;
715        for ( CFIterator i = f; i.hasTerms() && ! g.isOne(); i++ )
716            g = myicontent( i.coeff(), g );
717        return g;
718    }
719#else
720    return 1;
721#endif
722}
723
724CanonicalForm
725myicontent ( const CanonicalForm & f )
726{
727#ifdef HAVE_NTL
728    return myicontent( f, 0 );
729#else
730    return 1;
731#endif
732}
733
734CanonicalForm QGCD( const CanonicalForm & F, const CanonicalForm & G )
735{ // f,g in Q(a)[x1,...,xn]
736  if(F.isZero())
737  {
738    if(G.isZero())
739      return G; // G is zero
740    if(G.inCoeffDomain())
741      return CanonicalForm(1);
742    CanonicalForm lcinv= 1/Lc (G);
743    return G*lcinv; // return monic G
744  }
745  if(G.isZero()) // F is non-zero
746  {
747    if(F.inCoeffDomain())
748      return CanonicalForm(1);
749    CanonicalForm lcinv= 1/Lc (F);
750    return F*lcinv; // return monic F
751  }
752  if(F.inCoeffDomain() || G.inCoeffDomain())
753    return CanonicalForm(1);
754  // here: both NOT inCoeffDomain
755  CanonicalForm f, g, tmp, M, q, D, Dp, cl, newq, mipo;
756  int p, i;
757  int *bound, *other; // degree vectors
758  bool fail;
759  bool off_rational=!isOn(SW_RATIONAL);
760  On( SW_RATIONAL ); // needed by bCommonDen
761  f = F * bCommonDen(F);
762  g = G * bCommonDen(G);
763  TIMING_START (alg_content)
764  CanonicalForm contf= myicontent (f);
765  CanonicalForm contg= myicontent (g);
766  f /= contf;
767  g /= contg;
768  CanonicalForm gcdcfcg= myicontent (contf, contg);
769  TIMING_END_AND_PRINT (alg_content, "time for content in alg gcd: ")
770  Variable a, b;
771  if(hasFirstAlgVar(f,a))
772  {
773    if(hasFirstAlgVar(g,b))
774    {
775      if(b.level() > a.level())
776        a = b;
777    }
778  }
779  else
780  {
781    if(!hasFirstAlgVar(g,a))// both not in extension
782    {
783      Off( SW_RATIONAL );
784      Off( SW_USE_QGCD );
785      tmp = gcdcfcg*gcd( f, g );
786      On( SW_USE_QGCD );
787      if (off_rational) Off(SW_RATIONAL);
788      return tmp;
789    }
790  }
791  // here: a is the biggest alg. var in f and g AND some of f,g is in extension
792  setReduce(a,false); // do not reduce expressions modulo mipo
793  tmp = getMipo(a);
794  M = tmp * bCommonDen(tmp);
795  // here: f, g in Z[a][x1,...,xn], M in Z[a] not necessarily monic
796  Off( SW_RATIONAL ); // needed by mod
797  // calculate upper bound for degree vector of gcd
798  int mv = f.level(); i = g.level();
799  if(i > mv)
800    mv = i;
801  // here: mv is level of the largest variable in f, g
802  bound = new int[mv+1]; // 'bound' could be indexed from 0 to mv, but we will only use from 1 to mv
803  other = new int[mv+1];
804  for(int i=1; i<=mv; i++) // initialize 'bound', 'other' with zeros
805    bound[i] = other[i] = 0;
806  bound = leadDeg(f,bound); // 'bound' is set the leading degree vector of f
807  other = leadDeg(g,other);
808  for(int i=1; i<=mv; i++) // entry at i=0 not visited
809    if(other[i] < bound[i])
810      bound[i] = other[i];
811  // now 'bound' is the smaller vector
812  cl = lc(M) * lc(f) * lc(g);
813  q = 1;
814  D = 0;
815  CanonicalForm test= 0;
816  bool equal= false;
817  for( i=cf_getNumBigPrimes()-1; i>-1; i-- )
818  {
819    p = cf_getBigPrime(i);
820    if( mod( cl, p ).isZero() ) // skip lc-bad primes
821      continue;
822    fail = false;
823    setCharacteristic(p);
824    mipo = mapinto(M);
825    mipo /= mipo.lc();
826    // here: mipo is monic
827    TIMING_START (alg_gcd_p)
828    tryBrownGCD( mapinto(f), mapinto(g), mipo, Dp, fail );
829    TIMING_END_AND_PRINT (alg_gcd_p, "time for alg gcd mod p: ")
830    if( fail ) // mipo splits in char p
831      continue;
832    if( Dp.inCoeffDomain() ) // early termination
833    {
834      tryInvert(Dp,mipo,tmp,fail); // check if zero divisor
835      if(fail)
836        continue;
837      setReduce(a,true);
838      if (off_rational) Off(SW_RATIONAL); else On(SW_RATIONAL);
839      setCharacteristic(0);
840      return gcdcfcg;
841    }
842    setCharacteristic(0);
843    // here: Dp NOT inCoeffDomain
844    for(int i=1; i<=mv; i++)
845      other[i] = 0; // reset (this is necessary, because some entries may not be updated by call to leadDeg)
846    other = leadDeg(Dp,other);
847
848    if(isEqual(bound, other, 1, mv)) // equal
849    {
850      chineseRemainder( D, q, mapinto(Dp), p, tmp, newq );
851      // tmp = Dp mod p
852      // tmp = D mod q
853      // newq = p*q
854      q = newq;
855      if( D != tmp )
856        D = tmp;
857      On( SW_RATIONAL );
858      TIMING_START (alg_reconstruction)
859      tmp = Farey( D, q ); // Farey
860      tmp *= bCommonDen (tmp);
861      TIMING_END_AND_PRINT (alg_reconstruction,
862                            "time for rational reconstruction in alg gcd: ")
863      setReduce(a,true); // reduce expressions modulo mipo
864      On( SW_RATIONAL ); // needed by fdivides
865      if (test != tmp)
866        test= tmp;
867      else
868        equal= true; // modular image did not add any new information
869      TIMING_START (alg_termination)
870#ifdef HAVE_NTL
871#ifdef HAVE_FLINT
872      if (equal && tmp.isUnivariate() && f.isUnivariate() && g.isUnivariate()
873          && f.level() == tmp.level() && tmp.level() == g.level())
874      {
875        CanonicalForm Q, R, sf, sg, stmp;
876        Variable x= Variable (1);
877        sf= swapvar (f, f.mvar(), x);
878        sg= swapvar (g, f.mvar(), x);
879        stmp= swapvar (tmp, f.mvar(), x);
880        newtonDivrem (sf, stmp, Q, R);
881        if (R.isZero())
882        {
883          newtonDivrem (sg, stmp, Q, R);
884          if (R.isZero())
885          {
886            Off (SW_RATIONAL);
887            setReduce (a,true);
888            if (off_rational) Off(SW_RATIONAL); else On(SW_RATIONAL);
889            TIMING_END_AND_PRINT (alg_termination,
890                                 "time for successful termination test in alg gcd: ")
891            return tmp*gcdcfcg;
892          }
893        }
894      }
895      else
896#endif
897#endif
898      if(equal && fdivides( tmp, f ) && fdivides( tmp, g )) // trial division
899      {
900        Off( SW_RATIONAL );
901        setReduce(a,true);
902        if (off_rational) Off(SW_RATIONAL); else On(SW_RATIONAL);
903        TIMING_END_AND_PRINT (alg_termination,
904                            "time for successful termination test in alg gcd: ")
905        return tmp*gcdcfcg;
906      }
907      TIMING_END_AND_PRINT (alg_termination,
908                          "time for unsuccessful termination test in alg gcd: ")
909      Off( SW_RATIONAL );
910      setReduce(a,false); // do not reduce expressions modulo mipo
911      continue;
912    }
913    if( isLess(bound, other, 1, mv) ) // current prime unlucky
914      continue;
915    // here: isLess(other, bound, 1, mv) ) ==> all previous primes unlucky
916    q = p;
917    D = mapinto(Dp); // shortcut CRA // shortcut CRA
918    for(int i=1; i<=mv; i++) // tighten bound
919      bound[i] = other[i];
920  }
921  // hopefully, we never reach this point
922  setReduce(a,true);
923  Off( SW_USE_QGCD );
924  D = gcdcfcg*gcd( f, g );
925  On( SW_USE_QGCD );
926  if (off_rational) Off(SW_RATIONAL); else On(SW_RATIONAL);
927  return D;
928}
929
930
931int * leadDeg(const CanonicalForm & f, int *degs)
932{ // leading degree vector w.r.t. lex. monomial order x(i+1) > x(i)
933  // if f is in a coeff domain, the zero pointer is returned
934  // 'a' should point to an array of sufficient size level(f)+1
935  if(f.inCoeffDomain())
936    return 0;
937  CanonicalForm tmp = f;
938  do
939  {
940    degs[tmp.level()] = tmp.degree();
941    tmp = LC(tmp);
942  }
943  while(!tmp.inCoeffDomain());
944  return degs;
945}
946
947
948bool isLess(int *a, int *b, int lower, int upper)
949{ // compares the degree vectors a,b on the specified part. Note: x(i+1) > x(i)
950  for(int i=upper; i>=lower; i--)
951    if(a[i] == b[i])
952      continue;
953    else
954      return a[i] < b[i];
955  return true;
956}
957
958
959bool isEqual(int *a, int *b, int lower, int upper)
960{ // compares the degree vectors a,b on the specified part. Note: x(i+1) > x(i)
961  for(int i=lower; i<=upper; i++)
962    if(a[i] != b[i])
963      return false;
964  return true;
965}
966
967
968CanonicalForm firstLC(const CanonicalForm & f)
969{ // returns the leading coefficient (LC) of level <= 1
970  CanonicalForm ret = f;
971  while(ret.level() > 1)
972    ret = LC(ret);
973  return ret;
974}
975
976#ifndef HAVE_NTL
977void tryExtgcd( const CanonicalForm & F, const CanonicalForm & G, const CanonicalForm & M, CanonicalForm & result, CanonicalForm & s, CanonicalForm & t, bool & fail )
978{ // F, G are univariate polynomials (i.e. they have exactly one polynomial variable)
979  // F and G must have the same level AND level > 0
980  // we try to calculate gcd(F,G) = s*F + t*G
981  // if a zero divisor is encontered, 'fail' is set to one
982  // M is assumed to be monic
983  CanonicalForm P;
984  if(F.inCoeffDomain())
985  {
986    tryInvert( F, M, P, fail );
987    if(fail)
988      return;
989    result = 1;
990    s = P; t = 0;
991    return;
992  }
993  if(G.inCoeffDomain())
994  {
995    tryInvert( G, M, P, fail );
996    if(fail)
997      return;
998    result = 1;
999    s = 0; t = P;
1000    return;
1001  }
1002  // here: both not inCoeffDomain
1003  CanonicalForm inv, rem, tmp, u, v, q, sum=0;
1004  if( F.degree() > G.degree() )
1005  {
1006    P = F; result = G;  s=v=0; t=u=1;
1007  }
1008  else
1009  {
1010    P = G; result = F; s=v=1; t=u=0;
1011  }
1012  Variable x = P.mvar();
1013  // here: degree(P) >= degree(result)
1014  while(true)
1015  {
1016    tryDivrem (P, result, q, rem, inv, M, fail);
1017    if(fail)
1018      return;
1019    if( rem.isZero() )
1020    {
1021      s*=inv;
1022      s= reduce (s, M);
1023      t*=inv;
1024      t= reduce (t, M);
1025      result *= inv; // monify result
1026      result= reduce (result, M);
1027      return;
1028    }
1029    sum += q;
1030    if(result.degree(x) >= rem.degree(x))
1031    {
1032      P=result;
1033      result=rem;
1034      tmp=u-sum*s;
1035      u=s;
1036      s=tmp;
1037      tmp=v-sum*t;
1038      v=t;
1039      t=tmp;
1040      sum = 0; // reset
1041    }
1042    else
1043      P = rem;
1044  }
1045}
1046#endif
1047
1048static CanonicalForm trycontent ( const CanonicalForm & f, const Variable & x, const CanonicalForm & M, bool & fail )
1049{ // as 'content', but takes care of zero divisors
1050  ASSERT( x.level() > 0, "cannot calculate content with respect to algebraic variable" );
1051  Variable y = f.mvar();
1052  if ( y == x )
1053    return trycf_content( f, 0, M, fail );
1054  if ( y < x )
1055     return f;
1056  return swapvar( trycontent( swapvar( f, y, x ), y, M, fail ), y, x );
1057}
1058
1059
1060static CanonicalForm tryvcontent ( const CanonicalForm & f, const Variable & x, const CanonicalForm & M, bool & fail )
1061{ // as vcontent, but takes care of zero divisors
1062  ASSERT( x.level() > 0, "cannot calculate vcontent with respect to algebraic variable" );
1063  if ( f.mvar() <= x )
1064    return trycontent( f, x, M, fail );
1065  CFIterator i;
1066  CanonicalForm d = 0, e, ret;
1067  for ( i = f; i.hasTerms() && ! d.isOne() && ! fail; i++ )
1068  {
1069    e = tryvcontent( i.coeff(), x, M, fail );
1070    if(fail)
1071      break;
1072    tryBrownGCD( d, e, M, ret, fail );
1073    d = ret;
1074  }
1075  return d;
1076}
1077
1078
1079static CanonicalForm trycf_content ( const CanonicalForm & f, const CanonicalForm & g, const CanonicalForm & M, bool & fail )
1080{ // as cf_content, but takes care of zero divisors
1081  if ( f.inPolyDomain() || ( f.inExtension() && ! getReduce( f.mvar() ) ) )
1082  {
1083    CFIterator i = f;
1084    CanonicalForm tmp = g, result;
1085    while ( i.hasTerms() && ! tmp.isOne() && ! fail )
1086    {
1087      tryBrownGCD( i.coeff(), tmp, M, result, fail );
1088      tmp = result;
1089      i++;
1090    }
1091    return result;
1092  }
1093  return abs( f );
1094}
1095
Note: See TracBrowser for help on using the repository browser.