source: git/factory/algext.cc @ aba90e8

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