cfModGcd.cc
Go to the documentation of this file.
1 // -*- c++ -*-
2 //*****************************************************************************
3 /** @file cfModGcd.cc
4  *
5  * This file implements the GCD of two polynomials over \f$ F_{p} \f$ ,
6  * \f$ F_{p}(\alpha ) \f$, GF or Z based on Alg. 7.1. and 7.2. as described in
7  * "Algorithms for Computer Algebra" by Geddes, Czapor, Labahn via modular
8  * computations. And sparse modular variants as described in Zippel
9  * "Effective Polynomial Computation", deKleine, Monagan, Wittkopf
10  * "Algorithms for the non-monic case of the sparse modular GCD algorithm" and
11  * Javadi "A new solution to the normalization problem"
12  *
13  * @author Martin Lee
14  * @date 22.10.2009
15  *
16  * @par Copyright:
17  * (c) by The SINGULAR Team, see LICENSE file
18  *
19 **/
20 //*****************************************************************************
21 
22 
23 #include "config.h"
24 
25 
26 #include "cf_assert.h"
27 #include "debug.h"
28 #include "timing.h"
29 
30 #include "canonicalform.h"
31 #include "cfGcdUtil.h"
32 #include "cf_map.h"
33 #include "cf_util.h"
34 #include "cf_irred.h"
36 #include "cf_random.h"
37 #include "cf_reval.h"
38 #include "facHensel.h"
39 #include "cf_iter.h"
40 #include "cfNewtonPolygon.h"
41 #include "cf_algorithm.h"
42 #include "cf_primes.h"
43 
44 // inline helper functions:
45 #include "cf_map_ext.h"
46 
47 #ifdef HAVE_NTL
48 #include "NTLconvert.h"
49 #endif
50 
51 #ifdef HAVE_FLINT
52 #include "FLINTconvert.h"
53 #endif
54 
55 #include "cfModGcd.h"
56 
57 TIMING_DEFINE_PRINT(gcd_recursion)
58 TIMING_DEFINE_PRINT(newton_interpolation)
59 TIMING_DEFINE_PRINT(termination_test)
60 TIMING_DEFINE_PRINT(ez_p_compress)
61 TIMING_DEFINE_PRINT(ez_p_hensel_lift)
62 TIMING_DEFINE_PRINT(ez_p_eval)
63 TIMING_DEFINE_PRINT(ez_p_content)
64 
65 bool
67  const CanonicalForm& coF, const CanonicalForm& coG,
68  const CanonicalForm& cand)
69 {
70  CanonicalForm LCCand= abs (LC (cand));
71  if (LCCand*abs (LC (coF)) == abs (LC (F)))
72  {
73  if (LCCand*abs (LC (coG)) == abs (LC (G)))
74  {
75  if (abs (cand)*abs (coF) == abs (F))
76  {
77  if (abs (cand)*abs (coG) == abs (G))
78  return true;
79  }
80  return false;
81  }
82  return false;
83  }
84  return false;
85 }
86 
87 #ifdef HAVE_NTL
88 
89 static const double log2exp= 1.442695041;
90 
91 /// compressing two polynomials F and G, M is used for compressing,
92 /// N to reverse the compression
93 int myCompress (const CanonicalForm& F, const CanonicalForm& G, CFMap & M,
94  CFMap & N, bool topLevel)
95 {
96  int n= tmax (F.level(), G.level());
97  int * degsf= NEW_ARRAY(int,n + 1);
98  int * degsg= NEW_ARRAY(int,n + 1);
99 
100  for (int i = n; i >= 0; i--)
101  degsf[i]= degsg[i]= 0;
102 
103  degsf= degrees (F, degsf);
104  degsg= degrees (G, degsg);
105 
106  int both_non_zero= 0;
107  int f_zero= 0;
108  int g_zero= 0;
109  int both_zero= 0;
110 
111  if (topLevel)
112  {
113  for (int i= 1; i <= n; i++)
114  {
115  if (degsf[i] != 0 && degsg[i] != 0)
116  {
117  both_non_zero++;
118  continue;
119  }
120  if (degsf[i] == 0 && degsg[i] != 0 && i <= G.level())
121  {
122  f_zero++;
123  continue;
124  }
125  if (degsg[i] == 0 && degsf[i] && i <= F.level())
126  {
127  g_zero++;
128  continue;
129  }
130  }
131 
132  if (both_non_zero == 0)
133  {
134  DELETE_ARRAY(degsf);
135  DELETE_ARRAY(degsg);
136  return 0;
137  }
138 
139  // map Variables which do not occur in both polynomials to higher levels
140  int k= 1;
141  int l= 1;
142  for (int i= 1; i <= n; i++)
143  {
144  if (degsf[i] != 0 && degsg[i] == 0 && i <= F.level())
145  {
146  if (k + both_non_zero != i)
147  {
148  M.newpair (Variable (i), Variable (k + both_non_zero));
149  N.newpair (Variable (k + both_non_zero), Variable (i));
150  }
151  k++;
152  }
153  if (degsf[i] == 0 && degsg[i] != 0 && i <= G.level())
154  {
155  if (l + g_zero + both_non_zero != i)
156  {
157  M.newpair (Variable (i), Variable (l + g_zero + both_non_zero));
158  N.newpair (Variable (l + g_zero + both_non_zero), Variable (i));
159  }
160  l++;
161  }
162  }
163 
164  // sort Variables x_{i} in increasing order of
165  // min(deg_{x_{i}}(f),deg_{x_{i}}(g))
166  int m= tmax (F.level(), G.level());
167  int min_max_deg;
168  k= both_non_zero;
169  l= 0;
170  int i= 1;
171  while (k > 0)
172  {
173  if (degsf [i] != 0 && degsg [i] != 0)
174  min_max_deg= tmax (degsf[i], degsg[i]);
175  else
176  min_max_deg= 0;
177  while (min_max_deg == 0)
178  {
179  i++;
180  if (degsf [i] != 0 && degsg [i] != 0)
181  min_max_deg= tmax (degsf[i], degsg[i]);
182  else
183  min_max_deg= 0;
184  }
185  for (int j= i + 1; j <= m; j++)
186  {
187  if (degsf[j] != 0 && degsg [j] != 0 &&
188  tmax (degsf[j],degsg[j]) <= min_max_deg)
189  {
190  min_max_deg= tmax (degsf[j], degsg[j]);
191  l= j;
192  }
193  }
194  if (l != 0)
195  {
196  if (l != k)
197  {
198  M.newpair (Variable (l), Variable(k));
199  N.newpair (Variable (k), Variable(l));
200  degsf[l]= 0;
201  degsg[l]= 0;
202  l= 0;
203  }
204  else
205  {
206  degsf[l]= 0;
207  degsg[l]= 0;
208  l= 0;
209  }
210  }
211  else if (l == 0)
212  {
213  if (i != k)
214  {
215  M.newpair (Variable (i), Variable (k));
216  N.newpair (Variable (k), Variable (i));
217  degsf[i]= 0;
218  degsg[i]= 0;
219  }
220  else
221  {
222  degsf[i]= 0;
223  degsg[i]= 0;
224  }
225  i++;
226  }
227  k--;
228  }
229  }
230  else
231  {
232  //arrange Variables such that no gaps occur
233  for (int i= 1; i <= n; i++)
234  {
235  if (degsf[i] == 0 && degsg[i] == 0)
236  {
237  both_zero++;
238  continue;
239  }
240  else
241  {
242  if (both_zero != 0)
243  {
244  M.newpair (Variable (i), Variable (i - both_zero));
245  N.newpair (Variable (i - both_zero), Variable (i));
246  }
247  }
248  }
249  }
250 
251  DELETE_ARRAY(degsf);
252  DELETE_ARRAY(degsg);
253 
254  return 1;
255 }
256 
257 static inline CanonicalForm
258 uni_content (const CanonicalForm & F);
259 
262 {
263  if (F.inCoeffDomain())
264  return F.genOne();
265  if (F.level() == x.level() && F.isUnivariate())
266  return F;
267  if (F.level() != x.level() && F.isUnivariate())
268  return F.genOne();
269 
270  if (x.level() != 1)
271  {
272  CanonicalForm f= swapvar (F, x, Variable (1));
274  return swapvar (result, x, Variable (1));
275  }
276  else
277  return uni_content (F);
278 }
279 
280 /// compute the content of F, where F is considered as an element of
281 /// \f$ R[x_{1}][x_{2},\ldots ,x_{n}] \f$
282 static inline CanonicalForm
284 {
285  if (F.inBaseDomain())
286  return F.genOne();
287  if (F.level() == 1 && F.isUnivariate())
288  return F;
289  if (F.level() != 1 && F.isUnivariate())
290  return F.genOne();
291  if (degree (F,1) == 0) return F.genOne();
292 
293  int l= F.level();
294  if (l == 2)
295  return content(F);
296  else
297  {
298  CanonicalForm pol, c = 0;
299  CFIterator i = F;
300  for (; i.hasTerms(); i++)
301  {
302  pol= i.coeff();
303  pol= uni_content (pol);
304  c= gcd (c, pol);
305  if (c.isOne())
306  return c;
307  }
308  return c;
309  }
310 }
311 
314  CanonicalForm& contentF, CanonicalForm& contentG,
315  CanonicalForm& ppF, CanonicalForm& ppG, const int d)
316 {
317  CanonicalForm uniContentF, uniContentG, gcdcFcG;
318  contentF= 1;
319  contentG= 1;
320  ppF= F;
321  ppG= G;
323  for (int i= 1; i <= d; i++)
324  {
325  uniContentF= uni_content (F, Variable (i));
326  uniContentG= uni_content (G, Variable (i));
327  gcdcFcG= gcd (uniContentF, uniContentG);
328  contentF *= uniContentF;
329  contentG *= uniContentG;
330  ppF /= uniContentF;
331  ppG /= uniContentG;
332  result *= gcdcFcG;
333  }
334  return result;
335 }
336 
337 /// compute the leading coefficient of F, where F is considered as an element
338 /// of \f$ R[x_{1}][x_{2},\ldots ,x_{n}] \f$, order on
339 /// \f$ Mon (x_{2},\ldots ,x_{n}) \f$ is dp.
340 static inline
342 {
343  if (F.level() > 1)
344  {
345  Variable x= Variable (2);
346  int deg= totaldegree (F, x, F.mvar());
347  for (CFIterator i= F; i.hasTerms(); i++)
348  {
349  if (i.exp() + totaldegree (i.coeff(), x, i.coeff().mvar()) == deg)
350  return uni_lcoeff (i.coeff());
351  }
352  }
353  return F;
354 }
355 
356 /// Newton interpolation - Incremental algorithm.
357 /// Given a list of values alpha_i and a list of polynomials u_i, 1 <= i <= n,
358 /// computes the interpolation polynomial assuming that
359 /// the polynomials in u are the results of evaluating the variabe x
360 /// of the unknown polynomial at the alpha values.
361 /// This incremental version receives only the values of alpha_n and u_n and
362 /// the previous interpolation polynomial for points 1 <= i <= n-1, and computes
363 /// the polynomial interpolating in all the points.
364 /// newtonPoly must be equal to (x - alpha_1) * ... * (x - alpha_{n-1})
365 static inline CanonicalForm
367  const CanonicalForm & newtonPoly, const CanonicalForm & oldInterPoly,
368  const Variable & x)
369 {
370  CanonicalForm interPoly;
371 
372  interPoly= oldInterPoly + ((u - oldInterPoly(alpha, x))/newtonPoly(alpha, x))
373  *newtonPoly;
374  return interPoly;
375 }
376 
377 /// compute a random element a of \f$ F_{p}(\alpha ) \f$ ,
378 /// s.t. F(a) \f$ \neq 0 \f$ , F is a univariate polynomial, returns
379 /// fail if there are no field elements left which have not been used before
380 static inline CanonicalForm
381 randomElement (const CanonicalForm & F, const Variable & alpha, CFList & list,
382  bool & fail)
383 {
384  fail= false;
385  Variable x= F.mvar();
386  AlgExtRandomF genAlgExt (alpha);
387  FFRandom genFF;
388  CanonicalForm random, mipo;
389  mipo= getMipo (alpha);
390  int p= getCharacteristic ();
391  int d= degree (mipo);
392  double bound= pow ((double) p, (double) d);
393  do
394  {
395  if (list.length() == bound)
396  {
397  fail= true;
398  break;
399  }
400  if (list.length() < p)
401  {
402  random= genFF.generate();
403  while (find (list, random))
404  random= genFF.generate();
405  }
406  else
407  {
408  random= genAlgExt.generate();
409  while (find (list, random))
410  random= genAlgExt.generate();
411  }
412  if (F (random, x) == 0)
413  {
414  list.append (random);
415  continue;
416  }
417  } while (find (list, random));
418  return random;
419 }
420 
421 static inline
423 {
425  {
427  zz_p::init (getCharacteristic());
428  }
429  zz_pX NTLIrredpoly;
430  int i, m;
431  // extension of F_p needed
432  if (alpha.level() == 1)
433  {
434  i= 1;
435  m= 2;
436  } //extension of F_p(alpha)
437  if (alpha.level() != 1)
438  {
439  i= 4;
440  m= degree (getMipo (alpha));
441  }
442  BuildIrred (NTLIrredpoly, i*m);
443  CanonicalForm newMipo= convertNTLzzpX2CF (NTLIrredpoly, Variable (1));
444  return rootOf (newMipo);
445 }
446 
448 modGCDFq (const CanonicalForm& F, const CanonicalForm& G,
449  CanonicalForm& coF, CanonicalForm& coG,
450  Variable & alpha, CFList& l, bool& topLevel);
451 
454  Variable & alpha, CFList& l, bool& topLevel)
455 {
456  CanonicalForm dummy1, dummy2;
457  CanonicalForm result= modGCDFq (F, G, dummy1, dummy2, alpha, l,
458  topLevel);
459  return result;
460 }
461 
462 /// GCD of F and G over \f$ F_{p}(\alpha ) \f$ ,
463 /// l and topLevel are only used internally, output is monic
464 /// based on Alg. 7.2. as described in "Algorithms for
465 /// Computer Algebra" by Geddes, Czapor, Labahn
468  CanonicalForm& coF, CanonicalForm& coG,
469  Variable & alpha, CFList& l, bool& topLevel)
470 {
471  CanonicalForm A= F;
472  CanonicalForm B= G;
473  if (F.isZero() && degree(G) > 0)
474  {
475  coF= 0;
476  coG= Lc (G);
477  return G/Lc(G);
478  }
479  else if (G.isZero() && degree (F) > 0)
480  {
481  coF= Lc (F);
482  coG= 0;
483  return F/Lc(F);
484  }
485  else if (F.isZero() && G.isZero())
486  {
487  coF= coG= 0;
488  return F.genOne();
489  }
490  if (F.inBaseDomain() || G.inBaseDomain())
491  {
492  coF= F;
493  coG= G;
494  return F.genOne();
495  }
496  if (F.isUnivariate() && fdivides(F, G, coG))
497  {
498  coF= Lc (F);
499  return F/Lc(F);
500  }
501  if (G.isUnivariate() && fdivides(G, F, coF))
502  {
503  coG= Lc (G);
504  return G/Lc(G);
505  }
506  if (F == G)
507  {
508  coF= coG= Lc (F);
509  return F/Lc(F);
510  }
511 
512  CFMap M,N;
513  int best_level= myCompress (A, B, M, N, topLevel);
514 
515  if (best_level == 0)
516  {
517  coF= F;
518  coG= G;
519  return B.genOne();
520  }
521 
522  A= M(A);
523  B= M(B);
524 
525  Variable x= Variable(1);
526 
527  //univariate case
528  if (A.isUnivariate() && B.isUnivariate())
529  {
530  CanonicalForm result= gcd (A, B);
531  coF= N (A/result);
532  coG= N (B/result);
533  return N (result);
534  }
535 
536  CanonicalForm cA, cB; // content of A and B
537  CanonicalForm ppA, ppB; // primitive part of A and B
538  CanonicalForm gcdcAcB;
539 
540  cA = uni_content (A);
541  cB = uni_content (B);
542  gcdcAcB= gcd (cA, cB);
543  ppA= A/cA;
544  ppB= B/cB;
545 
546  CanonicalForm lcA, lcB; // leading coefficients of A and B
547  CanonicalForm gcdlcAlcB;
548 
549  lcA= uni_lcoeff (ppA);
550  lcB= uni_lcoeff (ppB);
551 
552  gcdlcAlcB= gcd (lcA, lcB);
553 
554  int d= totaldegree (ppA, Variable(2), Variable (ppA.level()));
555 
556  if (d == 0)
557  {
558  coF= N (ppA*(cA/gcdcAcB));
559  coG= N (ppB*(cB/gcdcAcB));
560  return N(gcdcAcB);
561  }
562 
563  int d0= totaldegree (ppB, Variable(2), Variable (ppB.level()));
564  if (d0 < d)
565  d= d0;
566  if (d == 0)
567  {
568  coF= N (ppA*(cA/gcdcAcB));
569  coG= N (ppB*(cB/gcdcAcB));
570  return N(gcdcAcB);
571  }
572 
573  CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH;
574  CanonicalForm newtonPoly, coF_random_element, coG_random_element, coF_m,
575  coG_m, ppCoF, ppCoG;
576 
577  newtonPoly= 1;
578  m= gcdlcAlcB;
579  G_m= 0;
580  coF= 0;
581  coG= 0;
582  H= 0;
583  bool fail= false;
584  topLevel= false;
585  bool inextension= false;
586  Variable V_buf= alpha, V_buf4= alpha;
587  CanonicalForm prim_elem, im_prim_elem;
588  CanonicalForm prim_elem_alpha, im_prim_elem_alpha;
589  CFList source, dest;
590  int bound1= degree (ppA, 1);
591  int bound2= degree (ppB, 1);
592  do
593  {
594  random_element= randomElement (m*lcA*lcB, V_buf, l, fail);
595  if (fail)
596  {
597  source= CFList();
598  dest= CFList();
599 
600  Variable V_buf3= V_buf;
601  V_buf= chooseExtension (V_buf);
602  bool prim_fail= false;
603  Variable V_buf2;
604  prim_elem= primitiveElement (V_buf4, V_buf2, prim_fail);
605  if (V_buf4 == alpha)
606  prim_elem_alpha= prim_elem;
607 
608  if (V_buf3 != V_buf4)
609  {
610  m= mapDown (m, prim_elem, im_prim_elem, V_buf4, source, dest);
611  G_m= mapDown (G_m, prim_elem, im_prim_elem, V_buf4, source, dest);
612  coF_m= mapDown (coF_m, prim_elem, im_prim_elem, V_buf4, source, dest);
613  coG_m= mapDown (coG_m, prim_elem, im_prim_elem, V_buf4, source, dest);
614  newtonPoly= mapDown (newtonPoly, prim_elem, im_prim_elem, V_buf4,
615  source, dest);
616  ppA= mapDown (ppA, prim_elem, im_prim_elem, V_buf4, source, dest);
617  ppB= mapDown (ppB, prim_elem, im_prim_elem, V_buf4, source, dest);
618  gcdlcAlcB= mapDown (gcdlcAlcB, prim_elem, im_prim_elem, V_buf4,
619  source, dest);
620  lcA= mapDown (lcA, prim_elem, im_prim_elem, V_buf4, source, dest);
621  lcB= mapDown (lcB, prim_elem, im_prim_elem, V_buf4, source, dest);
622  for (CFListIterator i= l; i.hasItem(); i++)
623  i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, V_buf4,
624  source, dest);
625  }
626 
627  ASSERT (!prim_fail, "failure in integer factorizer");
628  if (prim_fail)
629  ; //ERROR
630  else
631  im_prim_elem= mapPrimElem (prim_elem, V_buf4, V_buf);
632 
633  if (V_buf4 == alpha)
634  im_prim_elem_alpha= im_prim_elem;
635  else
636  im_prim_elem_alpha= mapUp (im_prim_elem_alpha, V_buf4, V_buf, prim_elem,
637  im_prim_elem, source, dest);
638  DEBOUTLN (cerr, "getMipo (V_buf4)= " << getMipo (V_buf4));
639  DEBOUTLN (cerr, "getMipo (V_buf2)= " << getMipo (V_buf2));
640  inextension= true;
641  for (CFListIterator i= l; i.hasItem(); i++)
642  i.getItem()= mapUp (i.getItem(), V_buf4, V_buf, prim_elem,
643  im_prim_elem, source, dest);
644  m= mapUp (m, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
645  G_m= mapUp (G_m, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
646  coF_m= mapUp (coF_m, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
647  coG_m= mapUp (coG_m, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
648  newtonPoly= mapUp (newtonPoly, V_buf4, V_buf, prim_elem, im_prim_elem,
649  source, dest);
650  ppA= mapUp (ppA, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
651  ppB= mapUp (ppB, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
652  gcdlcAlcB= mapUp (gcdlcAlcB, V_buf4, V_buf, prim_elem, im_prim_elem,
653  source, dest);
654  lcA= mapUp (lcA, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
655  lcB= mapUp (lcB, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
656 
657  fail= false;
658  random_element= randomElement (m*lcA*lcB, V_buf, l, fail );
659  DEBOUTLN (cerr, "fail= " << fail);
660  CFList list;
661  TIMING_START (gcd_recursion);
662  G_random_element=
663  modGCDFq (ppA (random_element, x), ppB (random_element, x),
664  coF_random_element, coG_random_element, V_buf,
665  list, topLevel);
666  TIMING_END_AND_PRINT (gcd_recursion,
667  "time for recursive call: ");
668  DEBOUTLN (cerr, "G_random_element= " << G_random_element);
669  V_buf4= V_buf;
670  }
671  else
672  {
673  CFList list;
674  TIMING_START (gcd_recursion);
675  G_random_element=
676  modGCDFq (ppA(random_element, x), ppB(random_element, x),
677  coF_random_element, coG_random_element, V_buf,
678  list, topLevel);
679  TIMING_END_AND_PRINT (gcd_recursion,
680  "time for recursive call: ");
681  DEBOUTLN (cerr, "G_random_element= " << G_random_element);
682  }
683 
684  if (!G_random_element.inCoeffDomain())
685  d0= totaldegree (G_random_element, Variable(2),
686  Variable (G_random_element.level()));
687  else
688  d0= 0;
689 
690  if (d0 == 0)
691  {
692  if (inextension)
693  {
694  CFList u, v;
695  ppA= mapDown (ppA, prim_elem_alpha, im_prim_elem_alpha, alpha, u, v);
696  ppB= mapDown (ppB, prim_elem_alpha, im_prim_elem_alpha, alpha, u, v);
697  prune1 (alpha);
698  }
699  coF= N (ppA*(cA/gcdcAcB));
700  coG= N (ppB*(cB/gcdcAcB));
701  return N(gcdcAcB);
702  }
703  if (d0 > d)
704  {
705  if (!find (l, random_element))
706  l.append (random_element);
707  continue;
708  }
709 
710  G_random_element=
711  (gcdlcAlcB(random_element, x)/uni_lcoeff (G_random_element))
712  * G_random_element;
713  coF_random_element= (lcA(random_element,x)/uni_lcoeff(coF_random_element))
714  *coF_random_element;
715  coG_random_element= (lcB(random_element,x)/uni_lcoeff(coG_random_element))
716  *coG_random_element;
717 
718  if (!G_random_element.inCoeffDomain())
719  d0= totaldegree (G_random_element, Variable(2),
720  Variable (G_random_element.level()));
721  else
722  d0= 0;
723 
724  if (d0 < d)
725  {
726  m= gcdlcAlcB;
727  newtonPoly= 1;
728  G_m= 0;
729  d= d0;
730  coF_m= 0;
731  coG_m= 0;
732  }
733 
734  TIMING_START (newton_interpolation);
735  H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
736  coF= newtonInterp (random_element, coF_random_element, newtonPoly, coF_m,x);
737  coG= newtonInterp (random_element, coG_random_element, newtonPoly, coG_m,x);
738  TIMING_END_AND_PRINT (newton_interpolation,
739  "time for newton interpolation: ");
740 
741  //termination test
742  if ((uni_lcoeff (H) == gcdlcAlcB) || (G_m == H))
743  {
744  TIMING_START (termination_test);
745  if (gcdlcAlcB.isOne())
746  cH= 1;
747  else
748  cH= uni_content (H);
749  ppH= H/cH;
750  ppH /= Lc (ppH);
751  CanonicalForm lcppH= gcdlcAlcB/cH;
752  CanonicalForm ccoF= lcppH/Lc (lcppH);
753  CanonicalForm ccoG= lcppH/Lc (lcppH);
754  ppCoF= coF/ccoF;
755  ppCoG= coG/ccoG;
756  if (inextension)
757  {
758  if (((degree (ppCoF,1)+degree (ppH,1) == bound1) &&
759  (degree (ppCoG,1)+degree (ppH,1) == bound2) &&
760  terminationTest (ppA, ppB, ppCoF, ppCoG, ppH)) ||
761  (fdivides (ppH, ppA, ppCoF) && fdivides (ppH, ppB, ppCoG)))
762  {
763  CFList u, v;
764  DEBOUTLN (cerr, "ppH before mapDown= " << ppH);
765  ppH= mapDown (ppH, prim_elem_alpha, im_prim_elem_alpha, alpha, u, v);
766  ppCoF= mapDown (ppCoF, prim_elem_alpha, im_prim_elem_alpha, alpha, u, v);
767  ppCoG= mapDown (ppCoG, prim_elem_alpha, im_prim_elem_alpha, alpha, u, v);
768  DEBOUTLN (cerr, "ppH after mapDown= " << ppH);
769  coF= N ((cA/gcdcAcB)*ppCoF);
770  coG= N ((cB/gcdcAcB)*ppCoG);
771  TIMING_END_AND_PRINT (termination_test,
772  "time for successful termination test Fq: ");
773  prune1 (alpha);
774  return N(gcdcAcB*ppH);
775  }
776  }
777  else if (((degree (ppCoF,1)+degree (ppH,1) == bound1) &&
778  (degree (ppCoG,1)+degree (ppH,1) == bound2) &&
779  terminationTest (ppA, ppB, ppCoF, ppCoG, ppH)) ||
780  (fdivides (ppH, ppA, ppCoF) && fdivides (ppH, ppB, ppCoG)))
781  {
782  coF= N ((cA/gcdcAcB)*ppCoF);
783  coG= N ((cB/gcdcAcB)*ppCoG);
784  TIMING_END_AND_PRINT (termination_test,
785  "time for successful termination test Fq: ");
786  return N(gcdcAcB*ppH);
787  }
788  TIMING_END_AND_PRINT (termination_test,
789  "time for unsuccessful termination test Fq: ");
790  }
791 
792  G_m= H;
793  coF_m= coF;
794  coG_m= coG;
795  newtonPoly= newtonPoly*(x - random_element);
796  m= m*(x - random_element);
797  if (!find (l, random_element))
798  l.append (random_element);
799  } while (1);
800 }
801 
802 /// compute a random element a of GF, s.t. F(a) \f$ \neq 0 \f$ , F is a
803 /// univariate polynomial, returns fail if there are no field elements left
804 /// which have not been used before
805 static inline
807 GFRandomElement (const CanonicalForm& F, CFList& list, bool& fail)
808 {
809  fail= false;
810  Variable x= F.mvar();
811  GFRandom genGF;
812  CanonicalForm random;
813  int p= getCharacteristic();
814  int d= getGFDegree();
815  int bound= ipower (p, d);
816  do
817  {
818  if (list.length() == bound)
819  {
820  fail= true;
821  break;
822  }
823  if (list.length() < 1)
824  random= 0;
825  else
826  {
827  random= genGF.generate();
828  while (find (list, random))
829  random= genGF.generate();
830  }
831  if (F (random, x) == 0)
832  {
833  list.append (random);
834  continue;
835  }
836  } while (find (list, random));
837  return random;
838 }
839 
841 modGCDGF (const CanonicalForm& F, const CanonicalForm& G,
842  CanonicalForm& coF, CanonicalForm& coG,
843  CFList& l, bool& topLevel);
844 
846 modGCDGF (const CanonicalForm& F, const CanonicalForm& G, CFList& l,
847  bool& topLevel)
848 {
849  CanonicalForm dummy1, dummy2;
850  CanonicalForm result= modGCDGF (F, G, dummy1, dummy2, l, topLevel);
851  return result;
852 }
853 
854 /// GCD of F and G over GF, based on Alg. 7.2. as described in "Algorithms for
855 /// Computer Algebra" by Geddes, Czapor, Labahn
856 /// Usually this algorithm will be faster than modGCDFq since GF has
857 /// faster field arithmetics, however it might fail if the input is large since
858 /// the size of the base field is bounded by 2^16, output is monic
861  CanonicalForm& coF, CanonicalForm& coG,
862  CFList& l, bool& topLevel)
863 {
864  CanonicalForm A= F;
865  CanonicalForm B= G;
866  if (F.isZero() && degree(G) > 0)
867  {
868  coF= 0;
869  coG= Lc (G);
870  return G/Lc(G);
871  }
872  else if (G.isZero() && degree (F) > 0)
873  {
874  coF= Lc (F);
875  coG= 0;
876  return F/Lc(F);
877  }
878  else if (F.isZero() && G.isZero())
879  {
880  coF= coG= 0;
881  return F.genOne();
882  }
883  if (F.inBaseDomain() || G.inBaseDomain())
884  {
885  coF= F;
886  coG= G;
887  return F.genOne();
888  }
889  if (F.isUnivariate() && fdivides(F, G, coG))
890  {
891  coF= Lc (F);
892  return F/Lc(F);
893  }
894  if (G.isUnivariate() && fdivides(G, F, coF))
895  {
896  coG= Lc (G);
897  return G/Lc(G);
898  }
899  if (F == G)
900  {
901  coF= coG= Lc (F);
902  return F/Lc(F);
903  }
904 
905  CFMap M,N;
906  int best_level= myCompress (A, B, M, N, topLevel);
907 
908  if (best_level == 0)
909  {
910  coF= F;
911  coG= G;
912  return B.genOne();
913  }
914 
915  A= M(A);
916  B= M(B);
917 
918  Variable x= Variable(1);
919 
920  //univariate case
921  if (A.isUnivariate() && B.isUnivariate())
922  {
923  CanonicalForm result= gcd (A, B);
924  coF= N (A/result);
925  coG= N (B/result);
926  return N (result);
927  }
928 
929  CanonicalForm cA, cB; // content of A and B
930  CanonicalForm ppA, ppB; // primitive part of A and B
931  CanonicalForm gcdcAcB;
932 
933  cA = uni_content (A);
934  cB = uni_content (B);
935  gcdcAcB= gcd (cA, cB);
936  ppA= A/cA;
937  ppB= B/cB;
938 
939  CanonicalForm lcA, lcB; // leading coefficients of A and B
940  CanonicalForm gcdlcAlcB;
941 
942  lcA= uni_lcoeff (ppA);
943  lcB= uni_lcoeff (ppB);
944 
945  gcdlcAlcB= gcd (lcA, lcB);
946 
947  int d= totaldegree (ppA, Variable(2), Variable (ppA.level()));
948  if (d == 0)
949  {
950  coF= N (ppA*(cA/gcdcAcB));
951  coG= N (ppB*(cB/gcdcAcB));
952  return N(gcdcAcB);
953  }
954  int d0= totaldegree (ppB, Variable(2), Variable (ppB.level()));
955  if (d0 < d)
956  d= d0;
957  if (d == 0)
958  {
959  coF= N (ppA*(cA/gcdcAcB));
960  coG= N (ppB*(cB/gcdcAcB));
961  return N(gcdcAcB);
962  }
963 
964  CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH;
965  CanonicalForm newtonPoly, coF_random_element, coG_random_element, coF_m,
966  coG_m, ppCoF, ppCoG;
967 
968  newtonPoly= 1;
969  m= gcdlcAlcB;
970  G_m= 0;
971  coF= 0;
972  coG= 0;
973  H= 0;
974  bool fail= false;
975  topLevel= false;
976  bool inextension= false;
977  int p=-1;
978  int k= getGFDegree();
979  int kk;
980  int expon;
981  char gf_name_buf= gf_name;
982  int bound1= degree (ppA, 1);
983  int bound2= degree (ppB, 1);
984  do
985  {
986  random_element= GFRandomElement (m*lcA*lcB, l, fail);
987  if (fail)
988  {
989  p= getCharacteristic();
990  expon= 2;
991  kk= getGFDegree();
992  if (ipower (p, kk*expon) < (1 << 16))
993  setCharacteristic (p, kk*(int)expon, 'b');
994  else
995  {
996  expon= (int) floor((log ((double)((1<<16) - 1)))/(log((double)p)*kk));
997  ASSERT (expon >= 2, "not enough points in modGCDGF");
998  setCharacteristic (p, (int)(kk*expon), 'b');
999  }
1000  inextension= true;
1001  fail= false;
1002  for (CFListIterator i= l; i.hasItem(); i++)
1003  i.getItem()= GFMapUp (i.getItem(), kk);
1004  m= GFMapUp (m, kk);
1005  G_m= GFMapUp (G_m, kk);
1006  newtonPoly= GFMapUp (newtonPoly, kk);
1007  coF_m= GFMapUp (coF_m, kk);
1008  coG_m= GFMapUp (coG_m, kk);
1009  ppA= GFMapUp (ppA, kk);
1010  ppB= GFMapUp (ppB, kk);
1011  gcdlcAlcB= GFMapUp (gcdlcAlcB, kk);
1012  lcA= GFMapUp (lcA, kk);
1013  lcB= GFMapUp (lcB, kk);
1014  random_element= GFRandomElement (m*lcA*lcB, l, fail);
1015  DEBOUTLN (cerr, "fail= " << fail);
1016  CFList list;
1017  TIMING_START (gcd_recursion);
1018  G_random_element= modGCDGF (ppA(random_element, x), ppB(random_element, x),
1019  coF_random_element, coG_random_element,
1020  list, topLevel);
1021  TIMING_END_AND_PRINT (gcd_recursion,
1022  "time for recursive call: ");
1023  DEBOUTLN (cerr, "G_random_element= " << G_random_element);
1024  }
1025  else
1026  {
1027  CFList list;
1028  TIMING_START (gcd_recursion);
1029  G_random_element= modGCDGF (ppA(random_element, x), ppB(random_element, x),
1030  coF_random_element, coG_random_element,
1031  list, topLevel);
1032  TIMING_END_AND_PRINT (gcd_recursion,
1033  "time for recursive call: ");
1034  DEBOUTLN (cerr, "G_random_element= " << G_random_element);
1035  }
1036 
1037  if (!G_random_element.inCoeffDomain())
1038  d0= totaldegree (G_random_element, Variable(2),
1039  Variable (G_random_element.level()));
1040  else
1041  d0= 0;
1042 
1043  if (d0 == 0)
1044  {
1045  if (inextension)
1046  {
1047  ppA= GFMapDown (ppA, k);
1048  ppB= GFMapDown (ppB, k);
1049  setCharacteristic (p, k, gf_name_buf);
1050  }
1051  coF= N (ppA*(cA/gcdcAcB));
1052  coG= N (ppB*(cB/gcdcAcB));
1053  return N(gcdcAcB);
1054  }
1055  if (d0 > d)
1056  {
1057  if (!find (l, random_element))
1058  l.append (random_element);
1059  continue;
1060  }
1061 
1062  G_random_element=
1063  (gcdlcAlcB(random_element, x)/uni_lcoeff(G_random_element)) *
1064  G_random_element;
1065 
1066  coF_random_element= (lcA(random_element,x)/uni_lcoeff(coF_random_element))
1067  *coF_random_element;
1068  coG_random_element= (lcB(random_element,x)/uni_lcoeff(coG_random_element))
1069  *coG_random_element;
1070 
1071  if (!G_random_element.inCoeffDomain())
1072  d0= totaldegree (G_random_element, Variable(2),
1073  Variable (G_random_element.level()));
1074  else
1075  d0= 0;
1076 
1077  if (d0 < d)
1078  {
1079  m= gcdlcAlcB;
1080  newtonPoly= 1;
1081  G_m= 0;
1082  d= d0;
1083  coF_m= 0;
1084  coG_m= 0;
1085  }
1086 
1087  TIMING_START (newton_interpolation);
1088  H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
1089  coF= newtonInterp (random_element, coF_random_element, newtonPoly, coF_m,x);
1090  coG= newtonInterp (random_element, coG_random_element, newtonPoly, coG_m,x);
1091  TIMING_END_AND_PRINT (newton_interpolation,
1092  "time for newton interpolation: ");
1093 
1094  //termination test
1095  if ((uni_lcoeff (H) == gcdlcAlcB) || (G_m == H))
1096  {
1097  TIMING_START (termination_test);
1098  if (gcdlcAlcB.isOne())
1099  cH= 1;
1100  else
1101  cH= uni_content (H);
1102  ppH= H/cH;
1103  ppH /= Lc (ppH);
1104  CanonicalForm lcppH= gcdlcAlcB/cH;
1105  CanonicalForm ccoF= lcppH/Lc (lcppH);
1106  CanonicalForm ccoG= lcppH/Lc (lcppH);
1107  ppCoF= coF/ccoF;
1108  ppCoG= coG/ccoG;
1109  if (inextension)
1110  {
1111  if (((degree (ppCoF,1)+degree (ppH,1) == bound1) &&
1112  (degree (ppCoG,1)+degree (ppH,1) == bound2) &&
1113  terminationTest (ppA, ppB, ppCoF, ppCoG, ppH)) ||
1114  (fdivides (ppH, ppA, ppCoF) && fdivides (ppH, ppB, ppCoG)))
1115  {
1116  DEBOUTLN (cerr, "ppH before mapDown= " << ppH);
1117  ppH= GFMapDown (ppH, k);
1118  ppCoF= GFMapDown (ppCoF, k);
1119  ppCoG= GFMapDown (ppCoG, k);
1120  DEBOUTLN (cerr, "ppH after mapDown= " << ppH);
1121  coF= N ((cA/gcdcAcB)*ppCoF);
1122  coG= N ((cB/gcdcAcB)*ppCoG);
1123  setCharacteristic (p, k, gf_name_buf);
1124  TIMING_END_AND_PRINT (termination_test,
1125  "time for successful termination GF: ");
1126  return N(gcdcAcB*ppH);
1127  }
1128  }
1129  else
1130  {
1131  if (((degree (ppCoF,1)+degree (ppH,1) == bound1) &&
1132  (degree (ppCoG,1)+degree (ppH,1) == bound2) &&
1133  terminationTest (ppA, ppB, ppCoF, ppCoG, ppH)) ||
1134  (fdivides (ppH, ppA, ppCoF) && fdivides (ppH, ppB, ppCoG)))
1135  {
1136  coF= N ((cA/gcdcAcB)*ppCoF);
1137  coG= N ((cB/gcdcAcB)*ppCoG);
1138  TIMING_END_AND_PRINT (termination_test,
1139  "time for successful termination GF: ");
1140  return N(gcdcAcB*ppH);
1141  }
1142  }
1143  TIMING_END_AND_PRINT (termination_test,
1144  "time for unsuccessful termination GF: ");
1145  }
1146 
1147  G_m= H;
1148  coF_m= coF;
1149  coG_m= coG;
1150  newtonPoly= newtonPoly*(x - random_element);
1151  m= m*(x - random_element);
1152  if (!find (l, random_element))
1153  l.append (random_element);
1154  } while (1);
1155 }
1156 
1157 static inline
1159 FpRandomElement (const CanonicalForm& F, CFList& list, bool& fail)
1160 {
1161  fail= false;
1162  Variable x= F.mvar();
1163  FFRandom genFF;
1164  CanonicalForm random;
1165  int p= getCharacteristic();
1166  int bound= p;
1167  do
1168  {
1169  if (list.length() == bound)
1170  {
1171  fail= true;
1172  break;
1173  }
1174  if (list.length() < 1)
1175  random= 0;
1176  else
1177  {
1178  random= genFF.generate();
1179  while (find (list, random))
1180  random= genFF.generate();
1181  }
1182  if (F (random, x) == 0)
1183  {
1184  list.append (random);
1185  continue;
1186  }
1187  } while (find (list, random));
1188  return random;
1189 }
1190 
1192 modGCDFp (const CanonicalForm& F, const CanonicalForm& G,
1193  CanonicalForm& coF, CanonicalForm& coG,
1194  bool& topLevel, CFList& l);
1195 
1198  bool& topLevel, CFList& l)
1199 {
1200  CanonicalForm dummy1, dummy2;
1201  CanonicalForm result= modGCDFp (F, G, dummy1, dummy2, topLevel, l);
1202  return result;
1203 }
1204 
1207  CanonicalForm& coF, CanonicalForm& coG,
1208  bool& topLevel, CFList& l)
1209 {
1210  CanonicalForm A= F;
1211  CanonicalForm B= G;
1212  if (F.isZero() && degree(G) > 0)
1213  {
1214  coF= 0;
1215  coG= Lc (G);
1216  return G/Lc(G);
1217  }
1218  else if (G.isZero() && degree (F) > 0)
1219  {
1220  coF= Lc (F);
1221  coG= 0;
1222  return F/Lc(F);
1223  }
1224  else if (F.isZero() && G.isZero())
1225  {
1226  coF= coG= 0;
1227  return F.genOne();
1228  }
1229  if (F.inBaseDomain() || G.inBaseDomain())
1230  {
1231  coF= F;
1232  coG= G;
1233  return F.genOne();
1234  }
1235  if (F.isUnivariate() && fdivides(F, G, coG))
1236  {
1237  coF= Lc (F);
1238  return F/Lc(F);
1239  }
1240  if (G.isUnivariate() && fdivides(G, F, coF))
1241  {
1242  coG= Lc (G);
1243  return G/Lc(G);
1244  }
1245  if (F == G)
1246  {
1247  coF= coG= Lc (F);
1248  return F/Lc(F);
1249  }
1250  CFMap M,N;
1251  int best_level= myCompress (A, B, M, N, topLevel);
1252 
1253  if (best_level == 0)
1254  {
1255  coF= F;
1256  coG= G;
1257  return B.genOne();
1258  }
1259 
1260  A= M(A);
1261  B= M(B);
1262 
1263  Variable x= Variable (1);
1264 
1265  //univariate case
1266  if (A.isUnivariate() && B.isUnivariate())
1267  {
1268  CanonicalForm result= gcd (A, B);
1269  coF= N (A/result);
1270  coG= N (B/result);
1271  return N (result);
1272  }
1273 
1274  CanonicalForm cA, cB; // content of A and B
1275  CanonicalForm ppA, ppB; // primitive part of A and B
1276  CanonicalForm gcdcAcB;
1277 
1278  cA = uni_content (A);
1279  cB = uni_content (B);
1280  gcdcAcB= gcd (cA, cB);
1281  ppA= A/cA;
1282  ppB= B/cB;
1283 
1284  CanonicalForm lcA, lcB; // leading coefficients of A and B
1285  CanonicalForm gcdlcAlcB;
1286  lcA= uni_lcoeff (ppA);
1287  lcB= uni_lcoeff (ppB);
1288 
1289  gcdlcAlcB= gcd (lcA, lcB);
1290 
1291  int d= totaldegree (ppA, Variable (2), Variable (ppA.level()));
1292  int d0;
1293 
1294  if (d == 0)
1295  {
1296  coF= N (ppA*(cA/gcdcAcB));
1297  coG= N (ppB*(cB/gcdcAcB));
1298  return N(gcdcAcB);
1299  }
1300 
1301  d0= totaldegree (ppB, Variable (2), Variable (ppB.level()));
1302 
1303  if (d0 < d)
1304  d= d0;
1305 
1306  if (d == 0)
1307  {
1308  coF= N (ppA*(cA/gcdcAcB));
1309  coG= N (ppB*(cB/gcdcAcB));
1310  return N(gcdcAcB);
1311  }
1312 
1313  CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH;
1314  CanonicalForm newtonPoly, coF_random_element, coG_random_element,
1315  coF_m, coG_m, ppCoF, ppCoG;
1316 
1317  newtonPoly= 1;
1318  m= gcdlcAlcB;
1319  H= 0;
1320  coF= 0;
1321  coG= 0;
1322  G_m= 0;
1323  Variable alpha, V_buf, cleanUp;
1324  bool fail= false;
1325  bool inextension= false;
1326  topLevel= false;
1327  CFList source, dest;
1328  int bound1= degree (ppA, 1);
1329  int bound2= degree (ppB, 1);
1330  do
1331  {
1332  if (inextension)
1333  random_element= randomElement (m*lcA*lcB, V_buf, l, fail);
1334  else
1335  random_element= FpRandomElement (m*lcA*lcB, l, fail);
1336 
1337  if (!fail && !inextension)
1338  {
1339  CFList list;
1340  TIMING_START (gcd_recursion);
1341  G_random_element=
1342  modGCDFp (ppA (random_element,x), ppB (random_element,x),
1343  coF_random_element, coG_random_element, topLevel,
1344  list);
1345  TIMING_END_AND_PRINT (gcd_recursion,
1346  "time for recursive call: ");
1347  DEBOUTLN (cerr, "G_random_element= " << G_random_element);
1348  }
1349  else if (!fail && inextension)
1350  {
1351  CFList list;
1352  TIMING_START (gcd_recursion);
1353  G_random_element=
1354  modGCDFq (ppA (random_element, x), ppB (random_element, x),
1355  coF_random_element, coG_random_element, V_buf,
1356  list, topLevel);
1357  TIMING_END_AND_PRINT (gcd_recursion,
1358  "time for recursive call: ");
1359  DEBOUTLN (cerr, "G_random_element= " << G_random_element);
1360  }
1361  else if (fail && !inextension)
1362  {
1363  source= CFList();
1364  dest= CFList();
1365  CFList list;
1367  int deg= 2;
1368  bool initialized= false;
1369  do
1370  {
1371  mipo= randomIrredpoly (deg, x);
1372  if (initialized)
1373  setMipo (alpha, mipo);
1374  else
1375  alpha= rootOf (mipo);
1376  inextension= true;
1377  initialized= true;
1378  fail= false;
1379  random_element= randomElement (m*lcA*lcB, alpha, l, fail);
1380  deg++;
1381  } while (fail);
1382  list= CFList();
1383  V_buf= alpha;
1384  cleanUp= alpha;
1385  TIMING_START (gcd_recursion);
1386  G_random_element=
1387  modGCDFq (ppA (random_element, x), ppB (random_element, x),
1388  coF_random_element, coG_random_element, alpha,
1389  list, topLevel);
1390  TIMING_END_AND_PRINT (gcd_recursion,
1391  "time for recursive call: ");
1392  DEBOUTLN (cerr, "G_random_element= " << G_random_element);
1393  }
1394  else if (fail && inextension)
1395  {
1396  source= CFList();
1397  dest= CFList();
1398 
1399  Variable V_buf3= V_buf;
1400  V_buf= chooseExtension (V_buf);
1401  bool prim_fail= false;
1402  Variable V_buf2;
1403  CanonicalForm prim_elem, im_prim_elem;
1404  prim_elem= primitiveElement (alpha, V_buf2, prim_fail);
1405 
1406  if (V_buf3 != alpha)
1407  {
1408  m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
1409  G_m= mapDown (G_m, prim_elem, im_prim_elem, alpha, source, dest);
1410  coF_m= mapDown (coF_m, prim_elem, im_prim_elem, alpha, source, dest);
1411  coG_m= mapDown (coG_m, prim_elem, im_prim_elem, alpha, source, dest);
1412  newtonPoly= mapDown (newtonPoly, prim_elem, im_prim_elem, alpha,
1413  source, dest);
1414  ppA= mapDown (ppA, prim_elem, im_prim_elem, alpha, source, dest);
1415  ppB= mapDown (ppB, prim_elem, im_prim_elem, alpha, source, dest);
1416  gcdlcAlcB= mapDown (gcdlcAlcB, prim_elem, im_prim_elem, alpha, source,
1417  dest);
1418  lcA= mapDown (lcA, prim_elem, im_prim_elem, alpha, source, dest);
1419  lcB= mapDown (lcB, prim_elem, im_prim_elem, alpha, source, dest);
1420  for (CFListIterator i= l; i.hasItem(); i++)
1421  i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, alpha,
1422  source, dest);
1423  }
1424 
1425  ASSERT (!prim_fail, "failure in integer factorizer");
1426  if (prim_fail)
1427  ; //ERROR
1428  else
1429  im_prim_elem= mapPrimElem (prim_elem, alpha, V_buf);
1430 
1431  DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (alpha));
1432  DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf2));
1433 
1434  for (CFListIterator i= l; i.hasItem(); i++)
1435  i.getItem()= mapUp (i.getItem(), alpha, V_buf, prim_elem,
1436  im_prim_elem, source, dest);
1437  m= mapUp (m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
1438  G_m= mapUp (G_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
1439  coF_m= mapUp (coF_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
1440  coG_m= mapUp (coG_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
1441  newtonPoly= mapUp (newtonPoly, alpha, V_buf, prim_elem, im_prim_elem,
1442  source, dest);
1443  ppA= mapUp (ppA, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
1444  ppB= mapUp (ppB, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
1445  gcdlcAlcB= mapUp (gcdlcAlcB, alpha, V_buf, prim_elem, im_prim_elem,
1446  source, dest);
1447  lcA= mapUp (lcA, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
1448  lcB= mapUp (lcB, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
1449  fail= false;
1450  random_element= randomElement (m*lcA*lcB, V_buf, l, fail );
1451  DEBOUTLN (cerr, "fail= " << fail);
1452  CFList list;
1453  TIMING_START (gcd_recursion);
1454  G_random_element=
1455  modGCDFq (ppA (random_element, x), ppB (random_element, x),
1456  coF_random_element, coG_random_element, V_buf,
1457  list, topLevel);
1458  TIMING_END_AND_PRINT (gcd_recursion,
1459  "time for recursive call: ");
1460  DEBOUTLN (cerr, "G_random_element= " << G_random_element);
1461  alpha= V_buf;
1462  }
1463 
1464  if (!G_random_element.inCoeffDomain())
1465  d0= totaldegree (G_random_element, Variable(2),
1466  Variable (G_random_element.level()));
1467  else
1468  d0= 0;
1469 
1470  if (d0 == 0)
1471  {
1472  if (inextension)
1473  prune (cleanUp);
1474  coF= N (ppA*(cA/gcdcAcB));
1475  coG= N (ppB*(cB/gcdcAcB));
1476  return N(gcdcAcB);
1477  }
1478 
1479  if (d0 > d)
1480  {
1481  if (!find (l, random_element))
1482  l.append (random_element);
1483  continue;
1484  }
1485 
1486  G_random_element= (gcdlcAlcB(random_element,x)/uni_lcoeff(G_random_element))
1487  *G_random_element;
1488 
1489  coF_random_element= (lcA(random_element,x)/uni_lcoeff(coF_random_element))
1490  *coF_random_element;
1491  coG_random_element= (lcB(random_element,x)/uni_lcoeff(coG_random_element))
1492  *coG_random_element;
1493 
1494  if (!G_random_element.inCoeffDomain())
1495  d0= totaldegree (G_random_element, Variable(2),
1496  Variable (G_random_element.level()));
1497  else
1498  d0= 0;
1499 
1500  if (d0 < d)
1501  {
1502  m= gcdlcAlcB;
1503  newtonPoly= 1;
1504  G_m= 0;
1505  d= d0;
1506  coF_m= 0;
1507  coG_m= 0;
1508  }
1509 
1510  TIMING_START (newton_interpolation);
1511  H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
1512  coF= newtonInterp (random_element, coF_random_element, newtonPoly, coF_m,x);
1513  coG= newtonInterp (random_element, coG_random_element, newtonPoly, coG_m,x);
1514  TIMING_END_AND_PRINT (newton_interpolation,
1515  "time for newton_interpolation: ");
1516 
1517  //termination test
1518  if ((uni_lcoeff (H) == gcdlcAlcB) || (G_m == H))
1519  {
1520  TIMING_START (termination_test);
1521  if (gcdlcAlcB.isOne())
1522  cH= 1;
1523  else
1524  cH= uni_content (H);
1525  ppH= H/cH;
1526  ppH /= Lc (ppH);
1527  CanonicalForm lcppH= gcdlcAlcB/cH;
1528  CanonicalForm ccoF= lcppH/Lc (lcppH);
1529  CanonicalForm ccoG= lcppH/Lc (lcppH);
1530  ppCoF= coF/ccoF;
1531  ppCoG= coG/ccoG;
1532  DEBOUTLN (cerr, "ppH= " << ppH);
1533  if (((degree (ppCoF,1)+degree (ppH,1) == bound1) &&
1534  (degree (ppCoG,1)+degree (ppH,1) == bound2) &&
1535  terminationTest (ppA, ppB, ppCoF, ppCoG, ppH)) ||
1536  (fdivides (ppH, ppA, ppCoF) && fdivides (ppH, ppB, ppCoG)))
1537  {
1538  if (inextension)
1539  prune (cleanUp);
1540  coF= N ((cA/gcdcAcB)*ppCoF);
1541  coG= N ((cB/gcdcAcB)*ppCoG);
1542  TIMING_END_AND_PRINT (termination_test,
1543  "time for successful termination Fp: ");
1544  return N(gcdcAcB*ppH);
1545  }
1546  TIMING_END_AND_PRINT (termination_test,
1547  "time for unsuccessful termination Fp: ");
1548  }
1549 
1550  G_m= H;
1551  coF_m= coF;
1552  coG_m= coG;
1553  newtonPoly= newtonPoly*(x - random_element);
1554  m= m*(x - random_element);
1555  if (!find (l, random_element))
1556  l.append (random_element);
1557  } while (1);
1558 }
1559 
1560 CFArray
1562 {
1563  int r= M.size();
1564  ASSERT (A.size() == r, "vector does not have right size");
1565 
1566  if (r == 1)
1567  {
1568  CFArray result= CFArray (1);
1569  result [0]= A [0] / M [0];
1570  return result;
1571  }
1572  // check solvability
1573  bool notDistinct= false;
1574  for (int i= 0; i < r - 1; i++)
1575  {
1576  for (int j= i + 1; j < r; j++)
1577  {
1578  if (M [i] == M [j])
1579  {
1580  notDistinct= true;
1581  break;
1582  }
1583  }
1584  }
1585  if (notDistinct)
1586  return CFArray();
1587 
1588  CanonicalForm master= 1;
1589  Variable x= Variable (1);
1590  for (int i= 0; i < r; i++)
1591  master *= x - M [i];
1592  CFList Pj;
1593  CanonicalForm tmp;
1594  for (int i= 0; i < r; i++)
1595  {
1596  tmp= master/(x - M [i]);
1597  tmp /= tmp (M [i], 1);
1598  Pj.append (tmp);
1599  }
1600  CFArray result= CFArray (r);
1601 
1602  CFListIterator j= Pj;
1603  for (int i= 1; i <= r; i++, j++)
1604  {
1605  tmp= 0;
1606  for (int l= 0; l < A.size(); l++)
1607  tmp += A[l]*j.getItem()[l];
1608  result[i - 1]= tmp;
1609  }
1610  return result;
1611 }
1612 
1613 CFArray
1615 {
1616  int r= M.size();
1617  ASSERT (A.size() == r, "vector does not have right size");
1618  if (r == 1)
1619  {
1620  CFArray result= CFArray (1);
1621  result [0]= A[0] / M [0];
1622  return result;
1623  }
1624  // check solvability
1625  bool notDistinct= false;
1626  for (int i= 0; i < r - 1; i++)
1627  {
1628  for (int j= i + 1; j < r; j++)
1629  {
1630  if (M [i] == M [j])
1631  {
1632  notDistinct= true;
1633  break;
1634  }
1635  }
1636  }
1637  if (notDistinct)
1638  return CFArray();
1639 
1640  CanonicalForm master= 1;
1641  Variable x= Variable (1);
1642  for (int i= 0; i < r; i++)
1643  master *= x - M [i];
1644  master *= x;
1645  CFList Pj;
1646  CanonicalForm tmp;
1647  for (int i= 0; i < r; i++)
1648  {
1649  tmp= master/(x - M [i]);
1650  tmp /= tmp (M [i], 1);
1651  Pj.append (tmp);
1652  }
1653 
1654  CFArray result= CFArray (r);
1655 
1656  CFListIterator j= Pj;
1657  for (int i= 1; i <= r; i++, j++)
1658  {
1659  tmp= 0;
1660 
1661  for (int l= 1; l <= A.size(); l++)
1662  tmp += A[l - 1]*j.getItem()[l];
1663  result[i - 1]= tmp;
1664  }
1665  return result;
1666 }
1667 
1668 /// M in row echolon form, rk rank of M
1669 CFArray
1670 readOffSolution (const CFMatrix& M, const long rk)
1671 {
1672  CFArray result= CFArray (rk);
1673  CanonicalForm tmp1, tmp2, tmp3;
1674  for (int i= rk; i >= 1; i--)
1675  {
1676  tmp3= 0;
1677  tmp1= M (i, M.columns());
1678  for (int j= M.columns() - 1; j >= 1; j--)
1679  {
1680  tmp2= M (i, j);
1681  if (j == i)
1682  break;
1683  else
1684  tmp3 += tmp2*result[j - 1];
1685  }
1686  result[i - 1]= (tmp1 - tmp3)/tmp2;
1687  }
1688  return result;
1689 }
1690 
1691 CFArray
1692 readOffSolution (const CFMatrix& M, const CFArray& L, const CFArray& partialSol)
1693 {
1694  CFArray result= CFArray (M.rows());
1695  CanonicalForm tmp1, tmp2, tmp3;
1696  int k;
1697  for (int i= M.rows(); i >= 1; i--)
1698  {
1699  tmp3= 0;
1700  tmp1= L[i - 1];
1701  k= 0;
1702  for (int j= M.columns(); j >= 1; j--, k++)
1703  {
1704  tmp2= M (i, j);
1705  if (j == i)
1706  break;
1707  else
1708  {
1709  if (k > partialSol.size() - 1)
1710  tmp3 += tmp2*result[j - 1];
1711  else
1712  tmp3 += tmp2*partialSol[partialSol.size() - k - 1];
1713  }
1714  }
1715  result [i - 1]= (tmp1 - tmp3)/tmp2;
1716  }
1717  return result;
1718 }
1719 
1720 long
1722 {
1723  ASSERT (L.size() <= M.rows(), "dimension exceeded");
1724  CFMatrix *N;
1725  N= new CFMatrix (M.rows(), M.columns() + 1);
1726 
1727  for (int i= 1; i <= M.rows(); i++)
1728  for (int j= 1; j <= M.columns(); j++)
1729  (*N) (i, j)= M (i, j);
1730 
1731  int j= 1;
1732  for (int i= 0; i < L.size(); i++, j++)
1733  (*N) (j, M.columns() + 1)= L[i];
1734 #ifdef HAVE_FLINT
1735  nmod_mat_t FLINTN;
1736  convertFacCFMatrix2nmod_mat_t (FLINTN, *N);
1737  long rk= nmod_mat_rref (FLINTN);
1738 
1739  delete N;
1740  N= convertNmod_mat_t2FacCFMatrix (FLINTN);
1741  nmod_mat_clear (FLINTN);
1742 #else
1743  int p= getCharacteristic ();
1744  if (fac_NTL_char != p)
1745  {
1746  fac_NTL_char= p;
1747  zz_p::init (p);
1748  }
1749  mat_zz_p *NTLN= convertFacCFMatrix2NTLmat_zz_p(*N);
1750  delete N;
1751  long rk= gauss (*NTLN);
1752 
1753  N= convertNTLmat_zz_p2FacCFMatrix (*NTLN);
1754  delete NTLN;
1755 #endif
1756 
1757  L= CFArray (M.rows());
1758  for (int i= 0; i < M.rows(); i++)
1759  L[i]= (*N) (i + 1, M.columns() + 1);
1760  M= (*N) (1, M.rows(), 1, M.columns());
1761  delete N;
1762  return rk;
1763 }
1764 
1765 long
1767 {
1768  ASSERT (L.size() <= M.rows(), "dimension exceeded");
1769  CFMatrix *N;
1770  N= new CFMatrix (M.rows(), M.columns() + 1);
1771 
1772  for (int i= 1; i <= M.rows(); i++)
1773  for (int j= 1; j <= M.columns(); j++)
1774  (*N) (i, j)= M (i, j);
1775 
1776  int j= 1;
1777  for (int i= 0; i < L.size(); i++, j++)
1778  (*N) (j, M.columns() + 1)= L[i];
1779  int p= getCharacteristic ();
1780  if (fac_NTL_char != p)
1781  {
1782  fac_NTL_char= p;
1783  zz_p::init (p);
1784  }
1785  zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
1786  zz_pE::init (NTLMipo);
1787  mat_zz_pE *NTLN= convertFacCFMatrix2NTLmat_zz_pE(*N);
1788  long rk= gauss (*NTLN);
1789 
1790  delete N;
1791  N= convertNTLmat_zz_pE2FacCFMatrix (*NTLN, alpha);
1792 
1793  delete NTLN;
1794 
1795  M= (*N) (1, M.rows(), 1, M.columns());
1796  L= CFArray (M.rows());
1797  for (int i= 0; i < M.rows(); i++)
1798  L[i]= (*N) (i + 1, M.columns() + 1);
1799 
1800  delete N;
1801  return rk;
1802 }
1803 
1804 CFArray
1805 solveSystemFp (const CFMatrix& M, const CFArray& L)
1806 {
1807  ASSERT (L.size() <= M.rows(), "dimension exceeded");
1808  CFMatrix *N;
1809  N= new CFMatrix (M.rows(), M.columns() + 1);
1810 
1811  for (int i= 1; i <= M.rows(); i++)
1812  for (int j= 1; j <= M.columns(); j++)
1813  (*N) (i, j)= M (i, j);
1814 
1815  int j= 1;
1816  for (int i= 0; i < L.size(); i++, j++)
1817  (*N) (j, M.columns() + 1)= L[i];
1818 
1819 #ifdef HAVE_FLINT
1820  nmod_mat_t FLINTN;
1821  convertFacCFMatrix2nmod_mat_t (FLINTN, *N);
1822  long rk= nmod_mat_rref (FLINTN);
1823 #else
1824  int p= getCharacteristic ();
1825  if (fac_NTL_char != p)
1826  {
1827  fac_NTL_char= p;
1828  zz_p::init (p);
1829  }
1830  mat_zz_p *NTLN= convertFacCFMatrix2NTLmat_zz_p(*N);
1831  long rk= gauss (*NTLN);
1832 #endif
1833  delete N;
1834  if (rk != M.columns())
1835  {
1836 #ifdef HAVE_FLINT
1837  nmod_mat_clear (FLINTN);
1838 #else
1839  delete NTLN;
1840 #endif
1841  return CFArray();
1842  }
1843 #ifdef HAVE_FLINT
1844  N= convertNmod_mat_t2FacCFMatrix (FLINTN);
1845  nmod_mat_clear (FLINTN);
1846 #else
1847  N= convertNTLmat_zz_p2FacCFMatrix (*NTLN);
1848  delete NTLN;
1849 #endif
1850  CFArray A= readOffSolution (*N, rk);
1851 
1852  delete N;
1853  return A;
1854 }
1855 
1856 CFArray
1857 solveSystemFq (const CFMatrix& M, const CFArray& L, const Variable& alpha)
1858 {
1859  ASSERT (L.size() <= M.rows(), "dimension exceeded");
1860  CFMatrix *N;
1861  N= new CFMatrix (M.rows(), M.columns() + 1);
1862 
1863  for (int i= 1; i <= M.rows(); i++)
1864  for (int j= 1; j <= M.columns(); j++)
1865  (*N) (i, j)= M (i, j);
1866  int j= 1;
1867  for (int i= 0; i < L.size(); i++, j++)
1868  (*N) (j, M.columns() + 1)= L[i];
1869  int p= getCharacteristic ();
1870  if (fac_NTL_char != p)
1871  {
1872  fac_NTL_char= p;
1873  zz_p::init (p);
1874  }
1875  zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
1876  zz_pE::init (NTLMipo);
1877  mat_zz_pE *NTLN= convertFacCFMatrix2NTLmat_zz_pE(*N);
1878  long rk= gauss (*NTLN);
1879 
1880  delete N;
1881  if (rk != M.columns())
1882  {
1883  delete NTLN;
1884  return CFArray();
1885  }
1886  N= convertNTLmat_zz_pE2FacCFMatrix (*NTLN, alpha);
1887 
1888  delete NTLN;
1889 
1890  CFArray A= readOffSolution (*N, rk);
1891 
1892  delete N;
1893  return A;
1894 }
1895 #endif
1896 
1897 CFArray
1899 {
1900  if (F.inCoeffDomain())
1901  {
1902  CFArray result= CFArray (1);
1903  result [0]= 1;
1904  return result;
1905  }
1906  if (F.isUnivariate())
1907  {
1908  CFArray result= CFArray (size(F));
1909  int j= 0;
1910  for (CFIterator i= F; i.hasTerms(); i++, j++)
1911  result[j]= power (F.mvar(), i.exp());
1912  return result;
1913  }
1914  int numMon= size (F);
1915  CFArray result= CFArray (numMon);
1916  int j= 0;
1917  CFArray recResult;
1918  Variable x= F.mvar();
1919  CanonicalForm powX;
1920  for (CFIterator i= F; i.hasTerms(); i++)
1921  {
1922  powX= power (x, i.exp());
1923  recResult= getMonoms (i.coeff());
1924  for (int k= 0; k < recResult.size(); k++)
1925  result[j+k]= powX*recResult[k];
1926  j += recResult.size();
1927  }
1928  return result;
1929 }
1930 
1931 #ifdef HAVE_NTL
1932 CFArray
1934 {
1935  if (F.inCoeffDomain())
1936  {
1937  CFArray result= CFArray (1);
1938  result [0]= F;
1939  return result;
1940  }
1941  if (F.isUnivariate())
1942  {
1943  ASSERT (evalPoints.length() == 1,
1944  "expected an eval point with only one component");
1945  CFArray result= CFArray (size(F));
1946  int j= 0;
1947  CanonicalForm evalPoint= evalPoints.getLast();
1948  for (CFIterator i= F; i.hasTerms(); i++, j++)
1949  result[j]= power (evalPoint, i.exp());
1950  return result;
1951  }
1952  int numMon= size (F);
1953  CFArray result= CFArray (numMon);
1954  int j= 0;
1955  CanonicalForm evalPoint= evalPoints.getLast();
1957  buf.removeLast();
1958  CFArray recResult;
1959  CanonicalForm powEvalPoint;
1960  for (CFIterator i= F; i.hasTerms(); i++)
1961  {
1962  powEvalPoint= power (evalPoint, i.exp());
1963  recResult= evaluateMonom (i.coeff(), buf);
1964  for (int k= 0; k < recResult.size(); k++)
1965  result[j+k]= powEvalPoint*recResult[k];
1966  j += recResult.size();
1967  }
1968  return result;
1969 }
1970 
1971 CFArray
1973 {
1974  CFArray result= A.size();
1975  CanonicalForm tmp;
1976  int k;
1977  for (int i= 0; i < A.size(); i++)
1978  {
1979  tmp= A[i];
1980  k= 1;
1981  for (CFListIterator j= evalPoints; j.hasItem(); j++, k++)
1982  tmp= tmp (j.getItem(), k);
1983  result[i]= tmp;
1984  }
1985  return result;
1986 }
1987 
1988 CFList
1991  const CanonicalForm& LCF, const bool& GF,
1992  const Variable& alpha, bool& fail, CFList& list
1993  )
1994 {
1995  int k= tmax (F.level(), G.level()) - 1;
1996  Variable x= Variable (1);
1997  CFList result;
1998  FFRandom genFF;
1999  GFRandom genGF;
2000  int p= getCharacteristic ();
2001  double bound;
2002  if (alpha != Variable (1))
2003  {
2004  bound= pow ((double) p, (double) degree (getMipo(alpha)));
2005  bound= pow (bound, (double) k);
2006  }
2007  else if (GF)
2008  {
2009  bound= pow ((double) p, (double) getGFDegree());
2010  bound= pow ((double) bound, (double) k);
2011  }
2012  else
2013  bound= pow ((double) p, (double) k);
2014 
2015  CanonicalForm random;
2016  int j;
2017  bool zeroOneOccured= false;
2018  bool allEqual= false;
2020  do
2021  {
2022  random= 0;
2023  // possible overflow if list.length() does not fit into a int
2024  if (list.length() >= bound)
2025  {
2026  fail= true;
2027  break;
2028  }
2029  for (int i= 0; i < k; i++)
2030  {
2031  if (GF)
2032  {
2033  result.append (genGF.generate());
2034  random += result.getLast()*power (x, i);
2035  }
2036  else if (alpha.level() != 1)
2037  {
2038  AlgExtRandomF genAlgExt (alpha);
2039  result.append (genAlgExt.generate());
2040  random += result.getLast()*power (x, i);
2041  }
2042  else
2043  {
2044  result.append (genFF.generate());
2045  random += result.getLast()*power (x, i);
2046  }
2047  if (result.getLast().isOne() || result.getLast().isZero())
2048  zeroOneOccured= true;
2049  }
2050  if (find (list, random))
2051  {
2052  zeroOneOccured= false;
2053  allEqual= false;
2054  result= CFList();
2055  continue;
2056  }
2057  if (zeroOneOccured)
2058  {
2059  list.append (random);
2060  zeroOneOccured= false;
2061  allEqual= false;
2062  result= CFList();
2063  continue;
2064  }
2065  // no zero at this point
2066  if (k > 1)
2067  {
2068  allEqual= true;
2069  CFIterator iter= random;
2070  buf= iter.coeff();
2071  iter++;
2072  for (; iter.hasTerms(); iter++)
2073  if (buf != iter.coeff())
2074  allEqual= false;
2075  }
2076  if (allEqual)
2077  {
2078  list.append (random);
2079  allEqual= false;
2080  zeroOneOccured= false;
2081  result= CFList();
2082  continue;
2083  }
2084 
2085  Feval= F;
2086  Geval= G;
2087  CanonicalForm LCeval= LCF;
2088  j= 1;
2089  for (CFListIterator i= result; i.hasItem(); i++, j++)
2090  {
2091  Feval= Feval (i.getItem(), j);
2092  Geval= Geval (i.getItem(), j);
2093  LCeval= LCeval (i.getItem(), j);
2094  }
2095 
2096  if (LCeval.isZero())
2097  {
2098  if (!find (list, random))
2099  list.append (random);
2100  zeroOneOccured= false;
2101  allEqual= false;
2102  result= CFList();
2103  continue;
2104  }
2105 
2106  if (list.length() >= bound)
2107  {
2108  fail= true;
2109  break;
2110  }
2111  } while (find (list, random));
2112 
2113  return result;
2114 }
2115 
2116 /// multiply two lists componentwise
2117 void mult (CFList& L1, const CFList& L2)
2118 {
2119  ASSERT (L1.length() == L2.length(), "lists of the same size expected");
2120 
2121  CFListIterator j= L2;
2122  for (CFListIterator i= L1; i.hasItem(); i++, j++)
2123  i.getItem() *= j.getItem();
2124 }
2125 
2127  CanonicalForm& Beval, const CFList& L)
2128 {
2129  Aeval= A;
2130  Beval= B;
2131  int j= 1;
2132  for (CFListIterator i= L; i.hasItem(); i++, j++)
2133  {
2134  Aeval= Aeval (i.getItem(), j);
2135  Beval= Beval (i.getItem(), j);
2136  }
2137 }
2138 
2141  const CanonicalForm& skeleton, const Variable& alpha,
2142  bool& fail, CFArray*& coeffMonoms, CFArray& Monoms
2143  )
2144 {
2145  CanonicalForm A= F;
2146  CanonicalForm B= G;
2147  if (F.isZero() && degree(G) > 0) return G/Lc(G);
2148  else if (G.isZero() && degree (F) > 0) return F/Lc(F);
2149  else if (F.isZero() && G.isZero()) return F.genOne();
2150  if (F.inBaseDomain() || G.inBaseDomain()) return F.genOne();
2151  if (F.isUnivariate() && fdivides(F, G)) return F/Lc(F);
2152  if (G.isUnivariate() && fdivides(G, F)) return G/Lc(G);
2153  if (F == G) return F/Lc(F);
2154 
2155  ASSERT (degree (A, 1) == 0, "expected degree (F, 1) == 0");
2156  ASSERT (degree (B, 1) == 0, "expected degree (G, 1) == 0");
2157 
2158  CFMap M,N;
2159  int best_level= myCompress (A, B, M, N, false);
2160 
2161  if (best_level == 0)
2162  return B.genOne();
2163 
2164  A= M(A);
2165  B= M(B);
2166 
2167  Variable x= Variable (1);
2168 
2169  //univariate case
2170  if (A.isUnivariate() && B.isUnivariate())
2171  return N (gcd (A, B));
2172 
2173  CanonicalForm skel= M(skeleton);
2174  CanonicalForm cA, cB; // content of A and B
2175  CanonicalForm ppA, ppB; // primitive part of A and B
2176  CanonicalForm gcdcAcB;
2177  cA = uni_content (A);
2178  cB = uni_content (B);
2179  gcdcAcB= gcd (cA, cB);
2180  ppA= A/cA;
2181  ppB= B/cB;
2182 
2183  CanonicalForm lcA, lcB; // leading coefficients of A and B
2184  CanonicalForm gcdlcAlcB;
2185  lcA= uni_lcoeff (ppA);
2186  lcB= uni_lcoeff (ppB);
2187 
2188  if (fdivides (lcA, lcB))
2189  {
2190  if (fdivides (A, B))
2191  return F/Lc(F);
2192  }
2193  if (fdivides (lcB, lcA))
2194  {
2195  if (fdivides (B, A))
2196  return G/Lc(G);
2197  }
2198 
2199  gcdlcAlcB= gcd (lcA, lcB);
2200  int skelSize= size (skel, skel.mvar());
2201 
2202  int j= 0;
2203  int biggestSize= 0;
2204 
2205  for (CFIterator i= skel; i.hasTerms(); i++, j++)
2206  biggestSize= tmax (biggestSize, size (i.coeff()));
2207 
2208  CanonicalForm g, Aeval, Beval;
2209 
2211  bool evalFail= false;
2212  CFList list;
2213  bool GF= false;
2214  CanonicalForm LCA= LC (A);
2215  CanonicalForm tmp;
2216  CFArray gcds= CFArray (biggestSize);
2217  CFList * pEvalPoints= new CFList [biggestSize];
2218  Variable V_buf= alpha, V_buf4= alpha;
2219  CFList source, dest;
2220  CanonicalForm prim_elem, im_prim_elem;
2221  CanonicalForm prim_elem_alpha, im_prim_elem_alpha;
2222  for (int i= 0; i < biggestSize; i++)
2223  {
2224  if (i == 0)
2225  evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf, evalFail,
2226  list);
2227  else
2228  {
2229  mult (evalPoints, pEvalPoints [0]);
2230  eval (A, B, Aeval, Beval, evalPoints);
2231  }
2232 
2233  if (evalFail)
2234  {
2235  if (V_buf.level() != 1)
2236  {
2237  do
2238  {
2239  Variable V_buf3= V_buf;
2240  V_buf= chooseExtension (V_buf);
2241  source= CFList();
2242  dest= CFList();
2243 
2244  bool prim_fail= false;
2245  Variable V_buf2;
2246  prim_elem= primitiveElement (V_buf4, V_buf2, prim_fail);
2247  if (V_buf4 == alpha && alpha.level() != 1)
2248  prim_elem_alpha= prim_elem;
2249 
2250  ASSERT (!prim_fail, "failure in integer factorizer");
2251  if (prim_fail)
2252  ; //ERROR
2253  else
2254  im_prim_elem= mapPrimElem (prim_elem, V_buf4, V_buf);
2255 
2256  DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf));
2257  DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf2));
2258 
2259  if (V_buf4 == alpha && alpha.level() != 1)
2260  im_prim_elem_alpha= im_prim_elem;
2261  else if (alpha.level() != 1)
2262  im_prim_elem_alpha= mapUp (im_prim_elem_alpha, V_buf4, V_buf,
2263  prim_elem, im_prim_elem, source, dest);
2264 
2265  for (CFListIterator j= list; j.hasItem(); j++)
2266  j.getItem()= mapUp (j.getItem(), V_buf4, V_buf, prim_elem,
2267  im_prim_elem, source, dest);
2268  for (int k= 0; k < i; k++)
2269  {
2270  for (CFListIterator j= pEvalPoints[k]; j.hasItem(); j++)
2271  j.getItem()= mapUp (j.getItem(), V_buf4, V_buf, prim_elem,
2272  im_prim_elem, source, dest);
2273  gcds[k]= mapUp (gcds[k], V_buf4, V_buf, prim_elem, im_prim_elem,
2274  source, dest);
2275  }
2276 
2277  if (alpha.level() != 1)
2278  {
2279  A= mapUp (A, V_buf4, V_buf, prim_elem, im_prim_elem, source,dest);
2280  B= mapUp (B, V_buf4, V_buf, prim_elem, im_prim_elem, source,dest);
2281  }
2282  V_buf4= V_buf;
2283  evalFail= false;
2284  evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf,
2285  evalFail, list);
2286  } while (evalFail);
2287  }
2288  else
2289  {
2291  int deg= 2;
2292  bool initialized= false;
2293  do
2294  {
2295  mipo= randomIrredpoly (deg, x);
2296  if (initialized)
2297  setMipo (V_buf, mipo);
2298  else
2299  V_buf= rootOf (mipo);
2300  evalFail= false;
2301  initialized= true;
2302  evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf,
2303  evalFail, list);
2304  deg++;
2305  } while (evalFail);
2306  V_buf4= V_buf;
2307  }
2308  }
2309 
2310  g= gcd (Aeval, Beval);
2311  g /= Lc (g);
2312 
2313  if (degree (g) != degree (skel) || (skelSize != size (g)))
2314  {
2315  delete[] pEvalPoints;
2316  fail= true;
2317  if (alpha.level() != 1 && V_buf != alpha)
2318  prune1 (alpha);
2319  return 0;
2320  }
2321  CFIterator l= skel;
2322  for (CFIterator k= g; k.hasTerms(); k++, l++)
2323  {
2324  if (k.exp() != l.exp())
2325  {
2326  delete[] pEvalPoints;
2327  fail= true;
2328  if (alpha.level() != 1 && V_buf != alpha)
2329  prune1 (alpha);
2330  return 0;
2331  }
2332  }
2333  pEvalPoints[i]= evalPoints;
2334  gcds[i]= g;
2335 
2336  tmp= 0;
2337  int j= 0;
2338  for (CFListIterator k= evalPoints; k.hasItem(); k++, j++)
2339  tmp += k.getItem()*power (x, j);
2340  list.append (tmp);
2341  }
2342 
2343  if (Monoms.size() == 0)
2344  Monoms= getMonoms (skel);
2345 
2346  coeffMonoms= new CFArray [skelSize];
2347  j= 0;
2348  for (CFIterator i= skel; i.hasTerms(); i++, j++)
2349  coeffMonoms[j]= getMonoms (i.coeff());
2350 
2351  CFArray* pL= new CFArray [skelSize];
2352  CFArray* pM= new CFArray [skelSize];
2353  for (int i= 0; i < biggestSize; i++)
2354  {
2355  CFIterator l= gcds [i];
2356  evalPoints= pEvalPoints [i];
2357  for (int k= 0; k < skelSize; k++, l++)
2358  {
2359  if (i == 0)
2360  pL[k]= CFArray (biggestSize);
2361  pL[k] [i]= l.coeff();
2362 
2363  if (i == 0)
2364  pM[k]= evaluate (coeffMonoms [k], evalPoints);
2365  }
2366  }
2367 
2368  CFArray solution;
2369  CanonicalForm result= 0;
2370  int ind= 0;
2371  CFArray bufArray;
2372  CFMatrix Mat;
2373  for (int k= 0; k < skelSize; k++)
2374  {
2375  if (biggestSize != coeffMonoms[k].size())
2376  {
2377  bufArray= CFArray (coeffMonoms[k].size());
2378  for (int i= 0; i < coeffMonoms[k].size(); i++)
2379  bufArray [i]= pL[k] [i];
2380  solution= solveGeneralVandermonde (pM [k], bufArray);
2381  }
2382  else
2383  solution= solveGeneralVandermonde (pM [k], pL[k]);
2384 
2385  if (solution.size() == 0)
2386  {
2387  delete[] pEvalPoints;
2388  delete[] pM;
2389  delete[] pL;
2390  delete[] coeffMonoms;
2391  fail= true;
2392  if (alpha.level() != 1 && V_buf != alpha)
2393  prune1 (alpha);
2394  return 0;
2395  }
2396  for (int l= 0; l < solution.size(); l++)
2397  result += solution[l]*Monoms [ind + l];
2398  ind += solution.size();
2399  }
2400 
2401  delete[] pEvalPoints;
2402  delete[] pM;
2403  delete[] pL;
2404  delete[] coeffMonoms;
2405 
2406  if (alpha.level() != 1 && V_buf != alpha)
2407  {
2408  CFList u, v;
2409  result= mapDown (result, prim_elem_alpha, im_prim_elem_alpha, alpha, u, v);
2410  prune1 (alpha);
2411  }
2412 
2413  result= N(result);
2414  if (fdivides (result, F) && fdivides (result, G))
2415  return result;
2416  else
2417  {
2418  fail= true;
2419  return 0;
2420  }
2421 }
2422 
2425  const CanonicalForm& skeleton, const Variable& alpha,
2426  bool& fail, CFArray*& coeffMonoms, CFArray& Monoms
2427  )
2428 {
2429  CanonicalForm A= F;
2430  CanonicalForm B= G;
2431  if (F.isZero() && degree(G) > 0) return G/Lc(G);
2432  else if (G.isZero() && degree (F) > 0) return F/Lc(F);
2433  else if (F.isZero() && G.isZero()) return F.genOne();
2434  if (F.inBaseDomain() || G.inBaseDomain()) return F.genOne();
2435  if (F.isUnivariate() && fdivides(F, G)) return F/Lc(F);
2436  if (G.isUnivariate() && fdivides(G, F)) return G/Lc(G);
2437  if (F == G) return F/Lc(F);
2438 
2439  ASSERT (degree (A, 1) == 0, "expected degree (F, 1) == 0");
2440  ASSERT (degree (B, 1) == 0, "expected degree (G, 1) == 0");
2441 
2442  CFMap M,N;
2443  int best_level= myCompress (A, B, M, N, false);
2444 
2445  if (best_level == 0)
2446  return B.genOne();
2447 
2448  A= M(A);
2449  B= M(B);
2450 
2451  Variable x= Variable (1);
2452 
2453  //univariate case
2454  if (A.isUnivariate() && B.isUnivariate())
2455  return N (gcd (A, B));
2456 
2457  CanonicalForm skel= M(skeleton);
2458 
2459  CanonicalForm cA, cB; // content of A and B
2460  CanonicalForm ppA, ppB; // primitive part of A and B
2461  CanonicalForm gcdcAcB;
2462  cA = uni_content (A);
2463  cB = uni_content (B);
2464  gcdcAcB= gcd (cA, cB);
2465  ppA= A/cA;
2466  ppB= B/cB;
2467 
2468  CanonicalForm lcA, lcB; // leading coefficients of A and B
2469  CanonicalForm gcdlcAlcB;
2470  lcA= uni_lcoeff (ppA);
2471  lcB= uni_lcoeff (ppB);
2472 
2473  if (fdivides (lcA, lcB))
2474  {
2475  if (fdivides (A, B))
2476  return F/Lc(F);
2477  }
2478  if (fdivides (lcB, lcA))
2479  {
2480  if (fdivides (B, A))
2481  return G/Lc(G);
2482  }
2483 
2484  gcdlcAlcB= gcd (lcA, lcB);
2485  int skelSize= size (skel, skel.mvar());
2486 
2487  int j= 0;
2488  int biggestSize= 0;
2489  int bufSize;
2490  int numberUni= 0;
2491  for (CFIterator i= skel; i.hasTerms(); i++, j++)
2492  {
2493  bufSize= size (i.coeff());
2494  biggestSize= tmax (biggestSize, bufSize);
2495  numberUni += bufSize;
2496  }
2497  numberUni--;
2498  numberUni= (int) ceil ( (double) numberUni / (skelSize - 1));
2499  biggestSize= tmax (biggestSize , numberUni);
2500 
2501  numberUni= biggestSize + size (LC(skel)) - 1;
2502  int biggestSize2= tmax (numberUni, biggestSize);
2503 
2504  CanonicalForm g, Aeval, Beval;
2505 
2507  CFArray coeffEval;
2508  bool evalFail= false;
2509  CFList list;
2510  bool GF= false;
2511  CanonicalForm LCA= LC (A);
2512  CanonicalForm tmp;
2513  CFArray gcds= CFArray (biggestSize);
2514  CFList * pEvalPoints= new CFList [biggestSize];
2515  Variable V_buf= alpha, V_buf4= alpha;
2516  CFList source, dest;
2517  CanonicalForm prim_elem, im_prim_elem;
2518  CanonicalForm prim_elem_alpha, im_prim_elem_alpha;
2519  for (int i= 0; i < biggestSize; i++)
2520  {
2521  if (i == 0)
2522  {
2523  if (getCharacteristic() > 3)
2524  evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf,
2525  evalFail, list);
2526  else
2527  evalFail= true;
2528 
2529  if (evalFail)
2530  {
2531  if (V_buf.level() != 1)
2532  {
2533  do
2534  {
2535  Variable V_buf3= V_buf;
2536  V_buf= chooseExtension (V_buf);
2537  source= CFList();
2538  dest= CFList();
2539 
2540  bool prim_fail= false;
2541  Variable V_buf2;
2542  prim_elem= primitiveElement (V_buf4, V_buf2, prim_fail);
2543  if (V_buf4 == alpha && alpha.level() != 1)
2544  prim_elem_alpha= prim_elem;
2545 
2546  ASSERT (!prim_fail, "failure in integer factorizer");
2547  if (prim_fail)
2548  ; //ERROR
2549  else
2550  im_prim_elem= mapPrimElem (prim_elem, V_buf4, V_buf);
2551 
2552  DEBOUTLN (cerr, "getMipo (V_buf)= " << getMipo (V_buf));
2553  DEBOUTLN (cerr, "getMipo (V_buf2)= " << getMipo (V_buf2));
2554 
2555  if (V_buf4 == alpha && alpha.level() != 1)
2556  im_prim_elem_alpha= im_prim_elem;
2557  else if (alpha.level() != 1)
2558  im_prim_elem_alpha= mapUp (im_prim_elem_alpha, V_buf4, V_buf,
2559  prim_elem, im_prim_elem, source, dest);
2560 
2561  for (CFListIterator i= list; i.hasItem(); i++)
2562  i.getItem()= mapUp (i.getItem(), V_buf4, V_buf, prim_elem,
2563  im_prim_elem, source, dest);
2564  if (alpha.level() != 1)
2565  {
2566  A= mapUp (A, V_buf4, V_buf, prim_elem, im_prim_elem, source,dest);
2567  B= mapUp (B, V_buf4, V_buf, prim_elem, im_prim_elem, source,dest);
2568  }
2569  evalFail= false;
2570  V_buf4= V_buf;
2571  evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf,
2572  evalFail, list);
2573  } while (evalFail);
2574  }
2575  else
2576  {
2578  int deg= 2;
2579  bool initialized= false;
2580  do
2581  {
2582  mipo= randomIrredpoly (deg, x);
2583  if (initialized)
2584  setMipo (V_buf, mipo);
2585  else
2586  V_buf= rootOf (mipo);
2587  evalFail= false;
2588  initialized= true;
2589  evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf,
2590  evalFail, list);
2591  deg++;
2592  } while (evalFail);
2593  V_buf4= V_buf;
2594  }
2595  }
2596  }
2597  else
2598  {
2599  mult (evalPoints, pEvalPoints[0]);
2600  eval (A, B, Aeval, Beval, evalPoints);
2601  }
2602 
2603  g= gcd (Aeval, Beval);
2604  g /= Lc (g);
2605 
2606  if (degree (g) != degree (skel) || (skelSize != size (g)))
2607  {
2608  delete[] pEvalPoints;
2609  fail= true;
2610  if (alpha.level() != 1 && V_buf != alpha)
2611  prune1 (alpha);
2612  return 0;
2613  }
2614  CFIterator l= skel;
2615  for (CFIterator k= g; k.hasTerms(); k++, l++)
2616  {
2617  if (k.exp() != l.exp())
2618  {
2619  delete[] pEvalPoints;
2620  fail= true;
2621  if (alpha.level() != 1 && V_buf != alpha)
2622  prune1 (alpha);
2623  return 0;
2624  }
2625  }
2626  pEvalPoints[i]= evalPoints;
2627  gcds[i]= g;
2628 
2629  tmp= 0;
2630  int j= 0;
2631  for (CFListIterator k= evalPoints; k.hasItem(); k++, j++)
2632  tmp += k.getItem()*power (x, j);
2633  list.append (tmp);
2634  }
2635 
2636  if (Monoms.size() == 0)
2637  Monoms= getMonoms (skel);
2638 
2639  coeffMonoms= new CFArray [skelSize];
2640 
2641  j= 0;
2642  for (CFIterator i= skel; i.hasTerms(); i++, j++)
2643  coeffMonoms[j]= getMonoms (i.coeff());
2644 
2645  int minimalColumnsIndex;
2646  if (skelSize > 1)
2647  minimalColumnsIndex= 1;
2648  else
2649  minimalColumnsIndex= 0;
2650  int minimalColumns=-1;
2651 
2652  CFArray* pM= new CFArray [skelSize];
2653  CFMatrix Mat;
2654  // find the Matrix with minimal number of columns
2655  for (int i= 0; i < skelSize; i++)
2656  {
2657  pM[i]= CFArray (coeffMonoms[i].size());
2658  if (i == 1)
2659  minimalColumns= coeffMonoms[i].size();
2660  if (i > 1)
2661  {
2662  if (minimalColumns > coeffMonoms[i].size())
2663  {
2664  minimalColumns= coeffMonoms[i].size();
2665  minimalColumnsIndex= i;
2666  }
2667  }
2668  }
2669  CFMatrix* pMat= new CFMatrix [2];
2670  pMat[0]= CFMatrix (biggestSize, coeffMonoms[0].size());
2671  pMat[1]= CFMatrix (biggestSize, coeffMonoms[minimalColumnsIndex].size());
2672  CFArray* pL= new CFArray [skelSize];
2673  for (int i= 0; i < biggestSize; i++)
2674  {
2675  CFIterator l= gcds [i];
2676  evalPoints= pEvalPoints [i];
2677  for (int k= 0; k < skelSize; k++, l++)
2678  {
2679  if (i == 0)
2680  pL[k]= CFArray (biggestSize);
2681  pL[k] [i]= l.coeff();
2682 
2683  if (i == 0 && k != 0 && k != minimalColumnsIndex)
2684  pM[k]= evaluate (coeffMonoms[k], evalPoints);
2685  else if (i == 0 && (k == 0 || k == minimalColumnsIndex))
2686  coeffEval= evaluate (coeffMonoms[k], evalPoints);
2687  if (i > 0 && (k == 0 || k == minimalColumnsIndex))
2688  coeffEval= evaluate (coeffMonoms[k], evalPoints);
2689 
2690  if (k == 0)
2691  {
2692  if (pMat[k].rows() >= i + 1)
2693  {
2694  for (int ind= 1; ind < coeffEval.size() + 1; ind++)
2695  pMat[k] (i + 1, ind)= coeffEval[ind - 1];
2696  }
2697  }
2698  if (k == minimalColumnsIndex)
2699  {
2700  if (pMat[1].rows() >= i + 1)
2701  {
2702  for (int ind= 1; ind < coeffEval.size() + 1; ind++)
2703  pMat[1] (i + 1, ind)= coeffEval[ind - 1];
2704  }
2705  }
2706  }
2707  }
2708 
2709  CFArray solution;
2710  CanonicalForm result= 0;
2711  int ind= 1;
2712  int matRows, matColumns;
2713  matRows= pMat[1].rows();
2714  matColumns= pMat[0].columns() - 1;
2715  matColumns += pMat[1].columns();
2716 
2717  Mat= CFMatrix (matRows, matColumns);
2718  for (int i= 1; i <= matRows; i++)
2719  for (int j= 1; j <= pMat[1].columns(); j++)
2720  Mat (i, j)= pMat[1] (i, j);
2721 
2722  for (int j= pMat[1].columns() + 1; j <= matColumns;
2723  j++, ind++)
2724  {
2725  for (int i= 1; i <= matRows; i++)
2726  Mat (i, j)= (-pMat [0] (i, ind + 1))*pL[minimalColumnsIndex] [i - 1];
2727  }
2728 
2729  CFArray firstColumn= CFArray (pMat[0].rows());
2730  for (int i= 0; i < pMat[0].rows(); i++)
2731  firstColumn [i]= pMat[0] (i + 1, 1);
2732 
2733  CFArray bufArray= pL[minimalColumnsIndex];
2734 
2735  for (int i= 0; i < biggestSize; i++)
2736  pL[minimalColumnsIndex] [i] *= firstColumn[i];
2737 
2738  CFMatrix bufMat= pMat[1];
2739  pMat[1]= Mat;
2740 
2741  if (V_buf.level() != 1)
2742  solution= solveSystemFq (pMat[1],
2743  pL[minimalColumnsIndex], V_buf);
2744  else
2745  solution= solveSystemFp (pMat[1],
2746  pL[minimalColumnsIndex]);
2747 
2748  if (solution.size() == 0)
2749  { //Javadi's method failed, try deKleine, Monagan, Wittkopf instead
2750  CFMatrix bufMat0= pMat[0];
2751  delete [] pMat;
2752  pMat= new CFMatrix [skelSize];
2753  pL[minimalColumnsIndex]= bufArray;
2754  CFList* bufpEvalPoints= NULL;
2755  CFArray bufGcds;
2756  if (biggestSize != biggestSize2)
2757  {
2758  bufpEvalPoints= pEvalPoints;
2759  pEvalPoints= new CFList [biggestSize2];
2760  bufGcds= gcds;
2761  gcds= CFArray (biggestSize2);
2762  for (int i= 0; i < biggestSize; i++)
2763  {
2764  pEvalPoints[i]= bufpEvalPoints [i];
2765  gcds[i]= bufGcds[i];
2766  }
2767  for (int i= 0; i < biggestSize2 - biggestSize; i++)
2768  {
2769  mult (evalPoints, pEvalPoints[0]);
2770  eval (A, B, Aeval, Beval, evalPoints);
2771  g= gcd (Aeval, Beval);
2772  g /= Lc (g);
2773  if (degree (g) != degree (skel) || (skelSize != size (g)))
2774  {
2775  delete[] pEvalPoints;
2776  delete[] pMat;
2777  delete[] pL;
2778  delete[] coeffMonoms;
2779  delete[] pM;
2780  if (bufpEvalPoints != NULL)
2781  delete [] bufpEvalPoints;
2782  fail= true;
2783  if (alpha.level() != 1 && V_buf != alpha)
2784  prune1 (alpha);
2785  return 0;
2786  }
2787  CFIterator l= skel;
2788  for (CFIterator k= g; k.hasTerms(); k++, l++)
2789  {
2790  if (k.exp() != l.exp())
2791  {
2792  delete[] pEvalPoints;
2793  delete[] pMat;
2794  delete[] pL;
2795  delete[] coeffMonoms;
2796  delete[] pM;
2797  if (bufpEvalPoints != NULL)
2798  delete [] bufpEvalPoints;
2799  fail= true;
2800  if (alpha.level() != 1 && V_buf != alpha)
2801  prune1 (alpha);
2802  return 0;
2803  }
2804  }
2805  pEvalPoints[i + biggestSize]= evalPoints;
2806  gcds[i + biggestSize]= g;
2807  }
2808  }
2809  for (int i= 0; i < biggestSize; i++)
2810  {
2811  CFIterator l= gcds [i];
2812  evalPoints= pEvalPoints [i];
2813  for (int k= 1; k < skelSize; k++, l++)
2814  {
2815  if (i == 0)
2816  pMat[k]= CFMatrix (biggestSize2,coeffMonoms[k].size()+biggestSize2-1);
2817  if (k == minimalColumnsIndex)
2818  continue;
2819  coeffEval= evaluate (coeffMonoms[k], evalPoints);
2820  if (pMat[k].rows() >= i + 1)
2821  {
2822  for (int ind= 1; ind < coeffEval.size() + 1; ind++)
2823  pMat[k] (i + 1, ind)= coeffEval[ind - 1];
2824  }
2825  }
2826  }
2827  Mat= bufMat0;
2828  pMat[0]= CFMatrix (biggestSize2, coeffMonoms[0].size() + biggestSize2 - 1);
2829  for (int j= 1; j <= Mat.rows(); j++)
2830  for (int k= 1; k <= Mat.columns(); k++)
2831  pMat [0] (j,k)= Mat (j,k);
2832  Mat= bufMat;
2833  for (int j= 1; j <= Mat.rows(); j++)
2834  for (int k= 1; k <= Mat.columns(); k++)
2835  pMat [minimalColumnsIndex] (j,k)= Mat (j,k);
2836  // write old matrix entries into new matrices
2837  for (int i= 0; i < skelSize; i++)
2838  {
2839  bufArray= pL[i];
2840  pL[i]= CFArray (biggestSize2);
2841  for (int j= 0; j < bufArray.size(); j++)
2842  pL[i] [j]= bufArray [j];
2843  }
2844  //write old vector entries into new and add new entries to old matrices
2845  for (int i= 0; i < biggestSize2 - biggestSize; i++)
2846  {
2847  CFIterator l= gcds [i + biggestSize];
2848  evalPoints= pEvalPoints [i + biggestSize];
2849  for (int k= 0; k < skelSize; k++, l++)
2850  {
2851  pL[k] [i + biggestSize]= l.coeff();
2852  coeffEval= evaluate (coeffMonoms[k], evalPoints);
2853  if (pMat[k].rows() >= i + biggestSize + 1)
2854  {
2855  for (int ind= 1; ind < coeffEval.size() + 1; ind++)
2856  pMat[k] (i + biggestSize + 1, ind)= coeffEval[ind - 1];
2857  }
2858  }
2859  }
2860  // begin new
2861  for (int i= 0; i < skelSize; i++)
2862  {
2863  if (pL[i].size() > 1)
2864  {
2865  for (int j= 2; j <= pMat[i].rows(); j++)
2866  pMat[i] (j, coeffMonoms[i].size() + j - 1)=
2867  -pL[i] [j - 1];
2868  }
2869  }
2870 
2871  matColumns= biggestSize2 - 1;
2872  matRows= 0;
2873  for (int i= 0; i < skelSize; i++)
2874  {
2875  if (V_buf.level() == 1)
2876  (void) gaussianElimFp (pMat[i], pL[i]);
2877  else
2878  (void) gaussianElimFq (pMat[i], pL[i], V_buf);
2879 
2880  if (pMat[i] (coeffMonoms[i].size(), coeffMonoms[i].size()) == 0)
2881  {
2882  delete[] pEvalPoints;
2883  delete[] pMat;
2884  delete[] pL;
2885  delete[] coeffMonoms;
2886  delete[] pM;
2887  if (bufpEvalPoints != NULL)
2888  delete [] bufpEvalPoints;
2889  fail= true;
2890  if (alpha.level() != 1 && V_buf != alpha)
2891  prune1 (alpha);
2892  return 0;
2893  }
2894  matRows += pMat[i].rows() - coeffMonoms[i].size();
2895  }
2896 
2897  CFMatrix bufMat;
2898  Mat= CFMatrix (matRows, matColumns);
2899  ind= 0;
2900  bufArray= CFArray (matRows);
2901  CFArray bufArray2;
2902  for (int i= 0; i < skelSize; i++)
2903  {
2904  if (coeffMonoms[i].size() + 1 >= pMat[i].rows() || coeffMonoms[i].size() + 1 >= pMat[i].columns())
2905  {
2906  delete[] pEvalPoints;
2907  delete[] pMat;
2908  delete[] pL;
2909  delete[] coeffMonoms;
2910  delete[] pM;
2911  if (bufpEvalPoints != NULL)
2912  delete [] bufpEvalPoints;
2913  fail= true;
2914  if (alpha.level() != 1 && V_buf != alpha)
2915  prune1 (alpha);
2916  return 0;
2917  }
2918  bufMat= pMat[i] (coeffMonoms[i].size() + 1, pMat[i].rows(),
2919  coeffMonoms[i].size() + 1, pMat[i].columns());
2920 
2921  for (int j= 1; j <= bufMat.rows(); j++)
2922  for (int k= 1; k <= bufMat.columns(); k++)
2923  Mat (j + ind, k)= bufMat(j, k);
2924  bufArray2= coeffMonoms[i].size();
2925  for (int j= 1; j <= pMat[i].rows(); j++)
2926  {
2927  if (j > coeffMonoms[i].size())
2928  bufArray [j-coeffMonoms[i].size() + ind - 1]= pL[i] [j - 1];
2929  else
2930  bufArray2 [j - 1]= pL[i] [j - 1];
2931  }
2932  pL[i]= bufArray2;
2933  ind += bufMat.rows();
2934  pMat[i]= pMat[i] (1, coeffMonoms[i].size(), 1, pMat[i].columns());
2935  }
2936 
2937  if (V_buf.level() != 1)
2938  solution= solveSystemFq (Mat, bufArray, V_buf);
2939  else
2940  solution= solveSystemFp (Mat, bufArray);
2941 
2942  if (solution.size() == 0)
2943  {
2944  delete[] pEvalPoints;
2945  delete[] pMat;
2946  delete[] pL;
2947  delete[] coeffMonoms;
2948  delete[] pM;
2949  if (bufpEvalPoints != NULL)
2950  delete [] bufpEvalPoints;
2951  fail= true;
2952  if (alpha.level() != 1 && V_buf != alpha)
2953  prune1 (alpha);
2954  return 0;
2955  }
2956 
2957  ind= 0;
2958  result= 0;
2959  CFArray bufSolution;
2960  for (int i= 0; i < skelSize; i++)
2961  {
2962  bufSolution= readOffSolution (pMat[i], pL[i], solution);
2963  for (int i= 0; i < bufSolution.size(); i++, ind++)
2964  result += Monoms [ind]*bufSolution[i];
2965  }
2966  if (alpha.level() != 1 && V_buf != alpha)
2967  {
2968  CFList u, v;
2969  result= mapDown (result,prim_elem_alpha, im_prim_elem_alpha, alpha, u, v);
2970  prune1 (alpha);
2971  }
2972  result= N(result);
2973  delete[] pEvalPoints;
2974  delete[] pMat;
2975  delete[] pL;
2976  delete[] coeffMonoms;
2977  delete[] pM;
2978 
2979  if (bufpEvalPoints != NULL)
2980  delete [] bufpEvalPoints;
2981  if (fdivides (result, F) && fdivides (result, G))
2982  return result;
2983  else
2984  {
2985  fail= true;
2986  return 0;
2987  }
2988  } // end of deKleine, Monagan & Wittkopf
2989 
2990  result += Monoms[0];
2991  int ind2= 0, ind3= 2;
2992  ind= 0;
2993  for (int l= 0; l < minimalColumnsIndex; l++)
2994  ind += coeffMonoms[l].size();
2995  for (int l= solution.size() - pMat[0].columns() + 1; l < solution.size();
2996  l++, ind2++, ind3++)
2997  {
2998  result += solution[l]*Monoms [1 + ind2];
2999  for (int i= 0; i < pMat[0].rows(); i++)
3000  firstColumn[i] += solution[l]*pMat[0] (i + 1, ind3);
3001  }
3002  for (int l= 0; l < solution.size() + 1 - pMat[0].columns(); l++)
3003  result += solution[l]*Monoms [ind + l];
3004  ind= coeffMonoms[0].size();
3005  for (int k= 1; k < skelSize; k++)
3006  {
3007  if (k == minimalColumnsIndex)
3008  {
3009  ind += coeffMonoms[k].size();
3010  continue;
3011  }
3012  if (k != minimalColumnsIndex)
3013  {
3014  for (int i= 0; i < biggestSize; i++)
3015  pL[k] [i] *= firstColumn [i];
3016  }
3017 
3018  if (biggestSize != coeffMonoms[k].size() && k != minimalColumnsIndex)
3019  {
3020  bufArray= CFArray (coeffMonoms[k].size());
3021  for (int i= 0; i < bufArray.size(); i++)
3022  bufArray [i]= pL[k] [i];
3023  solution= solveGeneralVandermonde (pM[k], bufArray);
3024  }
3025  else
3026  solution= solveGeneralVandermonde (pM[k], pL[k]);
3027 
3028  if (solution.size() == 0)
3029  {
3030  delete[] pEvalPoints;
3031  delete[] pMat;
3032  delete[] pL;
3033  delete[] coeffMonoms;
3034  delete[] pM;
3035  fail= true;
3036  if (alpha.level() != 1 && V_buf != alpha)
3037  prune1 (alpha);
3038  return 0;
3039  }
3040  if (k != minimalColumnsIndex)
3041  {
3042  for (int l= 0; l < solution.size(); l++)
3043  result += solution[l]*Monoms [ind + l];
3044  ind += solution.size();
3045  }
3046  }
3047 
3048  delete[] pEvalPoints;
3049  delete[] pMat;
3050  delete[] pL;
3051  delete[] pM;
3052  delete[] coeffMonoms;
3053 
3054  if (alpha.level() != 1 && V_buf != alpha)
3055  {
3056  CFList u, v;
3057  result= mapDown (result, prim_elem, im_prim_elem, alpha, u, v);
3058  prune1 (alpha);
3059  }
3060  result= N(result);
3061 
3062  if (fdivides (result, F) && fdivides (result, G))
3063  return result;
3064  else
3065  {
3066  fail= true;
3067  return 0;
3068  }
3069 }
3070 
3072  const Variable & alpha, CFList& l, bool& topLevel)
3073 {
3074  CanonicalForm A= F;
3075  CanonicalForm B= G;
3076  if (F.isZero() && degree(G) > 0) return G/Lc(G);
3077  else if (G.isZero() && degree (F) > 0) return F/Lc(F);
3078  else if (F.isZero() && G.isZero()) return F.genOne();
3079  if (F.inBaseDomain() || G.inBaseDomain()) return F.genOne();
3080  if (F.isUnivariate() && fdivides(F, G)) return F/Lc(F);
3081  if (G.isUnivariate() && fdivides(G, F)) return G/Lc(G);
3082  if (F == G) return F/Lc(F);
3083 
3084  CFMap M,N;
3085  int best_level= myCompress (A, B, M, N, topLevel);
3086 
3087  if (best_level == 0) return B.genOne();
3088 
3089  A= M(A);
3090  B= M(B);
3091 
3092  Variable x= Variable (1);
3093 
3094  //univariate case
3095  if (A.isUnivariate() && B.isUnivariate())
3096  return N (gcd (A, B));
3097 
3098  CanonicalForm cA, cB; // content of A and B
3099  CanonicalForm ppA, ppB; // primitive part of A and B
3100  CanonicalForm gcdcAcB;
3101 
3102  cA = uni_content (A);
3103  cB = uni_content (B);
3104  gcdcAcB= gcd (cA, cB);
3105  ppA= A/cA;
3106  ppB= B/cB;
3107 
3108  CanonicalForm lcA, lcB; // leading coefficients of A and B
3109  CanonicalForm gcdlcAlcB;
3110  lcA= uni_lcoeff (ppA);
3111  lcB= uni_lcoeff (ppB);
3112 
3113  if (fdivides (lcA, lcB))
3114  {
3115  if (fdivides (A, B))
3116  return F/Lc(F);
3117  }
3118  if (fdivides (lcB, lcA))
3119  {
3120  if (fdivides (B, A))
3121  return G/Lc(G);
3122  }
3123 
3124  gcdlcAlcB= gcd (lcA, lcB);
3125 
3126  int d= totaldegree (ppA, Variable (2), Variable (ppA.level()));
3127  int d0;
3128 
3129  if (d == 0)
3130  return N(gcdcAcB);
3131  d0= totaldegree (ppB, Variable (2), Variable (ppB.level()));
3132 
3133  if (d0 < d)
3134  d= d0;
3135 
3136  if (d == 0)
3137  return N(gcdcAcB);
3138 
3139  CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH, skeleton;
3140  CanonicalForm newtonPoly= 1;
3141  m= gcdlcAlcB;
3142  G_m= 0;
3143  H= 0;
3144  bool fail= false;
3145  topLevel= false;
3146  bool inextension= false;
3147  Variable V_buf= alpha, V_buf4= alpha;
3148  CanonicalForm prim_elem, im_prim_elem;
3149  CanonicalForm prim_elem_alpha, im_prim_elem_alpha;
3150  CFList source, dest;
3151  do // first do
3152  {
3153  random_element= randomElement (m, V_buf, l, fail);
3154  if (random_element == 0 && !fail)
3155  {
3156  if (!find (l, random_element))
3157  l.append (random_element);
3158  continue;
3159  }
3160  if (fail)
3161  {
3162  source= CFList();
3163  dest= CFList();
3164 
3165  Variable V_buf3= V_buf;
3166  V_buf= chooseExtension (V_buf);
3167  bool prim_fail= false;
3168  Variable V_buf2;
3169  prim_elem= primitiveElement (V_buf4, V_buf2, prim_fail);
3170  if (V_buf4 == alpha)
3171  prim_elem_alpha= prim_elem;
3172 
3173  if (V_buf3 != V_buf4)
3174  {
3175  m= mapDown (m, prim_elem, im_prim_elem, V_buf4, source, dest);
3176  G_m= mapDown (m, prim_elem, im_prim_elem, V_buf4, source, dest);
3177  newtonPoly= mapDown (newtonPoly, prim_elem, im_prim_elem, V_buf4,
3178  source, dest);
3179  ppA= mapDown (ppA, prim_elem, im_prim_elem, V_buf4, source, dest);
3180  ppB= mapDown (ppB, prim_elem, im_prim_elem, V_buf4, source, dest);
3181  gcdlcAlcB= mapDown (gcdlcAlcB, prim_elem, im_prim_elem, V_buf4, source,
3182  dest);
3183  for (CFListIterator i= l; i.hasItem(); i++)
3184  i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, V_buf4,
3185  source, dest);
3186  }
3187 
3188  ASSERT (!prim_fail, "failure in integer factorizer");
3189  if (prim_fail)
3190  ; //ERROR
3191  else
3192  im_prim_elem= mapPrimElem (prim_elem, V_buf4, V_buf);
3193 
3194  if (V_buf4 == alpha)
3195  im_prim_elem_alpha= im_prim_elem;
3196  else
3197  im_prim_elem_alpha= mapUp (im_prim_elem_alpha, V_buf4, V_buf, prim_elem,
3198  im_prim_elem, source, dest);
3199 
3200  DEBOUTLN (cerr, "getMipo (V_buf4)= " << getMipo (V_buf4));
3201  DEBOUTLN (cerr, "getMipo (V_buf2)= " << getMipo (V_buf2));
3202  inextension= true;
3203  for (CFListIterator i= l; i.hasItem(); i++)
3204  i.getItem()= mapUp (i.getItem(), V_buf4, V_buf, prim_elem,
3205  im_prim_elem, source, dest);
3206  m= mapUp (m, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
3207  G_m= mapUp (G_m, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
3208  newtonPoly= mapUp (newtonPoly, V_buf4, V_buf, prim_elem, im_prim_elem,
3209  source, dest);
3210  ppA= mapUp (ppA, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
3211  ppB= mapUp (ppB, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
3212  gcdlcAlcB= mapUp (gcdlcAlcB, V_buf4, V_buf, prim_elem, im_prim_elem,
3213  source, dest);
3214 
3215  fail= false;
3216  random_element= randomElement (m, V_buf, l, fail );
3217  DEBOUTLN (cerr, "fail= " << fail);
3218  CFList list;
3219  TIMING_START (gcd_recursion);
3220  G_random_element=
3221  sparseGCDFq (ppA (random_element, x), ppB (random_element, x), V_buf,
3222  list, topLevel);
3223  TIMING_END_AND_PRINT (gcd_recursion,
3224  "time for recursive call: ");
3225  DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3226  V_buf4= V_buf;
3227  }
3228  else
3229  {
3230  CFList list;
3231  TIMING_START (gcd_recursion);
3232  G_random_element=
3233  sparseGCDFq (ppA(random_element, x), ppB(random_element, x), V_buf,
3234  list, topLevel);
3235  TIMING_END_AND_PRINT (gcd_recursion,
3236  "time for recursive call: ");
3237  DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3238  }
3239 
3240  if (!G_random_element.inCoeffDomain())
3241  d0= totaldegree (G_random_element, Variable(2),
3242  Variable (G_random_element.level()));
3243  else
3244  d0= 0;
3245 
3246  if (d0 == 0)
3247  {
3248  if (inextension)
3249  prune1 (alpha);
3250  return N(gcdcAcB);
3251  }
3252  if (d0 > d)
3253  {
3254  if (!find (l, random_element))
3255  l.append (random_element);
3256  continue;
3257  }
3258 
3259  G_random_element=
3260  (gcdlcAlcB(random_element, x)/uni_lcoeff (G_random_element))
3261  * G_random_element;
3262 
3263  skeleton= G_random_element;
3264  if (!G_random_element.inCoeffDomain())
3265  d0= totaldegree (G_random_element, Variable(2),
3266  Variable (G_random_element.level()));
3267  else
3268  d0= 0;
3269 
3270  if (d0 < d)
3271  {
3272  m= gcdlcAlcB;
3273  newtonPoly= 1;
3274  G_m= 0;
3275  d= d0;
3276  }
3277 
3278  H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
3279  if (uni_lcoeff (H) == gcdlcAlcB)
3280  {
3281  cH= uni_content (H);
3282  ppH= H/cH;
3283  if (inextension)
3284  {
3285  CFList u, v;
3286  //maybe it's better to test if ppH is an element of F(\alpha) before
3287  //mapping down
3288  if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
3289  {
3290  DEBOUTLN (cerr, "ppH before mapDown= " << ppH);
3291  ppH= mapDown (ppH, prim_elem_alpha, im_prim_elem_alpha, alpha, u, v);
3292  ppH /= Lc(ppH);
3293  DEBOUTLN (cerr, "ppH after mapDown= " << ppH);
3294  prune1 (alpha);
3295  return N(gcdcAcB*ppH);
3296  }
3297  }
3298  else if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
3299  return N(gcdcAcB*ppH);
3300  }
3301  G_m= H;
3302  newtonPoly= newtonPoly*(x - random_element);
3303  m= m*(x - random_element);
3304  if (!find (l, random_element))
3305  l.append (random_element);
3306 
3307  if (getCharacteristic () > 3 && size (skeleton) < 100)
3308  {
3309  CFArray Monoms;
3310  CFArray *coeffMonoms;
3311  do //second do
3312  {
3313  random_element= randomElement (m, V_buf, l, fail);
3314  if (random_element == 0 && !fail)
3315  {
3316  if (!find (l, random_element))
3317  l.append (random_element);
3318  continue;
3319  }
3320  if (fail)
3321  {
3322  source= CFList();
3323  dest= CFList();
3324 
3325  Variable V_buf3= V_buf;
3326  V_buf= chooseExtension (V_buf);
3327  bool prim_fail= false;
3328  Variable V_buf2;
3329  prim_elem= primitiveElement (V_buf4, V_buf2, prim_fail);
3330  if (V_buf4 == alpha)
3331  prim_elem_alpha= prim_elem;
3332 
3333  if (V_buf3 != V_buf4)
3334  {
3335  m= mapDown (m, prim_elem, im_prim_elem, V_buf4, source, dest);
3336  G_m= mapDown (m, prim_elem, im_prim_elem, V_buf4, source, dest);
3337  newtonPoly= mapDown (newtonPoly, prim_elem, im_prim_elem, V_buf4,
3338  source, dest);
3339  ppA= mapDown (ppA, prim_elem, im_prim_elem, V_buf4, source, dest);
3340  ppB= mapDown (ppB, prim_elem, im_prim_elem, V_buf4, source, dest);
3341  gcdlcAlcB= mapDown (gcdlcAlcB, prim_elem, im_prim_elem, V_buf4,
3342  source, dest);
3343  for (CFListIterator i= l; i.hasItem(); i++)
3344  i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, V_buf4,
3345  source, dest);
3346  }
3347 
3348  ASSERT (!prim_fail, "failure in integer factorizer");
3349  if (prim_fail)
3350  ; //ERROR
3351  else
3352  im_prim_elem= mapPrimElem (prim_elem, V_buf4, V_buf);
3353 
3354  if (V_buf4 == alpha)
3355  im_prim_elem_alpha= im_prim_elem;
3356  else
3357  im_prim_elem_alpha= mapUp (im_prim_elem_alpha, V_buf4, V_buf,
3358  prim_elem, im_prim_elem, source, dest);
3359 
3360  DEBOUTLN (cerr, "getMipo (V_buf4)= " << getMipo (V_buf4));
3361  DEBOUTLN (cerr, "getMipo (V_buf2)= " << getMipo (V_buf2));
3362  inextension= true;
3363  for (CFListIterator i= l; i.hasItem(); i++)
3364  i.getItem()= mapUp (i.getItem(), V_buf4, V_buf, prim_elem,
3365  im_prim_elem, source, dest);
3366  m= mapUp (m, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
3367  G_m= mapUp (G_m, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
3368  newtonPoly= mapUp (newtonPoly, V_buf4, V_buf, prim_elem, im_prim_elem,
3369  source, dest);
3370  ppA= mapUp (ppA, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
3371  ppB= mapUp (ppB, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
3372 
3373  gcdlcAlcB= mapUp (gcdlcAlcB, V_buf4, V_buf, prim_elem, im_prim_elem,
3374  source, dest);
3375 
3376  fail= false;
3377  random_element= randomElement (m, V_buf, l, fail);
3378  DEBOUTLN (cerr, "fail= " << fail);
3379  CFList list;
3380  TIMING_START (gcd_recursion);
3381 
3382  V_buf4= V_buf;
3383 
3384  //sparseInterpolation
3385  bool sparseFail= false;
3386  if (LC (skeleton).inCoeffDomain())
3387  G_random_element=
3388  monicSparseInterpol (ppA (random_element, x), ppB(random_element,x),
3389  skeleton,V_buf, sparseFail, coeffMonoms,Monoms);
3390  else
3391  G_random_element=
3392  nonMonicSparseInterpol (ppA(random_element,x),ppB(random_element,x),
3393  skeleton, V_buf, sparseFail, coeffMonoms,
3394  Monoms);
3395  if (sparseFail)
3396  break;
3397 
3398  TIMING_END_AND_PRINT (gcd_recursion,
3399  "time for recursive call: ");
3400  DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3401  }
3402  else
3403  {
3404  CFList list;
3405  TIMING_START (gcd_recursion);
3406  bool sparseFail= false;
3407  if (LC (skeleton).inCoeffDomain())
3408  G_random_element=
3409  monicSparseInterpol (ppA (random_element, x),ppB(random_element, x),
3410  skeleton,V_buf, sparseFail,coeffMonoms, Monoms);
3411  else
3412  G_random_element=
3413  nonMonicSparseInterpol (ppA(random_element,x),ppB(random_element,x),
3414  skeleton, V_buf, sparseFail, coeffMonoms,
3415  Monoms);
3416  if (sparseFail)
3417  break;
3418 
3419  TIMING_END_AND_PRINT (gcd_recursion,
3420  "time for recursive call: ");
3421  DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3422  }
3423 
3424  if (!G_random_element.inCoeffDomain())
3425  d0= totaldegree (G_random_element, Variable(2),
3426  Variable (G_random_element.level()));
3427  else
3428  d0= 0;
3429 
3430  if (d0 == 0)
3431  {
3432  if (inextension)
3433  prune1 (alpha);
3434  return N(gcdcAcB);
3435  }
3436  if (d0 > d)
3437  {
3438  if (!find (l, random_element))
3439  l.append (random_element);
3440  continue;
3441  }
3442 
3443  G_random_element=
3444  (gcdlcAlcB(random_element, x)/uni_lcoeff (G_random_element))
3445  * G_random_element;
3446 
3447  if (!G_random_element.inCoeffDomain())
3448  d0= totaldegree (G_random_element, Variable(2),
3449  Variable (G_random_element.level()));
3450  else
3451  d0= 0;
3452 
3453  if (d0 < d)
3454  {
3455  m= gcdlcAlcB;
3456  newtonPoly= 1;
3457  G_m= 0;
3458  d= d0;
3459  }
3460 
3461  TIMING_START (newton_interpolation);
3462  H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
3463  TIMING_END_AND_PRINT (newton_interpolation,
3464  "time for newton interpolation: ");
3465 
3466  //termination test
3467  if (uni_lcoeff (H) == gcdlcAlcB)
3468  {
3469  cH= uni_content (H);
3470  ppH= H/cH;
3471  if (inextension)
3472  {
3473  CFList u, v;
3474  //maybe it's better to test if ppH is an element of F(\alpha) before
3475  //mapping down
3476  if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
3477  {
3478  DEBOUTLN (cerr, "ppH before mapDown= " << ppH);
3479  ppH= mapDown (ppH, prim_elem_alpha, im_prim_elem_alpha, alpha, u, v);
3480  ppH /= Lc(ppH);
3481  DEBOUTLN (cerr, "ppH after mapDown= " << ppH);
3482  prune1 (alpha);
3483  return N(gcdcAcB*ppH);
3484  }
3485  }
3486  else if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
3487  {
3488  return N(gcdcAcB*ppH);
3489  }
3490  }
3491 
3492  G_m= H;
3493  newtonPoly= newtonPoly*(x - random_element);
3494  m= m*(x - random_element);
3495  if (!find (l, random_element))
3496  l.append (random_element);
3497 
3498  } while (1);
3499  }
3500  } while (1);
3501 }
3502 
3504  bool& topLevel, CFList& l)
3505 {
3506  CanonicalForm A= F;
3507  CanonicalForm B= G;
3508  if (F.isZero() && degree(G) > 0) return G/Lc(G);
3509  else if (G.isZero() && degree (F) > 0) return F/Lc(F);
3510  else if (F.isZero() && G.isZero()) return F.genOne();
3511  if (F.inBaseDomain() || G.inBaseDomain()) return F.genOne();
3512  if (F.isUnivariate() && fdivides(F, G)) return F/Lc(F);
3513  if (G.isUnivariate() && fdivides(G, F)) return G/Lc(G);
3514  if (F == G) return F/Lc(F);
3515 
3516  CFMap M,N;
3517  int best_level= myCompress (A, B, M, N, topLevel);
3518 
3519  if (best_level == 0) return B.genOne();
3520 
3521  A= M(A);
3522  B= M(B);
3523 
3524  Variable x= Variable (1);
3525 
3526  //univariate case
3527  if (A.isUnivariate() && B.isUnivariate())
3528  return N (gcd (A, B));
3529 
3530  CanonicalForm cA, cB; // content of A and B
3531  CanonicalForm ppA, ppB; // primitive part of A and B
3532  CanonicalForm gcdcAcB;
3533 
3534  cA = uni_content (A);
3535  cB = uni_content (B);
3536  gcdcAcB= gcd (cA, cB);
3537  ppA= A/cA;
3538  ppB= B/cB;
3539 
3540  CanonicalForm lcA, lcB; // leading coefficients of A and B
3541  CanonicalForm gcdlcAlcB;
3542  lcA= uni_lcoeff (ppA);
3543  lcB= uni_lcoeff (ppB);
3544 
3545  if (fdivides (lcA, lcB))
3546  {
3547  if (fdivides (A, B))
3548  return F/Lc(F);
3549  }
3550  if (fdivides (lcB, lcA))
3551  {
3552  if (fdivides (B, A))
3553  return G/Lc(G);
3554  }
3555 
3556  gcdlcAlcB= gcd (lcA, lcB);
3557 
3558  int d= totaldegree (ppA, Variable (2), Variable (ppA.level()));
3559  int d0;
3560 
3561  if (d == 0)
3562  return N(gcdcAcB);
3563 
3564  d0= totaldegree (ppB, Variable (2), Variable (ppB.level()));
3565 
3566  if (d0 < d)
3567  d= d0;
3568 
3569  if (d == 0)
3570  return N(gcdcAcB);
3571 
3572  CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH, skeleton;
3573  CanonicalForm newtonPoly= 1;
3574  m= gcdlcAlcB;
3575  G_m= 0;
3576  H= 0;
3577  bool fail= false;
3578  topLevel= false;
3579  bool inextension= false;
3580  Variable V_buf, alpha, cleanUp;
3581  CanonicalForm prim_elem, im_prim_elem;
3582  CFList source, dest;
3583  do //first do
3584  {
3585  if (inextension)
3586  random_element= randomElement (m, V_buf, l, fail);
3587  else
3588  random_element= FpRandomElement (m, l, fail);
3589  if (random_element == 0 && !fail)
3590  {
3591  if (!find (l, random_element))
3592  l.append (random_element);
3593  continue;
3594  }
3595 
3596  if (!fail && !inextension)
3597  {
3598  CFList list;
3599  TIMING_START (gcd_recursion);
3600  G_random_element=
3601  sparseGCDFp (ppA (random_element,x), ppB (random_element,x), topLevel,
3602  list);
3603  TIMING_END_AND_PRINT (gcd_recursion,
3604  "time for recursive call: ");
3605  DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3606  }
3607  else if (!fail && inextension)
3608  {
3609  CFList list;
3610  TIMING_START (gcd_recursion);
3611  G_random_element=
3612  sparseGCDFq (ppA (random_element, x), ppB (random_element, x), alpha,
3613  list, topLevel);
3614  TIMING_END_AND_PRINT (gcd_recursion,
3615  "time for recursive call: ");
3616  DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3617  }
3618  else if (fail && !inextension)
3619  {
3620  source= CFList();
3621  dest= CFList();
3622  CFList list;
3624  int deg= 2;
3625  bool initialized= false;
3626  do
3627  {
3628  mipo= randomIrredpoly (deg, x);
3629  if (initialized)
3630  setMipo (alpha, mipo);
3631  else
3632  alpha= rootOf (mipo);
3633  inextension= true;
3634  fail= false;
3635  initialized= true;
3636  random_element= randomElement (m, alpha, l, fail);
3637  deg++;
3638  } while (fail);
3639  cleanUp= alpha;
3640  V_buf= alpha;
3641  list= CFList();
3642  TIMING_START (gcd_recursion);
3643  G_random_element=
3644  sparseGCDFq (ppA (random_element, x), ppB (random_element, x), alpha,
3645  list, topLevel);
3646  TIMING_END_AND_PRINT (gcd_recursion,
3647  "time for recursive call: ");
3648  DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3649  }
3650  else if (fail && inextension)
3651  {
3652  source= CFList();
3653  dest= CFList();
3654 
3655  Variable V_buf3= V_buf;
3656  V_buf= chooseExtension (V_buf);
3657  bool prim_fail= false;
3658  Variable V_buf2;
3659  CanonicalForm prim_elem, im_prim_elem;
3660  prim_elem= primitiveElement (alpha, V_buf2, prim_fail);
3661 
3662  if (V_buf3 != alpha)
3663  {
3664  m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
3665  G_m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
3666  newtonPoly= mapDown (newtonPoly, prim_elem, im_prim_elem, alpha, source,
3667  dest);
3668  ppA= mapDown (ppA, prim_elem, im_prim_elem, alpha, source, dest);
3669  ppB= mapDown (ppB, prim_elem, im_prim_elem, alpha, source, dest);
3670  gcdlcAlcB= mapDown (gcdlcAlcB, prim_elem, im_prim_elem, alpha, source,
3671  dest);
3672  for (CFListIterator i= l; i.hasItem(); i++)
3673  i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, alpha,
3674  source, dest);
3675  }
3676 
3677  ASSERT (!prim_fail, "failure in integer factorizer");
3678  if (prim_fail)
3679  ; //ERROR
3680  else
3681  im_prim_elem= mapPrimElem (prim_elem, alpha, V_buf);
3682 
3683  DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (alpha));
3684  DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf2));
3685 
3686  for (CFListIterator i= l; i.hasItem(); i++)
3687  i.getItem()= mapUp (i.getItem(), alpha, V_buf, prim_elem,
3688  im_prim_elem, source, dest);
3689  m= mapUp (m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3690  G_m= mapUp (G_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3691  newtonPoly= mapUp (newtonPoly, alpha, V_buf, prim_elem, im_prim_elem,
3692  source, dest);
3693  ppA= mapUp (ppA, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3694  ppB= mapUp (ppB, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3695  gcdlcAlcB= mapUp (gcdlcAlcB, alpha, V_buf, prim_elem, im_prim_elem,
3696  source, dest);
3697  fail= false;
3698  random_element= randomElement (m, V_buf, l, fail );
3699  DEBOUTLN (cerr, "fail= " << fail);
3700  CFList list;
3701  TIMING_START (gcd_recursion);
3702  G_random_element=
3703  sparseGCDFq (ppA (random_element, x), ppB (random_element, x), V_buf,
3704  list, topLevel);
3705  TIMING_END_AND_PRINT (gcd_recursion,
3706  "time for recursive call: ");
3707  DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3708  alpha= V_buf;
3709  }
3710 
3711  if (!G_random_element.inCoeffDomain())
3712  d0= totaldegree (G_random_element, Variable(2),
3713  Variable (G_random_element.level()));
3714  else
3715  d0= 0;
3716 
3717  if (d0 == 0)
3718  {
3719  if (inextension)
3720  prune (cleanUp);
3721  return N(gcdcAcB);
3722  }
3723  if (d0 > d)
3724  {
3725  if (!find (l, random_element))
3726  l.append (random_element);
3727  continue;
3728  }
3729 
3730  G_random_element=
3731  (gcdlcAlcB(random_element, x)/uni_lcoeff (G_random_element))
3732  * G_random_element;
3733 
3734  skeleton= G_random_element;
3735 
3736  if (!G_random_element.inCoeffDomain())
3737  d0= totaldegree (G_random_element, Variable(2),
3738  Variable (G_random_element.level()));
3739  else
3740  d0= 0;
3741 
3742  if (d0 < d)
3743  {
3744  m= gcdlcAlcB;
3745  newtonPoly= 1;
3746  G_m= 0;
3747  d= d0;
3748  }
3749 
3750  H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
3751 
3752  if (uni_lcoeff (H) == gcdlcAlcB)
3753  {
3754  cH= uni_content (H);
3755  ppH= H/cH;
3756  ppH /= Lc (ppH);
3757  DEBOUTLN (cerr, "ppH= " << ppH);
3758 
3759  if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
3760  {
3761  if (inextension)
3762  prune (cleanUp);
3763  return N(gcdcAcB*ppH);
3764  }
3765  }
3766  G_m= H;
3767  newtonPoly= newtonPoly*(x - random_element);
3768  m= m*(x - random_element);
3769  if (!find (l, random_element))
3770  l.append (random_element);
3771 
3772  if ((getCharacteristic() > 3 && size (skeleton) < 200))
3773  {
3774  CFArray Monoms;
3775  CFArray* coeffMonoms;
3776 
3777  do //second do
3778  {
3779  if (inextension)
3780  random_element= randomElement (m, alpha, l, fail);
3781  else
3782  random_element= FpRandomElement (m, l, fail);
3783  if (random_element == 0 && !fail)
3784  {
3785  if (!find (l, random_element))
3786  l.append (random_element);
3787  continue;
3788  }
3789 
3790  bool sparseFail= false;
3791  if (!fail && !inextension)
3792  {
3793  CFList list;
3794  TIMING_START (gcd_recursion);
3795  if (LC (skeleton).inCoeffDomain())
3796  G_random_element=
3797  monicSparseInterpol(ppA(random_element, x), ppB (random_element, x),
3798  skeleton, x, sparseFail, coeffMonoms,
3799  Monoms);
3800  else
3801  G_random_element=
3802  nonMonicSparseInterpol(ppA(random_element,x), ppB(random_element,x),
3803  skeleton, x, sparseFail,
3804  coeffMonoms, Monoms);
3805  TIMING_END_AND_PRINT (gcd_recursion,
3806  "time for recursive call: ");
3807  DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3808  }
3809  else if (!fail && inextension)
3810  {
3811  CFList list;
3812  TIMING_START (gcd_recursion);
3813  if (LC (skeleton).inCoeffDomain())
3814  G_random_element=
3815  monicSparseInterpol(ppA (random_element,x), ppB (random_element, x),
3816  skeleton, alpha, sparseFail, coeffMonoms,
3817  Monoms);
3818  else
3819  G_random_element=
3820  nonMonicSparseInterpol(ppA(random_element,x), ppB(random_element,x),
3821  skeleton, alpha, sparseFail, coeffMonoms,
3822  Monoms);
3823  TIMING_END_AND_PRINT (gcd_recursion,
3824  "time for recursive call: ");
3825  DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3826  }
3827  else if (fail && !inextension)
3828  {
3829  source= CFList();
3830  dest= CFList();
3831  CFList list;
3833  int deg= 2;
3834  bool initialized= false;
3835  do
3836  {
3837  mipo= randomIrredpoly (deg, x);
3838  if (initialized)
3839  setMipo (alpha, mipo);
3840  else
3841  alpha= rootOf (mipo);
3842  inextension= true;
3843  fail= false;
3844  initialized= true;
3845  random_element= randomElement (m, alpha, l, fail);
3846  deg++;
3847  } while (fail);
3848  cleanUp= alpha;
3849  V_buf= alpha;
3850  list= CFList();
3851  TIMING_START (gcd_recursion);
3852  if (LC (skeleton).inCoeffDomain())
3853  G_random_element=
3854  monicSparseInterpol (ppA (random_element,x), ppB (random_element,x),
3855  skeleton, alpha, sparseFail, coeffMonoms,
3856  Monoms);
3857  else
3858  G_random_element=
3859  nonMonicSparseInterpol(ppA(random_element,x), ppB(random_element,x),
3860  skeleton, alpha, sparseFail, coeffMonoms,
3861  Monoms);
3862  TIMING_END_AND_PRINT (gcd_recursion,
3863  "time for recursive call: ");
3864  DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3865  }
3866  else if (fail && inextension)
3867  {
3868  source= CFList();
3869  dest= CFList();
3870 
3871  Variable V_buf3= V_buf;
3872  V_buf= chooseExtension (V_buf);
3873  bool prim_fail= false;
3874  Variable V_buf2;
3875  CanonicalForm prim_elem, im_prim_elem;
3876  prim_elem= primitiveElement (alpha, V_buf2, prim_fail);
3877 
3878  if (V_buf3 != alpha)
3879  {
3880  m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
3881  G_m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
3882  newtonPoly= mapDown (newtonPoly, prim_elem, im_prim_elem, alpha,
3883  source, dest);
3884  ppA= mapDown (ppA, prim_elem, im_prim_elem, alpha, source, dest);
3885  ppB= mapDown (ppB, prim_elem, im_prim_elem, alpha, source, dest);
3886  gcdlcAlcB= mapDown (gcdlcAlcB, prim_elem, im_prim_elem, alpha,
3887  source, dest);
3888  for (CFListIterator i= l; i.hasItem(); i++)
3889  i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, alpha,
3890  source, dest);
3891  }
3892 
3893  ASSERT (!prim_fail, "failure in integer factorizer");
3894  if (prim_fail)
3895  ; //ERROR
3896  else
3897  im_prim_elem= mapPrimElem (prim_elem, alpha, V_buf);
3898 
3899  DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (alpha));
3900  DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf2));
3901 
3902  for (CFListIterator i= l; i.hasItem(); i++)
3903  i.getItem()= mapUp (i.getItem(), alpha, V_buf, prim_elem,
3904  im_prim_elem, source, dest);
3905  m= mapUp (m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3906  G_m= mapUp (G_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3907  newtonPoly= mapUp (newtonPoly, alpha, V_buf, prim_elem, im_prim_elem,
3908  source, dest);
3909  ppA= mapUp (ppA, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3910  ppB= mapUp (ppB, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3911  gcdlcAlcB= mapUp (gcdlcAlcB, alpha, V_buf, prim_elem, im_prim_elem,
3912  source, dest);
3913  fail= false;
3914  random_element= randomElement (m, V_buf, l, fail );
3915  DEBOUTLN (cerr, "fail= " << fail);
3916  CFList list;
3917  TIMING_START (gcd_recursion);
3918  if (LC (skeleton).inCoeffDomain())
3919  G_random_element=
3920  monicSparseInterpol (ppA (random_element, x), ppB (random_element, x),
3921  skeleton, V_buf, sparseFail, coeffMonoms,
3922  Monoms);
3923  else
3924  G_random_element=
3925  nonMonicSparseInterpol (ppA(random_element,x), ppB(random_element,x),
3926  skeleton, V_buf, sparseFail,
3927  coeffMonoms, Monoms);
3928  TIMING_END_AND_PRINT (gcd_recursion,
3929  "time for recursive call: ");
3930  DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3931  alpha= V_buf;
3932  }
3933 
3934  if (sparseFail)
3935  break;
3936 
3937  if (!G_random_element.inCoeffDomain())
3938  d0= totaldegree (G_random_element, Variable(2),
3939  Variable (G_random_element.level()));
3940  else
3941  d0= 0;
3942 
3943  if (d0 == 0)
3944  {
3945  if (inextension)
3946  prune (cleanUp);
3947  return N(gcdcAcB);
3948  }
3949  if (d0 > d)
3950  {
3951  if (!find (l, random_element))
3952  l.append (random_element);
3953  continue;
3954  }
3955 
3956  G_random_element=
3957  (gcdlcAlcB(random_element, x)/uni_lcoeff (G_random_element))
3958  * G_random_element;
3959 
3960  if (!G_random_element.inCoeffDomain())
3961  d0= totaldegree (G_random_element, Variable(2),
3962  Variable (G_random_element.level()));
3963  else
3964  d0= 0;
3965 
3966  if (d0 < d)
3967  {
3968  m= gcdlcAlcB;
3969  newtonPoly= 1;
3970  G_m= 0;
3971  d= d0;
3972  }
3973 
3974  TIMING_START (newton_interpolation);
3975  H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
3976  TIMING_END_AND_PRINT (newton_interpolation,
3977  "time for newton interpolation: ");
3978 
3979  //termination test
3980  if (uni_lcoeff (H) == gcdlcAlcB)
3981  {
3982  cH= uni_content (H);
3983  ppH= H/cH;
3984  ppH /= Lc (ppH);
3985  DEBOUTLN (cerr, "ppH= " << ppH);
3986  if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
3987  {
3988  if (inextension)
3989  prune (cleanUp);
3990  return N(gcdcAcB*ppH);
3991  }
3992  }
3993 
3994  G_m= H;
3995  newtonPoly= newtonPoly*(x - random_element);
3996  m= m*(x - random_element);
3997  if (!find (l, random_element))
3998  l.append (random_element);
3999 
4000  } while (1); //end of second do
4001  }
4002  else
4003  {
4004  if (inextension)
4005  prune (cleanUp);
4006  return N(gcdcAcB*modGCDFp (ppA, ppB));
4007  }
4008  } while (1); //end of first do
4009 }
4010 
4011 TIMING_DEFINE_PRINT(modZ_termination)
4012 TIMING_DEFINE_PRINT(modZ_recursion)
4013 
4014 /// modular gcd algorithm, see Keith, Czapor, Geddes "Algorithms for Computer
4015 /// Algebra", Algorithm 7.1
4016 CanonicalForm modGCDZ ( const CanonicalForm & FF, const CanonicalForm & GG )
4018  CanonicalForm f, g, cl, q(0), Dp, newD, D, newq, newqh;
4019  int p, i, dp_deg, d_deg=-1;
4020 
4021  CanonicalForm cd ( bCommonDen( FF ));
4022  f=cd*FF;
4025  cf= icontent (f);
4026  f /= cf;
4027  //cd = bCommonDen( f ); f *=cd;
4028  //f /=vcontent(f,Variable(1));
4029 
4030  cd = bCommonDen( GG );
4031  g=cd*GG;
4032  cg= icontent (g);
4033  g /= cg;
4034  //cd = bCommonDen( g ); g *=cd;
4035  //g /=vcontent(g,Variable(1));
4036 
4039  lcf= f.lc();
4040  lcg= g.lc();
4041  cl = gcd (f.lc(),g.lc());
4046  for (i= tmax (f.level(), g.level()); i > 0; i--)
4047  {
4048  if (degree (f, i) <= 0 || degree (g, i) <= 0)
4049  continue;
4050  else
4051  {
4052  minCommonDeg= tmin (degree (g, i), degree (f, i));
4053  break;
4054  }
4055  }
4056  if (i == 0)
4057  return gcdcfcg;
4058  for (; i > 0; i--)
4059  {
4060  if (degree (f, i) <= 0 || degree (g, i) <= 0)
4061  continue;
4062  else
4063  minCommonDeg= tmin (minCommonDeg, tmin (degree (g, i), degree (f, i)));
4064  }
4065  b= 2*tmin (maxNorm (f), maxNorm (g))*abs (cl)*
4066  power (CanonicalForm (2), minCommonDeg);
4067  bool equal= false;
4068  i = cf_getNumBigPrimes() - 1;
4069 
4072  //Off (SW_RATIONAL);
4073  while ( true )
4074  {
4075  p = cf_getBigPrime( i );
4076  i--;
4077  while ( i >= 0 && mod( cl*(lc(f)/cl)*(lc(g)/cl), p ) == 0 )
4078  {
4079  p = cf_getBigPrime( i );
4080  i--;
4081  }
4082  //printf("try p=%d\n",p);
4083  setCharacteristic( p );
4084  fp= mapinto (f);
4085  gp= mapinto (g);
4086  TIMING_START (modZ_recursion)
4087 #ifdef HAVE_NTL
4088  if (size (fp)/maxNumVars > 500 && size (gp)/maxNumVars > 500)
4089  Dp = modGCDFp (fp, gp, cofp, cogp);
4090  else
4091  {
4092  Dp= gcd_poly (fp, gp);
4093  cofp= fp/Dp;
4094  cogp= gp/Dp;
4095  }
4096 #else
4097  Dp= gcd_poly (fp, gp);
4098  cofp= fp/Dp;
4099  cogp= gp/Dp;
4100 #endif
4101  TIMING_END_AND_PRINT (modZ_recursion,
4102  "time for gcd mod p in modular gcd: ");
4103  Dp /=Dp.lc();
4104  Dp *= mapinto (cl);
4105  cofp /= lc (cofp);
4106  cofp *= mapinto (lcf);
4107  cogp /= lc (cogp);
4108  cogp *= mapinto (lcg);
4109  setCharacteristic( 0 );
4110  dp_deg=totaldegree(Dp);
4111  if ( dp_deg == 0 )
4112  {
4113  //printf(" -> 1\n");
4114  return CanonicalForm(gcdcfcg);
4115  }
4116  if ( q.isZero() )
4117  {
4118  D = mapinto( Dp );
4119  cof= mapinto (cofp);
4120  cog= mapinto (cogp);
4121  d_deg=dp_deg;
4122  q = p;
4123  Dn= balance_p (D, p);
4124  cofn= balance_p (cof, p);
4125  cogn= balance_p (cog, p);
4126  }
4127  else
4128  {
4129  if ( dp_deg == d_deg )
4130  {
4131  chineseRemainder( D, q, mapinto( Dp ), p, newD, newq );
4132  chineseRemainder( cof, q, mapinto (cofp), p, newCof, newq);
4133  chineseRemainder( cog, q, mapinto (cogp), p, newCog, newq);
4134  cof= newCof;
4135  cog= newCog;
4136  newqh= newq/2;
4137  Dn= balance_p (newD, newq, newqh);
4138  cofn= balance_p (newCof, newq, newqh);
4139  cogn= balance_p (newCog, newq, newqh);
4140  if (test != Dn) //balance_p (newD, newq))
4141  test= balance_p (newD, newq);
4142  else
4143  equal= true;
4144  q = newq;
4145  D = newD;
4146  }
4147  else if ( dp_deg < d_deg )
4148  {
4149  //printf(" were all bad, try more\n");
4150  // all previous p's are bad primes
4151  q = p;
4152  D = mapinto( Dp );
4153  cof= mapinto (cof);
4154  cog= mapinto (cog);
4155  d_deg=dp_deg;
4156  test= 0;
4157  equal= false;
4158  Dn= balance_p (D, p);
4159  cofn= balance_p (cof, p);
4160  cogn= balance_p (cog, p);
4161  }
4162  else
4163  {
4164  //printf(" was bad, try more\n");
4165  }
4166  //else dp_deg > d_deg: bad prime
4167  }
4168  if ( i >= 0 )
4169  {
4170  cDn= icontent (Dn);
4171  Dn /= cDn;
4172  cofn /= cl/cDn;
4173  //cofn /= icontent (cofn);
4174  cogn /= cl/cDn;
4175  //cogn /= icontent (cogn);
4176  TIMING_START (modZ_termination);
4177  if ((terminationTest (f,g, cofn, cogn, Dn)) ||
4178  ((equal || q > b) && fdivides (Dn, f) && fdivides (Dn, g)))
4179  {
4180  TIMING_END_AND_PRINT (modZ_termination,
4181  "time for successful termination in modular gcd: ");
4182  //printf(" -> success\n");
4183  return Dn*gcdcfcg;
4184  }
4185  TIMING_END_AND_PRINT (modZ_termination,
4186  "time for unsuccessful termination in modular gcd: ");
4187  equal= false;
4188  //else: try more primes
4189  }
4190  else
4191  { // try other method
4192  //printf("try other gcd\n");
4194  D=gcd_poly( f, g );
4196  return D*gcdcfcg;
4197  }
4198  }
4199 }
4200 
4201 
4202 #endif
int both_zero
Definition: cfGcdAlgExt.cc:71
TIMING_END_AND_PRINT(fac_alg_resultant, "time to compute resultant0: ")
CFMatrix * convertNmod_mat_t2FacCFMatrix(const nmod_mat_t m)
conversion of a FLINT matrix over Z/p to a factory matrix
CanonicalForm power(const CanonicalForm &f, int n)
exponentiation
List< CanonicalForm > CFList
int cf_getNumBigPrimes()
Definition: cf_primes.cc:45
CFArray solveSystemFq(const CFMatrix &M, const CFArray &L, const Variable &alpha)
Definition: cfModGcd.cc:1857
static CanonicalForm mapDown(const CanonicalForm &F, const Variable &alpha, const CanonicalForm &G, CFList &source, CFList &dest)
the CanonicalForm G is the output of map_up, returns F considered as an element over ...
Definition: cf_map_ext.cc:90
CanonicalForm monicSparseInterpol(const CanonicalForm &F, const CanonicalForm &G, const CanonicalForm &skeleton, const Variable &alpha, bool &fail, CFArray *&coeffMonoms, CFArray &Monoms)
Definition: cfModGcd.cc:2140
generate random elements in GF
Definition: cf_random.h:31
CanonicalForm icontent(const CanonicalForm &f)
CanonicalForm icontent ( const CanonicalForm & f )
Definition: cf_gcd.cc:71
int j
Definition: facHensel.cc:105
CF_NO_INLINE bool isOne() const
CF_INLINE bool CanonicalForm::isOne, isZero () const.
Definition: cf_inline.cc:354
#define D(A)
Definition: gentable.cc:129
Conversion to and from NTL.
This file defines functions for Hensel lifting.
static CanonicalForm bound(const CFMatrix &M)
Definition: cf_linsys.cc:460
CanonicalForm newCog
Definition: cfModGcd.cc:4070
This file defines functions for conversion to FLINT (www.flintlib.org) and back.
CanonicalForm fp
Definition: cfModGcd.cc:4043
CF_NO_INLINE CanonicalForm mod(const CanonicalForm &, const CanonicalForm &)
Definition: cf_inline.cc:564
CanonicalForm sparseGCDFq(const CanonicalForm &F, const CanonicalForm &G, const Variable &alpha, CFList &l, bool &topLevel)
Definition: cfModGcd.cc:3071
#define DELETE_ARRAY(P)
Definition: cf_defs.h:49
void convertFacCFMatrix2nmod_mat_t(nmod_mat_t M, const CFMatrix &m)
conversion of a factory matrix over Z/p to a nmod_mat_t
void prune1(const Variable &alpha)
Definition: variable.cc:289
void Off(int sw)
switches
Matrix< CanonicalForm > CFMatrix
CanonicalForm cd(bCommonDen(FF))
Definition: cfModGcd.cc:4030
CanonicalForm LCF
Definition: facFactorize.cc:53
This file provides functions to compute the Newton polygon of a bivariate polynomial.
some useful template functions.
functions to print debug output
template CanonicalForm tmax(const CanonicalForm &, const CanonicalForm &)
TIMING_START(fac_alg_resultant)
CFList evalPoints(const CanonicalForm &F, CFList &eval, const Variable &alpha, CFList &list, const bool &GF, bool &fail)
evaluation point search for multivariate factorization, looks for a (F.level() - 1)-tuple such that t...
static Variable chooseExtension(const Variable &alpha)
Definition: cfModGcd.cc:422
CanonicalForm primitiveElement(const Variable &alpha, Variable &beta, bool &fail)
determine a primitive element of , is a primitive element of a field which is isomorphic to ...
Definition: cf_map_ext.cc:310
f
Definition: cfModGcd.cc:4022
factory&#39;s class for variables
Definition: factory.h:117
void mult(CFList &L1, const CFList &L2)
multiply two lists componentwise
Definition: cfModGcd.cc:2117
int dp_deg
Definition: cfModGcd.cc:4019
static CanonicalForm randomElement(const CanonicalForm &F, const Variable &alpha, CFList &list, bool &fail)
compute a random element a of , s.t. F(a) , F is a univariate polynomial, returns fail if there are...
Definition: cfModGcd.cc:381
CF_NO_INLINE CanonicalForm coeff() const
get the current coefficient
void setMipo(const Variable &alpha, const CanonicalForm &mipo)
Definition: variable.cc:219
CF_NO_INLINE bool isZero() const
Definition: cf_inline.cc:372
CFList *& Aeval
<[in] poly
Definition: facFactorize.h:31
cl
Definition: cfModGcd.cc:4041
int * degsg
Definition: cfEzgcd.cc:53
CFArray evaluateMonom(const CanonicalForm &F, const CFList &evalPoints)
Definition: cfModGcd.cc:1933
int * degrees(const CanonicalForm &f, int *degs=0)
int * degrees ( const CanonicalForm & f, int * degs )
Definition: cf_ops.cc:493
CFArray solveGeneralVandermonde(const CFArray &M, const CFArray &A)
Definition: cfModGcd.cc:1614
generate random evaluation points
static CanonicalForm GFRandomElement(const CanonicalForm &F, CFList &list, bool &fail)
compute a random element a of GF, s.t. F(a) , F is a univariate polynomial, returns fail if there ar...
Definition: cfModGcd.cc:807
CanonicalForm sparseGCDFp(const CanonicalForm &F, const CanonicalForm &G, bool &topLevel, CFList &l)
Definition: cfModGcd.cc:3503
CFFListIterator iter
Definition: facAbsBiFact.cc:54
CanonicalForm nonMonicSparseInterpol(const CanonicalForm &F, const CanonicalForm &G, const CanonicalForm &skeleton, const Variable &alpha, bool &fail, CFArray *&coeffMonoms, CFArray &Monoms)
Definition: cfModGcd.cc:2424
gmp_float log(const gmp_float &a)
Definition: mpr_complex.cc:344
factory&#39;s main class
Definition: canonicalform.h:77
char gf_name
Definition: gfops.cc:52
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
static const double log2exp
Definition: cfModGcd.cc:89
int d_deg
Definition: cfModGcd.cc:4019
static CanonicalForm newtonInterp(const CanonicalForm &alpha, const CanonicalForm &u, const CanonicalForm &newtonPoly, const CanonicalForm &oldInterPoly, const Variable &x)
Newton interpolation - Incremental algorithm. Given a list of values alpha_i and a list of polynomial...
Definition: cfModGcd.cc:366
CanonicalForm convertNTLzzpX2CF(const zz_pX &poly, const Variable &x)
Definition: NTLconvert.cc:248
Array< CanonicalForm > CFArray
CanonicalForm cDn
Definition: cfModGcd.cc:4070
g
Definition: cfModGcd.cc:4031
const CanonicalForm CFMap CFMap bool topLevel
Definition: cfGcdAlgExt.cc:57
CanonicalForm modGCDGF(const CanonicalForm &F, const CanonicalForm &G, CanonicalForm &coF, CanonicalForm &coG, CFList &l, bool &topLevel)
GCD of F and G over GF, based on Alg. 7.2. as described in "Algorithms for Computer Algebra" by Gedde...
Definition: cfModGcd.cc:860
int k
Definition: cfEzgcd.cc:92
Variable alpha
Definition: facAbsBiFact.cc:52
coprimality check and change of representation mod n
void eval(const CanonicalForm &A, const CanonicalForm &B, CanonicalForm &Aeval, CanonicalForm &Beval, const CFList &L)
Definition: cfModGcd.cc:2126
CanonicalForm Lc(const CanonicalForm &f)
CanonicalForm getMipo(const Variable &alpha, const Variable &x)
Definition: variable.cc:207
CanonicalForm generate() const
Definition: cf_random.cc:145
void setCharacteristic(int c)
Definition: cf_char.cc:23
CFMatrix * convertNTLmat_zz_p2FacCFMatrix(const mat_zz_p &m)
Definition: NTLconvert.cc:1182
int getCharacteristic()
Definition: cf_char.cc:51
template bool find(const List< CanonicalForm > &, const CanonicalForm &)
Rational abs(const Rational &a)
Definition: GMPrat.cc:439
void prune(Variable &alpha)
Definition: variable.cc:261
map polynomials
int size() const
Definition: ftmpl_array.cc:92
CanonicalForm content(const CanonicalForm &)
CanonicalForm content ( const CanonicalForm & f )
Definition: cf_gcd.cc:180
This file implements functions to map between extensions of finite fields.
const CanonicalForm & GG
Definition: cfModGcd.cc:4017
int minCommonDeg
Definition: cfModGcd.cc:4045
CanonicalForm lc(const CanonicalForm &f)
CanonicalForm lcg
Definition: cfModGcd.cc:4038
const signed long floor(const ampf< Precision > &x)
Definition: amp.h:774
Variable rootOf(const CanonicalForm &, char name='@')
returns a symbolic root of polynomial with name name Use it to define algebraic variables ...
Definition: variable.cc:162
CFArray readOffSolution(const CFMatrix &M, const long rk)
M in row echolon form, rk rank of M.
Definition: cfModGcd.cc:1670
#define NEW_ARRAY(T, N)
Definition: cf_defs.h:48
#define M
Definition: sirandom.c:24
CanonicalForm cogp
Definition: cfModGcd.cc:4070
bool equal
Definition: cfModGcd.cc:4067
CanonicalForm b
Definition: cfModGcd.cc:4044
CanonicalForm swapvar(const CanonicalForm &, const Variable &, const Variable &)
swapvar() - swap variables x1 and x2 in f.
Definition: cf_ops.cc:168
CFArray getMonoms(const CanonicalForm &F)
extract monomials of F, parts in algebraic variable are considered coefficients
Definition: cfModGcd.cc:1898
int i
Definition: cfModGcd.cc:4019
void removeLast()
Definition: ftmpl_list.cc:317
static CanonicalForm mapUp(const Variable &alpha, const Variable &beta)
and is a primitive element, returns the image of
Definition: cf_map_ext.cc:67
static CanonicalForm uni_lcoeff(const CanonicalForm &F)
compute the leading coefficient of F, where F is considered as an element of , order on is dp...
Definition: cfModGcd.cc:341
bool inBaseDomain() const
const CanonicalForm CFMap CFMap & N
Definition: cfEzgcd.cc:48
mat_zz_p * convertFacCFMatrix2NTLmat_zz_p(const CFMatrix &m)
Definition: NTLconvert.cc:1166
long gaussianElimFq(CFMatrix &M, CFArray &L, const Variable &alpha)
Definition: cfModGcd.cc:1766
CanonicalForm generate() const
Definition: cf_random.cc:66
const signed long ceil(const ampf< Precision > &x)
Definition: amp.h:789
int * degsf
Definition: cfEzgcd.cc:52
int status int void * buf
Definition: si_signals.h:59
CanonicalForm modGCDFp(const CanonicalForm &F, const CanonicalForm &G, CanonicalForm &coF, CanonicalForm &coG, bool &topLevel, CFList &l)
Definition: cfModGcd.cc:1206
int level() const
Definition: factory.h:134
CanonicalForm GFMapDown(const CanonicalForm &F, int k)
maps a polynomial over to a polynomial over , d needs to be a multiple of k
Definition: cf_map_ext.cc:243
bool isUnivariate() const
T & getItem() const
Definition: ftmpl_list.cc:431
CanonicalForm lcf
Definition: cfModGcd.cc:4038
#define A
Definition: sirandom.c:23
int rows() const
Definition: ftmpl_matrix.h:43
CanonicalForm generate() const
Definition: cf_random.cc:56
int columns() const
Definition: ftmpl_matrix.h:44
CanonicalForm newCof
Definition: cfModGcd.cc:4070
int m
Definition: cfEzgcd.cc:121
void On(int sw)
switches
Iterators for CanonicalForm&#39;s.
generate random elements in F_p
Definition: cf_random.h:43
CanonicalForm mapinto(const CanonicalForm &f)
CanonicalForm extractContents(const CanonicalForm &F, const CanonicalForm &G, CanonicalForm &contentF, CanonicalForm &contentG, CanonicalForm &ppF, CanonicalForm &ppG, const int d)
Definition: cfModGcd.cc:313
CanonicalForm cogn
Definition: cfModGcd.cc:4070
long gaussianElimFp(CFMatrix &M, CFArray &L)
Definition: cfModGcd.cc:1721
T getLast() const
Definition: ftmpl_list.cc:309
CFArray solveVandermonde(const CFArray &M, const CFArray &A)
Definition: cfModGcd.cc:1561
generate random integers, random elements of finite fields
CFArray evaluate(const CFArray &A, const CFList &evalPoints)
Definition: cfModGcd.cc:1972
CanonicalForm H
Definition: facAbsFact.cc:64
static const int SW_USE_CHINREM_GCD
set to 1 to use modular gcd over Z
Definition: cf_defs.h:38
Variable mvar() const
mvar() returns the main variable of CO or Variable() if CO is in a base domain.
CFList tmp2
Definition: facFqBivar.cc:70
CFList evaluationPoints(const CanonicalForm &F, const CanonicalForm &G, CanonicalForm &Feval, CanonicalForm &Geval, const CanonicalForm &LCF, const bool &GF, const Variable &alpha, bool &fail, CFList &list)
Definition: cfModGcd.cc:1989
declarations of higher level algorithms.
CanonicalForm bCommonDen(const CanonicalForm &f)
CanonicalForm bCommonDen ( const CanonicalForm & f )
CanonicalForm gp
Definition: cfModGcd.cc:4043
CanonicalForm cof
Definition: cfModGcd.cc:4070
CanonicalForm modGCDFq(const CanonicalForm &F, const CanonicalForm &G, CanonicalForm &coF, CanonicalForm &coG, Variable &alpha, CFList &l, bool &topLevel)
GCD of F and G over , l and topLevel are only used internally, output is monic based on Alg...
Definition: cfModGcd.cc:467
int size(const CanonicalForm &f, const Variable &v)
int size ( const CanonicalForm & f, const Variable & v )
Definition: cf_ops.cc:600
CanonicalForm maxNorm(const CanonicalForm &f)
CanonicalForm maxNorm ( const CanonicalForm & f )
const CanonicalForm & G
Definition: cfModGcd.cc:66
class to iterate through CanonicalForm&#39;s
Definition: cf_iter.h:44
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 )
int getNumVars(const CanonicalForm &f)
int getNumVars ( const CanonicalForm & f )
Definition: cf_ops.cc:314
const Variable & v
< [in] a sqrfree bivariate poly
Definition: facBivar.h:37
CanonicalForm modGCDZ(const CanonicalForm &FF, const CanonicalForm &GG)
modular GCD over Z
CanonicalForm Feval
Definition: facAbsFact.cc:64
generate random irreducible univariate polynomials
CanonicalForm cf
Definition: cfModGcd.cc:4024
CanonicalForm balance_p(const CanonicalForm &f, const CanonicalForm &q, const CanonicalForm &qh)
same as balance_p ( const CanonicalForm & f, const CanonicalForm & q ) but qh= q/2 is provided...
Definition: cfGcdUtil.cc:227
void newpair(const Variable &v, const CanonicalForm &s)
void CFMap::newpair ( const Variable & v, const CanonicalForm & s )
Definition: cf_map.cc:120
#define NULL
Definition: omList.c:10
generate random elements in F_p(alpha)
Definition: cf_random.h:70
int length() const
Definition: ftmpl_list.cc:273
int ipower(int b, int m)
int ipower ( int b, int m )
Definition: cf_util.cc:25
CFArray solveSystemFp(const CFMatrix &M, const CFArray &L)
Definition: cfModGcd.cc:1805
access to prime tables
zz_pX convertFacCF2NTLzzpX(const CanonicalForm &f)
Definition: NTLconvert.cc:100
CanonicalForm Dn
Definition: cfModGcd.cc:4037
CanonicalForm mipo
Definition: facAlgExt.cc:57
CF_NO_INLINE int hasTerms() const
check if iterator has reached the end of CanonicalForm
int f_zero
Definition: cfEzgcd.cc:62
CanonicalForm cg
Definition: cfModGcd.cc:4024
CanonicalForm GFMapUp(const CanonicalForm &F, int k)
maps a polynomial over to a polynomial over , d needs to be a multiple of k
Definition: cf_map_ext.cc:207
b *CanonicalForm B
Definition: facBivar.cc:52
int gcd(int a, int b)
Definition: walkSupport.cc:836
int cf_getBigPrime(int i)
Definition: cf_primes.cc:39
mat_zz_pE * convertFacCFMatrix2NTLmat_zz_pE(const CFMatrix &m)
Definition: NTLconvert.cc:1195
CFList tmp1
Definition: facFqBivar.cc:70
static CanonicalForm uni_content(const CanonicalForm &F)
compute the content of F, where F is considered as an element of
Definition: cfModGcd.cc:283
int getGFDegree()
Definition: cf_char.cc:56
int g_zero
Definition: cfEzgcd.cc:63
Variable x
Definition: cfModGcd.cc:4023
const CanonicalForm const CanonicalForm & coF
Definition: cfModGcd.cc:66
int totaldegree(const CanonicalForm &f)
int totaldegree ( const CanonicalForm & f )
Definition: cf_ops.cc:523
modular and sparse modular GCD algorithms over finite fields and Z.
static CanonicalForm FpRandomElement(const CanonicalForm &F, CFList &list, bool &fail)
Definition: cfModGcd.cc:1159
const CanonicalForm CFMap CFMap int & both_non_zero
Definition: cfEzgcd.cc:50
int maxNumVars
Definition: cfModGcd.cc:4071
int level() const
level() returns the level of CO.
CanonicalForm randomIrredpoly(int i, const Variable &x)
computes a random monic irreducible univariate polynomial in x over Fp of degree i via NTL ...
Definition: cf_irred.cc:42
#define ASSERT(expression, message)
Definition: cf_assert.h:99
#define DEBOUTLN(stream, objects)
Definition: debug.h:49
const CanonicalForm const CanonicalForm const CanonicalForm & coG
Definition: cfModGcd.cc:66
long ind2(long arg)
Definition: kutil.cc:4044
CanonicalForm lc() const
CanonicalForm CanonicalForm::lc (), Lc (), LC (), LC ( v ) const.
CanonicalForm genOne() const
bool terminationTest(const CanonicalForm &F, const CanonicalForm &G, const CanonicalForm &coF, const CanonicalForm &coG, const CanonicalForm &cand)
CanonicalForm gcdcfcg
Definition: cfModGcd.cc:4042
void append(const T &)
Definition: ftmpl_list.cc:256
CanonicalForm gcd_poly(const CanonicalForm &f, const CanonicalForm &g)
CanonicalForm gcd_poly ( const CanonicalForm & f, const CanonicalForm & g )
Definition: cf_gcd.cc:94
long fac_NTL_char
Definition: NTLconvert.cc:41
int p
Definition: cfModGcd.cc:4019
int degree(const CanonicalForm &f)
CanonicalForm mapPrimElem(const CanonicalForm &primElem, const Variable &alpha, const Variable &beta)
compute the image of a primitive element of in . We assume .
Definition: cf_map_ext.cc:377
CF_NO_INLINE int exp() const
get the current exponent
CanonicalForm LC(const CanonicalForm &f)
CanonicalForm cofp
Definition: cfModGcd.cc:4070
Rational pow(const Rational &a, int e)
Definition: GMPrat.cc:414
CFMatrix * convertNTLmat_zz_pE2FacCFMatrix(const mat_zz_pE &m, const Variable &alpha)
Definition: NTLconvert.cc:1211
static void evalPoint(const CanonicalForm &F, const CanonicalForm &G, CanonicalForm &FEval, CanonicalForm &GEval, CFGenerator &evalPoint)
CanonicalForm cofn
Definition: cfModGcd.cc:4070
return result
Definition: facAbsBiFact.cc:76
int l
Definition: cfEzgcd.cc:93
Header for factory&#39;s main class CanonicalForm.
const CanonicalForm const CanonicalForm const CanonicalForm const CanonicalForm & cand
Definition: cfModGcd.cc:69
template CanonicalForm tmin(const CanonicalForm &, const CanonicalForm &)
bool inCoeffDomain() const
TIMING_DEFINE_PRINT(gcd_recursion) TIMING_DEFINE_PRINT(newton_interpolation) TIMING_DEFINE_PRINT(termination_test) TIMING_DEFINE_PRINT(ez_p_compress) TIMING_DEFINE_PRINT(ez_p_hensel_lift) TIMING_DEFINE_PRINT(ez_p_eval) TIMING_DEFINE_PRINT(ez_p_content) bool terminationTest(const CanonicalForm &F
CanonicalForm cog
Definition: cfModGcd.cc:4070