source: git/factory/algext.cc @ 16f511

spielwiese
Last change on this file since 16f511 was 16f511, checked in by Oleksandr Motsak <motsak@…>, 11 years ago
Fixed the usage of "config.h" (if defined HAVE_CONFIG_H)
  • Property mode set to 100644
File size: 32.0 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
241void tryDivrem (const CanonicalForm& F, const CanonicalForm& G, CanonicalForm& Q,
242                CanonicalForm& R, CanonicalForm& inv, const CanonicalForm& mipo,
243                bool& fail)
244{
245  if (F.inCoeffDomain())
246  {
247    Q= 0;
248    R= F;
249    return;
250  }
251
252  CanonicalForm A, B;
253  Variable x= F.mvar();
254  A= F;
255  B= G;
256  int degA= degree (A, x);
257  int degB= degree (B, x);
258
259  if (degA < degB)
260  {
261    R= A;
262    Q= 0;
263    return;
264  }
265
266  tryInvert (Lc (B), mipo, inv, fail);
267  if (fail)
268    return;
269
270  R= A;
271  Q= 0;
272  CanonicalForm Qi;
273  for (int i= degA -degB; i >= 0; i--)
274  {
275    if (degree (R, x) == i + degB)
276    {
277      Qi= Lc (R)*inv*power (x, i);
278      Qi= reduce (Qi, mipo);
279      R -= Qi*B;
280      R= reduce (R, mipo);
281      Q += Qi;
282    }
283  }
284}
285
286void tryEuclid( const CanonicalForm & A, const CanonicalForm & B, const CanonicalForm & M, CanonicalForm & result, bool & fail )
287{
288  CanonicalForm P;
289  if(A.inCoeffDomain())
290  {
291    tryInvert( A, M, P, fail );
292    if(fail)
293      return;
294    result = 1;
295    return;
296  }
297  if(B.inCoeffDomain())
298  {
299    tryInvert( B, M, P, fail );
300    if(fail)
301      return;
302    result = 1;
303    return;
304  }
305  // here: both not inCoeffDomain
306  if( A.degree() > B.degree() )
307  {
308    P = A; result = B;
309  }
310  else
311  {
312    P = B; result = A;
313  }
314  CanonicalForm inv;
315  if( result.isZero() )
316  {
317    tryInvert( Lc(P), M, inv, fail );
318    if(fail)
319      return;
320    result = inv*P; // monify result (not reduced, yet)
321    result= reduce (result, M);
322    return;
323  }
324  Variable x = P.mvar();
325  CanonicalForm rem, Q;
326  // here: degree(P) >= degree(result)
327  while(true)
328  {
329    tryDivrem (P, result, Q, rem, inv, M, fail);
330    if (fail)
331      return;
332    if( rem.isZero() )
333    {
334      result *= inv;
335      result= reduce (result, M);
336      return;
337    }
338    if(result.degree(x) >= rem.degree(x))
339    {
340      P = result;
341      result = rem;
342    }
343    else
344      P = rem;
345  }
346}
347
348bool hasFirstAlgVar( const CanonicalForm & f, Variable & a )
349{
350  if( f.inBaseDomain() ) // f has NO alg. variable
351    return false;
352  if( f.level()<0 ) // f has only alg. vars, so take the first one
353  {
354    a = f.mvar();
355    return true;
356  }
357  for(CFIterator i=f; i.hasTerms(); i++)
358    if( hasFirstAlgVar( i.coeff(), a ))
359      return true; // 'a' is already set
360  return false;
361}
362
363CanonicalForm QGCD( const CanonicalForm & F, const CanonicalForm & G );
364int * leadDeg(const CanonicalForm & f, int *degs);
365bool isLess(int *a, int *b, int lower, int upper);
366bool isEqual(int *a, int *b, int lower, int upper);
367CanonicalForm firstLC(const CanonicalForm & f);
368static CanonicalForm trycontent ( const CanonicalForm & f, const Variable & x, const CanonicalForm & M, bool & fail );
369static CanonicalForm tryvcontent ( const CanonicalForm & f, const Variable & x, const CanonicalForm & M, bool & fail );
370static CanonicalForm trycf_content ( const CanonicalForm & f, const CanonicalForm & g, const CanonicalForm & M, bool & fail );
371
372static inline CanonicalForm
373tryNewtonInterp (const CanonicalForm alpha, const CanonicalForm u,
374              const CanonicalForm newtonPoly, const CanonicalForm oldInterPoly,
375              const Variable & x, const CanonicalForm& M, bool& fail)
376{
377  CanonicalForm interPoly;
378
379  CanonicalForm inv;
380  tryInvert (newtonPoly (alpha, x), M, inv, fail);
381  if (fail)
382    return 0;
383
384  interPoly= oldInterPoly+reduce ((u - oldInterPoly (alpha, x))*inv*newtonPoly, M);
385  return interPoly;
386}
387
388void tryBrownGCD( const CanonicalForm & F, const CanonicalForm & G, const CanonicalForm & M, CanonicalForm & result, bool & fail, bool topLevel )
389{ // assume F,G are multivariate polys over Z/p(a) for big prime p, M "univariate" polynomial in an algebraic variable
390  // M is assumed to be monic
391  if(F.isZero())
392  {
393    if(G.isZero())
394    {
395      result = G; // G is zero
396      return;
397    }
398    if(G.inCoeffDomain())
399    {
400      tryInvert(G,M,result,fail);
401      if(fail)
402        return;
403      result = 1;
404      return;
405    }
406    // try to make G monic modulo M
407    CanonicalForm inv;
408    tryInvert(Lc(G),M,inv,fail);
409    if(fail)
410      return;
411    result = inv*G;
412    result= reduce (result, M);
413    return;
414  }
415  if(G.isZero()) // F is non-zero
416  {
417    if(F.inCoeffDomain())
418    {
419      tryInvert(F,M,result,fail);
420      if(fail)
421        return;
422      result = 1;
423      return;
424    }
425    // try to make F monic modulo M
426    CanonicalForm inv;
427    tryInvert(Lc(F),M,inv,fail);
428    if(fail)
429      return;
430    result = inv*F;
431    result= reduce (result, M);
432    return;
433  }
434  // here: F,G both nonzero
435  if(F.inCoeffDomain())
436  {
437    tryInvert(F,M,result,fail);
438    if(fail)
439      return;
440    result = 1;
441    return;
442  }
443  if(G.inCoeffDomain())
444  {
445    tryInvert(G,M,result,fail);
446    if(fail)
447      return;
448    result = 1;
449    return;
450  }
451  TIMING_START (alg_compress)
452  CFMap MM,NN;
453  int lev= myCompress (F, G, MM, NN, topLevel);
454  if (lev == 0)
455  {
456    result= 1;
457    return;
458  }
459  CanonicalForm f=MM(F);
460  CanonicalForm g=MM(G);
461  TIMING_END_AND_PRINT (alg_compress, "time to compress in alg gcd: ")
462  // here: f,g are compressed
463  // compute largest variable in f or g (least one is Variable(1))
464  int mv = f.level();
465  if(g.level() > mv)
466    mv = g.level();
467  // here: mv is level of the largest variable in f, g
468  Variable v1= Variable (1);
469#ifdef HAVE_NTL
470  Variable v= M.mvar();
471  if (fac_NTL_char != getCharacteristic())
472  {
473    fac_NTL_char= getCharacteristic();
474    zz_p::init (getCharacteristic());
475  }
476  zz_pX NTLMipo= convertFacCF2NTLzzpX (M);
477  zz_pE::init (NTLMipo);
478  zz_pEX NTLResult;
479  zz_pEX NTLF;
480  zz_pEX NTLG;
481#endif
482  if(mv == 1) // f,g univariate
483  {
484    TIMING_START (alg_euclid_p)
485#ifdef HAVE_NTL
486    NTLF= convertFacCF2NTLzz_pEX (f, NTLMipo);
487    NTLG= convertFacCF2NTLzz_pEX (g, NTLMipo);
488    tryNTLGCD (NTLResult, NTLF, NTLG, fail);
489    if (fail)
490      return;
491    result= convertNTLzz_pEX2CF (NTLResult, f.mvar(), v);
492#else
493    tryEuclid(f,g,M,result,fail);
494    if(fail)
495      return;
496#endif
497    TIMING_END_AND_PRINT (alg_euclid_p, "time for euclidean alg mod p: ")
498    result= NN (reduce (result, M)); // do not forget to map back
499    return;
500  }
501  TIMING_START (alg_content_p)
502  // here: mv > 1
503  CanonicalForm cf = tryvcontent(f, Variable(2), M, fail); // cf is univariate poly in var(1)
504  if(fail)
505    return;
506  CanonicalForm cg = tryvcontent(g, Variable(2), M, fail);
507  if(fail)
508    return;
509  CanonicalForm c;
510#ifdef HAVE_NTL
511  NTLF= convertFacCF2NTLzz_pEX (cf, NTLMipo);
512  NTLG= convertFacCF2NTLzz_pEX (cg, NTLMipo);
513  tryNTLGCD (NTLResult, NTLF, NTLG, fail);
514  if (fail)
515    return;
516  c= convertNTLzz_pEX2CF (NTLResult, v1, v);
517#else
518  tryEuclid(cf,cg,M,c,fail);
519  if(fail)
520    return;
521#endif
522  // f /= cf
523  f.tryDiv (cf, M, fail);
524  if(fail)
525    return;
526  // g /= cg
527  g.tryDiv (cg, M, fail);
528  if(fail)
529    return;
530  TIMING_END_AND_PRINT (alg_content_p, "time for content in alg gcd mod p: ")
531  if(f.inCoeffDomain())
532  {
533    tryInvert(f,M,result,fail);
534    if(fail)
535      return;
536    result = NN(c);
537    return;
538  }
539  if(g.inCoeffDomain())
540  {
541    tryInvert(g,M,result,fail);
542    if(fail)
543      return;
544    result = NN(c);
545    return;
546  }
547  int *L = new int[mv+1]; // L is addressed by i from 2 to mv
548  int *N = new int[mv+1];
549  for(int i=2; i<=mv; i++)
550    L[i] = N[i] = 0;
551  L = leadDeg(f, L);
552  N = leadDeg(g, N);
553  CanonicalForm gamma;
554  TIMING_START (alg_euclid_p)
555#ifdef HAVE_NTL
556  NTLF= convertFacCF2NTLzz_pEX (firstLC (f), NTLMipo);
557  NTLG= convertFacCF2NTLzz_pEX (firstLC (g), NTLMipo);
558  tryNTLGCD (NTLResult, NTLF, NTLG, fail);
559  if (fail)
560    return;
561  gamma= convertNTLzz_pEX2CF (NTLResult, v1, v);
562#else
563  tryEuclid( firstLC(f), firstLC(g), M, gamma, fail );
564  if(fail)
565    return;
566#endif
567  TIMING_END_AND_PRINT (alg_euclid_p, "time for gcd of lcs in alg mod p: ")
568  for(int i=2; i<=mv; i++) // entries at i=0,1 not visited
569    if(N[i] < L[i])
570      L[i] = N[i];
571  // L is now upper bound for degrees of gcd
572  int *dg_im = new int[mv+1]; // for the degree vector of the image we don't need any entry at i=1
573  for(int i=2; i<=mv; i++)
574    dg_im[i] = 0; // initialize
575  CanonicalForm gamma_image, m=1;
576  CanonicalForm gm=0;
577  CanonicalForm g_image, alpha, gnew;
578  FFGenerator gen = FFGenerator();
579  Variable x= Variable (1);
580  bool divides= true;
581  for(FFGenerator gen = FFGenerator(); gen.hasItems(); gen.next())
582  {
583    alpha = gen.item();
584    gamma_image = reduce(gamma(alpha, x),M); // plug in alpha for var(1)
585    if(gamma_image.isZero()) // skip lc-bad points var(1)-alpha
586      continue;
587    TIMING_START (alg_recursion_p)
588    tryBrownGCD( f(alpha, x), g(alpha, x), M, g_image, fail, false ); // recursive call with one var less
589    TIMING_END_AND_PRINT (alg_recursion_p,
590                         "time for recursive calls in alg gcd mod p: ")
591    if(fail)
592      return;
593    g_image = reduce(g_image, M);
594    if(g_image.inCoeffDomain()) // early termination
595    {
596      tryInvert(g_image,M,result,fail);
597      if(fail)
598        return;
599      result = NN(c);
600      return;
601    }
602    for(int i=2; i<=mv; i++)
603      dg_im[i] = 0; // reset (this is necessary, because some entries may not be updated by call to leadDeg)
604    dg_im = leadDeg(g_image, dg_im); // dg_im cannot be NIL-pointer
605    if(isEqual(dg_im, L, 2, mv))
606    {
607      CanonicalForm inv;
608      tryInvert (firstLC (g_image), M, inv, fail);
609      if (fail)
610        return;
611      g_image *= inv;
612      g_image *= gamma_image; // multiply by multiple of image lc(gcd)
613      g_image= reduce (g_image, M);
614      TIMING_START (alg_newton_p)
615      gnew= tryNewtonInterp (alpha, g_image, m, gm, x, M, fail);
616      TIMING_END_AND_PRINT (alg_newton_p,
617                            "time for Newton interpolation in alg gcd mod p: ")
618      // gnew = gm mod m
619      // gnew = g_image mod var(1)-alpha
620      // mnew = m * (var(1)-alpha)
621      if(fail)
622        return;
623      m *= (x - alpha);
624      if((firstLC(gnew) == gamma) || (gnew == gm)) // gnew did not change
625      {
626        TIMING_START (alg_termination_p)
627        cf = tryvcontent(gnew, Variable(2), M, fail);
628        if(fail)
629          return;
630        divides = true;
631        g_image= gnew;
632        g_image.tryDiv (cf, M, fail);
633        if(fail)
634          return;
635        divides= tryFdivides (g_image,f, M, fail); // trial division (f)
636        if(fail)
637          return;
638        if(divides)
639        {
640          bool divides2= tryFdivides (g_image,g, M, fail); // trial division (g)
641          if(fail)
642            return;
643          if(divides2)
644          {
645            result = NN(reduce (c*g_image, M));
646            TIMING_END_AND_PRINT (alg_termination_p,
647                      "time for successful termination test in alg gcd mod p: ")
648            return;
649          }
650        }
651        TIMING_END_AND_PRINT (alg_termination_p,
652                    "time for unsuccessful termination test in alg gcd mod p: ")
653      }
654      gm = gnew;
655      continue;
656    }
657
658    if(isLess(L, dg_im, 2, mv)) // dg_im > L --> current point unlucky
659      continue;
660
661    // here: isLess(dg_im, L, 2, mv) --> all previous points were unlucky
662    m = CanonicalForm(1); // reset
663    gm = 0; // reset
664    for(int i=2; i<=mv; i++) // tighten bound
665      L[i] = dg_im[i];
666  }
667  // we are out of evaluation points
668  fail = true;
669}
670
671static CanonicalForm
672myicontent ( const CanonicalForm & f, const CanonicalForm & c )
673{
674#ifdef HAVE_NTL
675    if (f.isOne() || c.isOne())
676      return 1;
677    if ( f.inBaseDomain() && c.inBaseDomain())
678    {
679      if (c.isZero()) return abs(f);
680      return bgcd( f, c );
681    }
682    else if ( (f.inCoeffDomain() && c.inCoeffDomain()) ||
683              (f.inCoeffDomain() && c.inBaseDomain()) ||
684              (f.inBaseDomain() && c.inCoeffDomain()))
685    {
686      if (c.isZero()) return abs (f);
687#ifdef HAVE_FLINT
688      fmpz_poly_t FLINTf, FLINTc;
689      convertFacCF2Fmpz_poly_t (FLINTf, f);
690      convertFacCF2Fmpz_poly_t (FLINTc, c);
691      fmpz_poly_gcd (FLINTc, FLINTc, FLINTf);
692      CanonicalForm result;
693      if (f.inCoeffDomain())
694        result= convertFmpz_poly_t2FacCF (FLINTc, f.mvar());
695      else
696        result= convertFmpz_poly_t2FacCF (FLINTc, c.mvar());
697      fmpz_poly_clear (FLINTc);
698      fmpz_poly_clear (FLINTf);
699      return result;
700#else
701      ZZX NTLf= convertFacCF2NTLZZX (f);
702      ZZX NTLc= convertFacCF2NTLZZX (c);
703      NTLc= GCD (NTLc, NTLf);
704      if (f.inCoeffDomain())
705        return convertNTLZZX2CF(NTLc,f.mvar());
706      else
707        return convertNTLZZX2CF(NTLc,c.mvar());
708#endif
709    }
710    else
711    {
712        CanonicalForm g = c;
713        for ( CFIterator i = f; i.hasTerms() && ! g.isOne(); i++ )
714            g = myicontent( i.coeff(), g );
715        return g;
716    }
717#else
718    return 1;
719#endif
720}
721
722CanonicalForm
723myicontent ( const CanonicalForm & f )
724{
725#ifdef HAVE_NTL
726    return myicontent( f, 0 );
727#else
728    return 1;
729#endif
730}
731
732CanonicalForm QGCD( const CanonicalForm & F, const CanonicalForm & G )
733{ // f,g in Q(a)[x1,...,xn]
734  if(F.isZero())
735  {
736    if(G.isZero())
737      return G; // G is zero
738    if(G.inCoeffDomain())
739      return CanonicalForm(1);
740    CanonicalForm lcinv= 1/Lc (G);
741    return G*lcinv; // return monic G
742  }
743  if(G.isZero()) // F is non-zero
744  {
745    if(F.inCoeffDomain())
746      return CanonicalForm(1);
747    CanonicalForm lcinv= 1/Lc (F);
748    return F*lcinv; // return monic F
749  }
750  if(F.inCoeffDomain() || G.inCoeffDomain())
751    return CanonicalForm(1);
752  // here: both NOT inCoeffDomain
753  CanonicalForm f, g, tmp, M, q, D, Dp, cl, newq, mipo;
754  int p, i;
755  int *bound, *other; // degree vectors
756  bool fail;
757  bool off_rational=!isOn(SW_RATIONAL);
758  On( SW_RATIONAL ); // needed by bCommonDen
759  f = F * bCommonDen(F);
760  g = G * bCommonDen(G);
761  TIMING_START (alg_content)
762  CanonicalForm contf= myicontent (f);
763  CanonicalForm contg= myicontent (g);
764  f /= contf;
765  g /= contg;
766  CanonicalForm gcdcfcg= myicontent (contf, contg);
767  TIMING_END_AND_PRINT (alg_content, "time for content in alg gcd: ")
768  Variable a, b;
769  if(hasFirstAlgVar(f,a))
770  {
771    if(hasFirstAlgVar(g,b))
772    {
773      if(b.level() > a.level())
774        a = b;
775    }
776  }
777  else
778  {
779    if(!hasFirstAlgVar(g,a))// both not in extension
780    {
781      Off( SW_RATIONAL );
782      Off( SW_USE_QGCD );
783      tmp = gcdcfcg*gcd( f, g );
784      On( SW_USE_QGCD );
785      if (off_rational) Off(SW_RATIONAL);
786      return tmp;
787    }
788  }
789  // here: a is the biggest alg. var in f and g AND some of f,g is in extension
790  setReduce(a,false); // do not reduce expressions modulo mipo
791  tmp = getMipo(a);
792  M = tmp * bCommonDen(tmp);
793  // here: f, g in Z[a][x1,...,xn], M in Z[a] not necessarily monic
794  Off( SW_RATIONAL ); // needed by mod
795  // calculate upper bound for degree vector of gcd
796  int mv = f.level(); i = g.level();
797  if(i > mv)
798    mv = i;
799  // here: mv is level of the largest variable in f, g
800  bound = new int[mv+1]; // 'bound' could be indexed from 0 to mv, but we will only use from 1 to mv
801  other = new int[mv+1];
802  for(int i=1; i<=mv; i++) // initialize 'bound', 'other' with zeros
803    bound[i] = other[i] = 0;
804  bound = leadDeg(f,bound); // 'bound' is set the leading degree vector of f
805  other = leadDeg(g,other);
806  for(int i=1; i<=mv; i++) // entry at i=0 not visited
807    if(other[i] < bound[i])
808      bound[i] = other[i];
809  // now 'bound' is the smaller vector
810  cl = lc(M) * lc(f) * lc(g);
811  q = 1;
812  D = 0;
813  CanonicalForm test= 0;
814  bool equal= false;
815  for( i=cf_getNumBigPrimes()-1; i>-1; i-- )
816  {
817    p = cf_getBigPrime(i);
818    if( mod( cl, p ).isZero() ) // skip lc-bad primes
819      continue;
820    fail = false;
821    setCharacteristic(p);
822    mipo = mapinto(M);
823    mipo /= mipo.lc();
824    // here: mipo is monic
825    TIMING_START (alg_gcd_p)
826    tryBrownGCD( mapinto(f), mapinto(g), mipo, Dp, fail );
827    TIMING_END_AND_PRINT (alg_gcd_p, "time for alg gcd mod p: ")
828    if( fail ) // mipo splits in char p
829      continue;
830    if( Dp.inCoeffDomain() ) // early termination
831    {
832      tryInvert(Dp,mipo,tmp,fail); // check if zero divisor
833      if(fail)
834        continue;
835      setReduce(a,true);
836      if (off_rational) Off(SW_RATIONAL); else On(SW_RATIONAL);
837      setCharacteristic(0);
838      return gcdcfcg;
839    }
840    setCharacteristic(0);
841    // here: Dp NOT inCoeffDomain
842    for(int i=1; i<=mv; i++)
843      other[i] = 0; // reset (this is necessary, because some entries may not be updated by call to leadDeg)
844    other = leadDeg(Dp,other);
845
846    if(isEqual(bound, other, 1, mv)) // equal
847    {
848      chineseRemainder( D, q, mapinto(Dp), p, tmp, newq );
849      // tmp = Dp mod p
850      // tmp = D mod q
851      // newq = p*q
852      q = newq;
853      if( D != tmp )
854        D = tmp;
855      On( SW_RATIONAL );
856      TIMING_START (alg_reconstruction)
857      tmp = Farey( D, q ); // Farey
858      tmp *= bCommonDen (tmp);
859      TIMING_END_AND_PRINT (alg_reconstruction,
860                            "time for rational reconstruction in alg gcd: ")
861      setReduce(a,true); // reduce expressions modulo mipo
862      On( SW_RATIONAL ); // needed by fdivides
863      if (test != tmp)
864        test= tmp;
865      else
866        equal= true; // modular image did not add any new information
867      TIMING_START (alg_termination)
868#ifdef HAVE_NTL
869#ifdef HAVE_FLINT
870      if (equal && tmp.isUnivariate() && f.isUnivariate() && g.isUnivariate()
871          && f.level() == tmp.level() && tmp.level() == g.level())
872      {
873        CanonicalForm Q, R, sf, sg, stmp;
874        Variable x= Variable (1);
875        sf= swapvar (f, f.mvar(), x);
876        sg= swapvar (g, f.mvar(), x);
877        stmp= swapvar (tmp, f.mvar(), x);
878        newtonDivrem (sf, stmp, Q, R);
879        if (R.isZero())
880        {
881          newtonDivrem (sg, stmp, Q, R);
882          if (R.isZero())
883          {
884            Off (SW_RATIONAL);
885            setReduce (a,true);
886            if (off_rational) Off(SW_RATIONAL); else On(SW_RATIONAL);
887            TIMING_END_AND_PRINT (alg_termination,
888                                 "time for successful termination test in alg gcd: ")
889            return tmp*gcdcfcg;
890          }
891        }
892      }
893      else
894#endif
895#endif
896      if(equal && fdivides( tmp, f ) && fdivides( tmp, g )) // trial division
897      {
898        Off( SW_RATIONAL );
899        setReduce(a,true);
900        if (off_rational) Off(SW_RATIONAL); else On(SW_RATIONAL);
901        TIMING_END_AND_PRINT (alg_termination,
902                            "time for successful termination test in alg gcd: ")
903        return tmp*gcdcfcg;
904      }
905      TIMING_END_AND_PRINT (alg_termination,
906                          "time for unsuccessful termination test in alg gcd: ")
907      Off( SW_RATIONAL );
908      setReduce(a,false); // do not reduce expressions modulo mipo
909      continue;
910    }
911    if( isLess(bound, other, 1, mv) ) // current prime unlucky
912      continue;
913    // here: isLess(other, bound, 1, mv) ) ==> all previous primes unlucky
914    q = p;
915    D = mapinto(Dp); // shortcut CRA // shortcut CRA
916    for(int i=1; i<=mv; i++) // tighten bound
917      bound[i] = other[i];
918  }
919  // hopefully, we never reach this point
920  setReduce(a,true);
921  Off( SW_USE_QGCD );
922  D = gcdcfcg*gcd( f, g );
923  On( SW_USE_QGCD );
924  if (off_rational) Off(SW_RATIONAL); else On(SW_RATIONAL);
925  return D;
926}
927
928
929int * leadDeg(const CanonicalForm & f, int *degs)
930{ // leading degree vector w.r.t. lex. monomial order x(i+1) > x(i)
931  // if f is in a coeff domain, the zero pointer is returned
932  // 'a' should point to an array of sufficient size level(f)+1
933  if(f.inCoeffDomain())
934    return 0;
935  CanonicalForm tmp = f;
936  do
937  {
938    degs[tmp.level()] = tmp.degree();
939    tmp = LC(tmp);
940  }
941  while(!tmp.inCoeffDomain());
942  return degs;
943}
944
945
946bool isLess(int *a, int *b, int lower, int upper)
947{ // compares the degree vectors a,b on the specified part. Note: x(i+1) > x(i)
948  for(int i=upper; i>=lower; i--)
949    if(a[i] == b[i])
950      continue;
951    else
952      return a[i] < b[i];
953  return true;
954}
955
956
957bool isEqual(int *a, int *b, int lower, int upper)
958{ // compares the degree vectors a,b on the specified part. Note: x(i+1) > x(i)
959  for(int i=lower; i<=upper; i++)
960    if(a[i] != b[i])
961      return false;
962  return true;
963}
964
965
966CanonicalForm firstLC(const CanonicalForm & f)
967{ // returns the leading coefficient (LC) of level <= 1
968  CanonicalForm ret = f;
969  while(ret.level() > 1)
970    ret = LC(ret);
971  return ret;
972}
973
974void tryExtgcd( const CanonicalForm & F, const CanonicalForm & G, const CanonicalForm & M, CanonicalForm & result, CanonicalForm & s, CanonicalForm & t, bool & fail )
975{ // F, G are univariate polynomials (i.e. they have exactly one polynomial variable)
976  // F and G must have the same level AND level > 0
977  // we try to calculate gcd(F,G) = s*F + t*G
978  // if a zero divisor is encontered, 'fail' is set to one
979  // M is assumed to be monic
980  CanonicalForm P;
981  if(F.inCoeffDomain())
982  {
983    tryInvert( F, M, P, fail );
984    if(fail)
985      return;
986    result = 1;
987    s = P; t = 0;
988    return;
989  }
990  if(G.inCoeffDomain())
991  {
992    tryInvert( G, M, P, fail );
993    if(fail)
994      return;
995    result = 1;
996    s = 0; t = P;
997    return;
998  }
999  // here: both not inCoeffDomain
1000  CanonicalForm inv, rem, tmp, u, v, q, sum=0;
1001  if( F.degree() > G.degree() )
1002  {
1003    P = F; result = G;  s=v=0; t=u=1;
1004  }
1005  else
1006  {
1007    P = G; result = F; s=v=1; t=u=0;
1008  }
1009  Variable x = P.mvar();
1010  // here: degree(P) >= degree(result)
1011  while(true)
1012  {
1013    tryDivrem (P, result, q, rem, inv, M, fail);
1014    if(fail)
1015      return;
1016    if( rem.isZero() )
1017    {
1018      s*=inv;
1019      s= reduce (s, M);
1020      t*=inv;
1021      t= reduce (t, M);
1022      result *= inv; // monify result
1023      result= reduce (result, M);
1024      return;
1025    }
1026    sum += q;
1027    if(result.degree(x) >= rem.degree(x))
1028    {
1029      P=result;
1030      result=rem;
1031      tmp=u-sum*s;
1032      u=s;
1033      s=tmp;
1034      tmp=v-sum*t;
1035      v=t;
1036      t=tmp;
1037      sum = 0; // reset
1038    }
1039    else
1040      P = rem;
1041  }
1042}
1043
1044
1045static CanonicalForm trycontent ( const CanonicalForm & f, const Variable & x, const CanonicalForm & M, bool & fail )
1046{ // as 'content', but takes care of zero divisors
1047  ASSERT( x.level() > 0, "cannot calculate content with respect to algebraic variable" );
1048  Variable y = f.mvar();
1049  if ( y == x )
1050    return trycf_content( f, 0, M, fail );
1051  if ( y < x )
1052     return f;
1053  return swapvar( trycontent( swapvar( f, y, x ), y, M, fail ), y, x );
1054}
1055
1056
1057static CanonicalForm tryvcontent ( const CanonicalForm & f, const Variable & x, const CanonicalForm & M, bool & fail )
1058{ // as vcontent, but takes care of zero divisors
1059  ASSERT( x.level() > 0, "cannot calculate vcontent with respect to algebraic variable" );
1060  if ( f.mvar() <= x )
1061    return trycontent( f, x, M, fail );
1062  CFIterator i;
1063  CanonicalForm d = 0, e, ret;
1064  for ( i = f; i.hasTerms() && ! d.isOne() && ! fail; i++ )
1065  {
1066    e = tryvcontent( i.coeff(), x, M, fail );
1067    if(fail)
1068      break;
1069    tryBrownGCD( d, e, M, ret, fail );
1070    d = ret;
1071  }
1072  return d;
1073}
1074
1075
1076static CanonicalForm trycf_content ( const CanonicalForm & f, const CanonicalForm & g, const CanonicalForm & M, bool & fail )
1077{ // as cf_content, but takes care of zero divisors
1078  if ( f.inPolyDomain() || ( f.inExtension() && ! getReduce( f.mvar() ) ) )
1079  {
1080    CFIterator i = f;
1081    CanonicalForm tmp = g, result;
1082    while ( i.hasTerms() && ! tmp.isOne() && ! fail )
1083    {
1084      tryBrownGCD( i.coeff(), tmp, M, result, fail );
1085      tmp = result;
1086      i++;
1087    }
1088    return result;
1089  }
1090  return abs( f );
1091}
1092
1093void tryExtgcd( const CanonicalForm & F, const CanonicalForm & G, CanonicalForm & result, CanonicalForm & s, CanonicalForm & t, bool & fail )
1094{
1095  // F, G are univariate polynomials (i.e. they have exactly one polynomial variable)
1096  // F and G must have the same level AND level > 0
1097  // we try to calculate gcd(f,g) = s*F + t*G
1098  // if a zero divisor is encontered, 'fail' is set to one
1099  Variable a, b;
1100  if( !hasFirstAlgVar(F,a) && !hasFirstAlgVar(G,b) ) // note lazy evaluation
1101  {
1102    result = extgcd( F, G, s, t ); // no zero divisors possible
1103    return;
1104  }
1105  if( b.level() > a.level() )
1106    a = b;
1107  // here: a is the biggest alg. var in F and G
1108  CanonicalForm M = getMipo(a);
1109  CanonicalForm P;
1110  if( degree(F) > degree(G) )
1111  {
1112    P=F; result=G; s=0; t=1;
1113  }
1114  else
1115  {
1116    P=G; result=F; s=1; t=0;
1117  }
1118  CanonicalForm inv, rem, q, u, v;
1119  // here: degree(P) >= degree(result)
1120  while(true)
1121  {
1122    tryInvert( Lc(result), M, inv, fail );
1123    if(fail)
1124      return;
1125    // here: Lc(result) is invertible modulo M
1126    q = Lc(P)*inv * power( P.mvar(), degree(P)-degree(result) );
1127    rem = P - q*result;
1128    // here: s*F + t*G = result
1129    if( rem.isZero() )
1130    {
1131      s*=inv;
1132      t*=inv;
1133      result *= inv; // monify result
1134      return;
1135    }
1136    P=result;
1137    result=rem;
1138    rem=u-q*s;
1139    u=s;
1140    s=rem;
1141    rem=v-q*t;
1142    v=t;
1143    t=rem;
1144  }
1145}
1146
1147void tryCRA( const CanonicalForm & x1, const CanonicalForm & q1, const CanonicalForm & x2, const CanonicalForm & q2, CanonicalForm & xnew, CanonicalForm & qnew, bool & fail )
1148{ // polys of level <= 1 are considered coefficients. q1,q2 are assumed to be coprime
1149  // xnew = x1 mod q1 (coefficientwise in the above sense)
1150  // xnew = x2 mod q2
1151  // qnew = q1*q2
1152  CanonicalForm tmp;
1153  if(x1.level() <= 1 && x2.level() <= 1) // base case
1154  {
1155    tryExtgcd(q1,q2,tmp,xnew,qnew,fail);
1156    if(fail)
1157      return;
1158    xnew = x1 + (x2-x1) * xnew * q1;
1159    qnew = q1*q2;
1160    xnew = mod(xnew,qnew);
1161    return;
1162  }
1163  CanonicalForm tmp2;
1164  xnew = 0;
1165  qnew = q1 * q2;
1166  // here: x1.level() > 1 || x2.level() > 1
1167  if(x1.level() > x2.level())
1168  {
1169    for(CFIterator i=x1; i.hasTerms(); i++)
1170    {
1171      if(i.exp() == 0) // const. term
1172      {
1173        tryCRA(i.coeff(),q1,x2,q2,tmp,tmp2,fail);
1174        if(fail)
1175          return;
1176        xnew += tmp;
1177      }
1178      else
1179      {
1180        tryCRA(i.coeff(),q1,0,q2,tmp,tmp2,fail);
1181        if(fail)
1182          return;
1183        xnew += tmp * power(x1.mvar(),i.exp());
1184      }
1185    }
1186    return;
1187  }
1188  // here: x1.level() <= x2.level() && ( x1.level() > 1 || x2.level() > 1 )
1189  if(x2.level() > x1.level())
1190  {
1191    for(CFIterator j=x2; j.hasTerms(); j++)
1192    {
1193      if(j.exp() == 0) // const. term
1194      {
1195        tryCRA(x1,q1,j.coeff(),q2,tmp,tmp2,fail);
1196        if(fail)
1197          return;
1198        xnew += tmp;
1199      }
1200      else
1201      {
1202        tryCRA(0,q1,j.coeff(),q2,tmp,tmp2,fail);
1203        if(fail)
1204          return;
1205        xnew += tmp * power(x2.mvar(),j.exp());
1206      }
1207    }
1208    return;
1209  }
1210  // here: x1.level() == x2.level() && x1.level() > 1 && x2.level() > 1
1211  CFIterator i = x1;
1212  CFIterator j = x2;
1213  while(i.hasTerms() || j.hasTerms())
1214  {
1215    if(i.hasTerms())
1216    {
1217      if(j.hasTerms())
1218      {
1219        if(i.exp() == j.exp())
1220        {
1221          tryCRA(i.coeff(),q1,j.coeff(),q2,tmp,tmp2,fail);
1222          if(fail)
1223            return;
1224          xnew += tmp * power(x1.mvar(),i.exp());
1225          i++; j++;
1226        }
1227        else
1228        {
1229          if(i.exp() < j.exp())
1230          {
1231            tryCRA(i.coeff(),q1,0,q2,tmp,tmp2,fail);
1232            if(fail)
1233              return;
1234            xnew += tmp * power(x1.mvar(),i.exp());
1235            i++;
1236          }
1237          else // i.exp() > j.exp()
1238          {
1239            tryCRA(0,q1,j.coeff(),q2,tmp,tmp2,fail);
1240            if(fail)
1241              return;
1242            xnew += tmp * power(x1.mvar(),j.exp());
1243            j++;
1244          }
1245        }
1246      }
1247      else // j is out of terms
1248      {
1249        tryCRA(i.coeff(),q1,0,q2,tmp,tmp2,fail);
1250        if(fail)
1251          return;
1252        xnew += tmp * power(x1.mvar(),i.exp());
1253        i++;
1254      }
1255    }
1256    else // i is out of terms
1257    {
1258      tryCRA(0,q1,j.coeff(),q2,tmp,tmp2,fail);
1259      if(fail)
1260        return;
1261      xnew += tmp * power(x1.mvar(),j.exp());
1262      j++;
1263    }
1264  }
1265}
1266
Note: See TracBrowser for help on using the repository browser.