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