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