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