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