source: git/factory/cfModGcd.cc @ 03c742

spielwiese
Last change on this file since 03c742 was 03c742, checked in by Hans Schoenemann <hannes@…>, 4 years ago
factory: BuildIrred, FLINT/NTL seperation
  • Property mode set to 100644
File size: 115.8 KB
Line 
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"
35#include "templates/ftmpl_functions.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
57TIMING_DEFINE_PRINT(gcd_recursion)
58TIMING_DEFINE_PRINT(newton_interpolation)
59TIMING_DEFINE_PRINT(termination_test)
60TIMING_DEFINE_PRINT(ez_p_compress)
61TIMING_DEFINE_PRINT(ez_p_hensel_lift)
62TIMING_DEFINE_PRINT(ez_p_eval)
63TIMING_DEFINE_PRINT(ez_p_content)
64
65bool
66terminationTest (const CanonicalForm& F, const CanonicalForm& G,
67                 const CanonicalForm& coF, const CanonicalForm& coG,
68                 const CanonicalForm& cand)
69{
70  CanonicalForm LCCand= abs (LC (cand));
71  if (LCCand*abs (LC (coF)) == abs (LC (F)))
72  {
73    if (LCCand*abs (LC (coG)) == abs (LC (G)))
74    {
75      if (abs (cand)*abs (coF) == abs (F))
76      {
77        if (abs (cand)*abs (coG) == abs (G))
78          return true;
79      }
80      return false;
81    }
82    return false;
83  }
84  return false;
85}
86
87#if defined(HAVE_NTL)|| defined(HAVE_FLINT)
88
89static const double log2exp= 1.442695041;
90
91/// compressing two polynomials F and G, M is used for compressing,
92/// N to reverse the compression
93int myCompress (const CanonicalForm& F, const CanonicalForm& G, CFMap & M,
94                CFMap & N, bool topLevel)
95{
96  int n= tmax (F.level(), G.level());
97  int * degsf= NEW_ARRAY(int,n + 1);
98  int * degsg= NEW_ARRAY(int,n + 1);
99
100  for (int i = n; i >= 0; i--)
101    degsf[i]= degsg[i]= 0;
102
103  degsf= degrees (F, degsf);
104  degsg= degrees (G, degsg);
105
106  int both_non_zero= 0;
107  int f_zero= 0;
108  int g_zero= 0;
109  int both_zero= 0;
110
111  if (topLevel)
112  {
113    for (int i= 1; i <= n; i++)
114    {
115      if (degsf[i] != 0 && degsg[i] != 0)
116      {
117        both_non_zero++;
118        continue;
119      }
120      if (degsf[i] == 0 && degsg[i] != 0 && i <= G.level())
121      {
122        f_zero++;
123        continue;
124      }
125      if (degsg[i] == 0 && degsf[i] && i <= F.level())
126      {
127        g_zero++;
128        continue;
129      }
130    }
131
132    if (both_non_zero == 0)
133    {
134      DELETE_ARRAY(degsf);
135      DELETE_ARRAY(degsg);
136      return 0;
137    }
138
139    // map Variables which do not occur in both polynomials to higher levels
140    int k= 1;
141    int l= 1;
142    for (int i= 1; i <= n; i++)
143    {
144      if (degsf[i] != 0 && degsg[i] == 0 && i <= F.level())
145      {
146        if (k + both_non_zero != i)
147        {
148          M.newpair (Variable (i), Variable (k + both_non_zero));
149          N.newpair (Variable (k + both_non_zero), Variable (i));
150        }
151        k++;
152      }
153      if (degsf[i] == 0 && degsg[i] != 0 && i <= G.level())
154      {
155        if (l + g_zero + both_non_zero != i)
156        {
157          M.newpair (Variable (i), Variable (l + g_zero + both_non_zero));
158          N.newpair (Variable (l + g_zero + both_non_zero), Variable (i));
159        }
160        l++;
161      }
162    }
163
164    // sort Variables x_{i} in increasing order of
165    // min(deg_{x_{i}}(f),deg_{x_{i}}(g))
166    int m= tmax (F.level(), G.level());
167    int min_max_deg;
168    k= both_non_zero;
169    l= 0;
170    int i= 1;
171    while (k > 0)
172    {
173      if (degsf [i] != 0 && degsg [i] != 0)
174        min_max_deg= tmax (degsf[i], degsg[i]);
175      else
176        min_max_deg= 0;
177      while (min_max_deg == 0)
178      {
179        i++;
180        if (degsf [i] != 0 && degsg [i] != 0)
181          min_max_deg= tmax (degsf[i], degsg[i]);
182        else
183          min_max_deg= 0;
184      }
185      for (int j= i + 1; j <=  m; j++)
186      {
187        if (degsf[j] != 0 && degsg [j] != 0 &&
188            tmax (degsf[j],degsg[j]) <= min_max_deg)
189        {
190          min_max_deg= tmax (degsf[j], degsg[j]);
191          l= j;
192        }
193      }
194      if (l != 0)
195      {
196        if (l != k)
197        {
198          M.newpair (Variable (l), Variable(k));
199          N.newpair (Variable (k), Variable(l));
200          degsf[l]= 0;
201          degsg[l]= 0;
202          l= 0;
203        }
204        else
205        {
206          degsf[l]= 0;
207          degsg[l]= 0;
208          l= 0;
209        }
210      }
211      else if (l == 0)
212      {
213        if (i != k)
214        {
215          M.newpair (Variable (i), Variable (k));
216          N.newpair (Variable (k), Variable (i));
217          degsf[i]= 0;
218          degsg[i]= 0;
219        }
220        else
221        {
222          degsf[i]= 0;
223          degsg[i]= 0;
224        }
225        i++;
226      }
227      k--;
228    }
229  }
230  else
231  {
232    //arrange Variables such that no gaps occur
233    for (int i= 1; i <= n; i++)
234    {
235      if (degsf[i] == 0 && degsg[i] == 0)
236      {
237        both_zero++;
238        continue;
239      }
240      else
241      {
242        if (both_zero != 0)
243        {
244          M.newpair (Variable (i), Variable (i - both_zero));
245          N.newpair (Variable (i - both_zero), Variable (i));
246        }
247      }
248    }
249  }
250
251  DELETE_ARRAY(degsf);
252  DELETE_ARRAY(degsg);
253
254  return 1;
255}
256
257static inline CanonicalForm
258uni_content (const CanonicalForm & F);
259
260CanonicalForm
261uni_content (const CanonicalForm& F, const Variable& x)
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));
273    CanonicalForm result= uni_content (f);
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$
282static inline CanonicalForm
283uni_content (const CanonicalForm & F)
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
312CanonicalForm
313extractContents (const CanonicalForm& F, const CanonicalForm& G,
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;
322  CanonicalForm result= 1;
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.
340static inline
341CanonicalForm uni_lcoeff (const CanonicalForm& F)
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})
365static inline CanonicalForm
366newtonInterp(const CanonicalForm & alpha, const CanonicalForm & u,
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
380static inline CanonicalForm
381randomElement (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
421static inline
422Variable chooseExtension (const Variable & alpha)
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
443  if (fac_NTL_char != getCharacteristic())
444  {
445    fac_NTL_char= getCharacteristic();
446    zz_p::init (getCharacteristic());
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#ifdef HAVE_NTL
456CanonicalForm
457modGCDFq (const CanonicalForm& F, const CanonicalForm& G,
458                  CanonicalForm& coF, CanonicalForm& coG,
459                  Variable & alpha, CFList& l, bool& topLevel);
460#endif
461
462#ifdef HAVE_NTL // modGCDFq
463CanonicalForm
464modGCDFq (const CanonicalForm& F, const CanonicalForm& G,
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#ifdef HAVE_NTL // primitiveElement, FindRoot
479CanonicalForm
480modGCDFq (const CanonicalForm& F, const CanonicalForm& G,
481                  CanonicalForm& coF, CanonicalForm& coG,
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
819static inline
820CanonicalForm
821GFRandomElement (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
854CanonicalForm
855modGCDGF (const CanonicalForm& F, const CanonicalForm& G,
856        CanonicalForm& coF, CanonicalForm& coG,
857        CFList& l, bool& topLevel);
858
859CanonicalForm
860modGCDGF (const CanonicalForm& F, const CanonicalForm& G, CFList& l,
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
873CanonicalForm
874modGCDGF (const CanonicalForm& F, const CanonicalForm& G,
875        CanonicalForm& coF, CanonicalForm& coG,
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
1171static inline
1172CanonicalForm
1173FpRandomElement (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#ifdef HAVE_NTL // primitiveElement, FindRoot
1206CanonicalForm
1207modGCDFp (const CanonicalForm& F, const CanonicalForm&  G,
1208             CanonicalForm& coF, CanonicalForm& coG,
1209             bool& topLevel, CFList& l);
1210#endif
1211
1212#ifdef HAVE_NTL // primitiveElement, FindRoot
1213CanonicalForm
1214modGCDFp (const CanonicalForm& F, const CanonicalForm& G,
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#ifdef HAVE_NTL // primitiveElement, FindRoot
1224CanonicalForm
1225modGCDFp (const CanonicalForm& F, const CanonicalForm&  G,
1226             CanonicalForm& coF, CanonicalForm& coG,
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;
1385      CanonicalForm mipo;
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
1580CFArray
1581solveVandermonde (const CFArray& M, const CFArray& A)
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
1633CFArray
1634solveGeneralVandermonde (const CFArray& M, const CFArray& A)
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
1689CFArray
1690readOffSolution (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
1711CFArray
1712readOffSolution (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
1740long
1741gaussianElimFp (CFMatrix& M, CFArray& L)
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
1773  N= convertNTLmat_zz_p2FacCFMatrix (*NTLN);
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#ifdef HAVE_NTL //gauss, zz_pE
1786long
1787gaussianElimFq (CFMatrix& M, CFArray& L, const Variable& alpha)
1788{
1789  ASSERT (L.size() <= M.rows(), "dimension exceeded");
1790  CFMatrix *N;
1791  N= new CFMatrix (M.rows(), M.columns() + 1);
1792
1793  for (int i= 1; i <= M.rows(); i++)
1794    for (int j= 1; j <= M.columns(); j++)
1795      (*N) (i, j)= M (i, j);
1796
1797  int j= 1;
1798  for (int i= 0; i < L.size(); i++, j++)
1799    (*N) (j, M.columns() + 1)= L[i];
1800  int p= getCharacteristic ();
1801  if (fac_NTL_char != p)
1802  {
1803    fac_NTL_char= p;
1804    zz_p::init (p);
1805  }
1806  zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
1807  zz_pE::init (NTLMipo);
1808  mat_zz_pE *NTLN= convertFacCFMatrix2NTLmat_zz_pE(*N);
1809  long rk= gauss (*NTLN);
1810
1811  delete N;
1812  N= convertNTLmat_zz_pE2FacCFMatrix (*NTLN, alpha);
1813
1814  delete NTLN;
1815
1816  M= (*N) (1, M.rows(), 1, M.columns());
1817  L= CFArray (M.rows());
1818  for (int i= 0; i < M.rows(); i++)
1819    L[i]= (*N) (i + 1, M.columns() + 1);
1820
1821  delete N;
1822  return rk;
1823}
1824#endif
1825
1826CFArray
1827solveSystemFp (const CFMatrix& M, const CFArray& L)
1828{
1829  ASSERT (L.size() <= M.rows(), "dimension exceeded");
1830  CFMatrix *N;
1831  N= new CFMatrix (M.rows(), M.columns() + 1);
1832
1833  for (int i= 1; i <= M.rows(); i++)
1834    for (int j= 1; j <= M.columns(); j++)
1835      (*N) (i, j)= M (i, j);
1836
1837  int j= 1;
1838  for (int i= 0; i < L.size(); i++, j++)
1839    (*N) (j, M.columns() + 1)= L[i];
1840
1841#ifdef HAVE_FLINT
1842  nmod_mat_t FLINTN;
1843  convertFacCFMatrix2nmod_mat_t (FLINTN, *N);
1844  long rk= nmod_mat_rref (FLINTN);
1845#else
1846  int p= getCharacteristic ();
1847  if (fac_NTL_char != p)
1848  {
1849    fac_NTL_char= p;
1850    zz_p::init (p);
1851  }
1852  mat_zz_p *NTLN= convertFacCFMatrix2NTLmat_zz_p(*N);
1853  long rk= gauss (*NTLN);
1854#endif
1855  delete N;
1856  if (rk != M.columns())
1857  {
1858#ifdef HAVE_FLINT
1859    nmod_mat_clear (FLINTN);
1860#else
1861    delete NTLN;
1862#endif
1863    return CFArray();
1864  }
1865#ifdef HAVE_FLINT
1866  N= convertNmod_mat_t2FacCFMatrix (FLINTN);
1867  nmod_mat_clear (FLINTN);
1868#else
1869  N= convertNTLmat_zz_p2FacCFMatrix (*NTLN);
1870  delete NTLN;
1871#endif
1872  CFArray A= readOffSolution (*N, rk);
1873
1874  delete N;
1875  return A;
1876}
1877
1878#ifdef HAVE_NTL //gauss, zz_pE
1879CFArray
1880solveSystemFq (const CFMatrix& M, const CFArray& L, const Variable& alpha)
1881{
1882  ASSERT (L.size() <= M.rows(), "dimension exceeded");
1883  CFMatrix *N;
1884  N= new CFMatrix (M.rows(), M.columns() + 1);
1885
1886  for (int i= 1; i <= M.rows(); i++)
1887    for (int j= 1; j <= M.columns(); j++)
1888      (*N) (i, j)= M (i, j);
1889  int j= 1;
1890  for (int i= 0; i < L.size(); i++, j++)
1891    (*N) (j, M.columns() + 1)= L[i];
1892  int p= getCharacteristic ();
1893  if (fac_NTL_char != p)
1894  {
1895    fac_NTL_char= p;
1896    zz_p::init (p);
1897  }
1898  zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
1899  zz_pE::init (NTLMipo);
1900  mat_zz_pE *NTLN= convertFacCFMatrix2NTLmat_zz_pE(*N);
1901  long rk= gauss (*NTLN);
1902
1903  delete N;
1904  if (rk != M.columns())
1905  {
1906    delete NTLN;
1907    return CFArray();
1908  }
1909  N= convertNTLmat_zz_pE2FacCFMatrix (*NTLN, alpha);
1910
1911  delete NTLN;
1912
1913  CFArray A= readOffSolution (*N, rk);
1914
1915  delete N;
1916  return A;
1917}
1918#endif
1919#endif
1920
1921CFArray
1922getMonoms (const CanonicalForm& F)
1923{
1924  if (F.inCoeffDomain())
1925  {
1926    CFArray result= CFArray (1);
1927    result [0]= 1;
1928    return result;
1929  }
1930  if (F.isUnivariate())
1931  {
1932    CFArray result= CFArray (size(F));
1933    int j= 0;
1934    for (CFIterator i= F; i.hasTerms(); i++, j++)
1935      result[j]= power (F.mvar(), i.exp());
1936    return result;
1937  }
1938  int numMon= size (F);
1939  CFArray result= CFArray (numMon);
1940  int j= 0;
1941  CFArray recResult;
1942  Variable x= F.mvar();
1943  CanonicalForm powX;
1944  for (CFIterator i= F; i.hasTerms(); i++)
1945  {
1946    powX= power (x, i.exp());
1947    recResult= getMonoms (i.coeff());
1948    for (int k= 0; k < recResult.size(); k++)
1949      result[j+k]= powX*recResult[k];
1950    j += recResult.size();
1951  }
1952  return result;
1953}
1954
1955#ifdef HAVE_NTL
1956CFArray
1957evaluateMonom (const CanonicalForm& F, const CFList& evalPoints)
1958{
1959  if (F.inCoeffDomain())
1960  {
1961    CFArray result= CFArray (1);
1962    result [0]= F;
1963    return result;
1964  }
1965  if (F.isUnivariate())
1966  {
1967    ASSERT (evalPoints.length() == 1,
1968            "expected an eval point with only one component");
1969    CFArray result= CFArray (size(F));
1970    int j= 0;
1971    CanonicalForm evalPoint= evalPoints.getLast();
1972    for (CFIterator i= F; i.hasTerms(); i++, j++)
1973      result[j]= power (evalPoint, i.exp());
1974    return result;
1975  }
1976  int numMon= size (F);
1977  CFArray result= CFArray (numMon);
1978  int j= 0;
1979  CanonicalForm evalPoint= evalPoints.getLast();
1980  CFList buf= evalPoints;
1981  buf.removeLast();
1982  CFArray recResult;
1983  CanonicalForm powEvalPoint;
1984  for (CFIterator i= F; i.hasTerms(); i++)
1985  {
1986    powEvalPoint= power (evalPoint, i.exp());
1987    recResult= evaluateMonom (i.coeff(), buf);
1988    for (int k= 0; k < recResult.size(); k++)
1989      result[j+k]= powEvalPoint*recResult[k];
1990    j += recResult.size();
1991  }
1992  return result;
1993}
1994
1995CFArray
1996evaluate (const CFArray& A, const CFList& evalPoints)
1997{
1998  CFArray result= A.size();
1999  CanonicalForm tmp;
2000  int k;
2001  for (int i= 0; i < A.size(); i++)
2002  {
2003    tmp= A[i];
2004    k= 1;
2005    for (CFListIterator j= evalPoints; j.hasItem(); j++, k++)
2006      tmp= tmp (j.getItem(), k);
2007    result[i]= tmp;
2008  }
2009  return result;
2010}
2011
2012CFList
2013evaluationPoints (const CanonicalForm& F, const CanonicalForm& G,
2014                  CanonicalForm& Feval, CanonicalForm& Geval,
2015                  const CanonicalForm& LCF, const bool& GF,
2016                  const Variable& alpha, bool& fail, CFList& list
2017                 )
2018{
2019  int k= tmax (F.level(), G.level()) - 1;
2020  Variable x= Variable (1);
2021  CFList result;
2022  FFRandom genFF;
2023  GFRandom genGF;
2024  int p= getCharacteristic ();
2025  double bound;
2026  if (alpha != Variable (1))
2027  {
2028    bound= pow ((double) p, (double) degree (getMipo(alpha)));
2029    bound= pow (bound, (double) k);
2030  }
2031  else if (GF)
2032  {
2033    bound= pow ((double) p, (double) getGFDegree());
2034    bound= pow ((double) bound, (double) k);
2035  }
2036  else
2037    bound= pow ((double) p, (double) k);
2038
2039  CanonicalForm random;
2040  int j;
2041  bool zeroOneOccured= false;
2042  bool allEqual= false;
2043  CanonicalForm buf;
2044  do
2045  {
2046    random= 0;
2047    // possible overflow if list.length() does not fit into a int
2048    if (list.length() >= bound)
2049    {
2050      fail= true;
2051      break;
2052    }
2053    for (int i= 0; i < k; i++)
2054    {
2055      if (GF)
2056      {
2057        result.append (genGF.generate());
2058        random += result.getLast()*power (x, i);
2059      }
2060      else if (alpha.level() != 1)
2061      {
2062        AlgExtRandomF genAlgExt (alpha);
2063        result.append (genAlgExt.generate());
2064        random += result.getLast()*power (x, i);
2065      }
2066      else
2067      {
2068        result.append (genFF.generate());
2069        random += result.getLast()*power (x, i);
2070      }
2071      if (result.getLast().isOne() || result.getLast().isZero())
2072        zeroOneOccured= true;
2073    }
2074    if (find (list, random))
2075    {
2076      zeroOneOccured= false;
2077      allEqual= false;
2078      result= CFList();
2079      continue;
2080    }
2081    if (zeroOneOccured)
2082    {
2083      list.append (random);
2084      zeroOneOccured= false;
2085      allEqual= false;
2086      result= CFList();
2087      continue;
2088    }
2089    // no zero at this point
2090    if (k > 1)
2091    {
2092      allEqual= true;
2093      CFIterator iter= random;
2094      buf= iter.coeff();
2095      iter++;
2096      for (; iter.hasTerms(); iter++)
2097        if (buf != iter.coeff())
2098          allEqual= false;
2099    }
2100    if (allEqual)
2101    {
2102      list.append (random);
2103      allEqual= false;
2104      zeroOneOccured= false;
2105      result= CFList();
2106      continue;
2107    }
2108
2109    Feval= F;
2110    Geval= G;
2111    CanonicalForm LCeval= LCF;
2112    j= 1;
2113    for (CFListIterator i= result; i.hasItem(); i++, j++)
2114    {
2115      Feval= Feval (i.getItem(), j);
2116      Geval= Geval (i.getItem(), j);
2117      LCeval= LCeval (i.getItem(), j);
2118    }
2119
2120    if (LCeval.isZero())
2121    {
2122      if (!find (list, random))
2123        list.append (random);
2124      zeroOneOccured= false;
2125      allEqual= false;
2126      result= CFList();
2127      continue;
2128    }
2129
2130    if (list.length() >= bound)
2131    {
2132      fail= true;
2133      break;
2134    }
2135  } while (find (list, random));
2136
2137  return result;
2138}
2139
2140/// multiply two lists componentwise
2141void mult (CFList& L1, const CFList& L2)
2142{
2143  ASSERT (L1.length() == L2.length(), "lists of the same size expected");
2144
2145  CFListIterator j= L2;
2146  for (CFListIterator i= L1; i.hasItem(); i++, j++)
2147    i.getItem() *= j.getItem();
2148}
2149
2150void eval (const CanonicalForm& A, const CanonicalForm& B, CanonicalForm& Aeval,
2151           CanonicalForm& Beval, const CFList& L)
2152{
2153  Aeval= A;
2154  Beval= B;
2155  int j= 1;
2156  for (CFListIterator i= L; i.hasItem(); i++, j++)
2157  {
2158    Aeval= Aeval (i.getItem(), j);
2159    Beval= Beval (i.getItem(), j);
2160  }
2161}
2162
2163#ifdef HAVE_NTL // primitiveElement,FindRoot
2164CanonicalForm
2165monicSparseInterpol (const CanonicalForm& F, const CanonicalForm& G,
2166                     const CanonicalForm& skeleton, const Variable& alpha,
2167                     bool& fail, CFArray*& coeffMonoms, CFArray& Monoms
2168                    )
2169{
2170  CanonicalForm A= F;
2171  CanonicalForm B= G;
2172  if (F.isZero() && degree(G) > 0) return G/Lc(G);
2173  else if (G.isZero() && degree (F) > 0) return F/Lc(F);
2174  else if (F.isZero() && G.isZero()) return F.genOne();
2175  if (F.inBaseDomain() || G.inBaseDomain()) return F.genOne();
2176  if (F.isUnivariate() && fdivides(F, G)) return F/Lc(F);
2177  if (G.isUnivariate() && fdivides(G, F)) return G/Lc(G);
2178  if (F == G) return F/Lc(F);
2179
2180  ASSERT (degree (A, 1) == 0, "expected degree (F, 1) == 0");
2181  ASSERT (degree (B, 1) == 0, "expected degree (G, 1) == 0");
2182
2183  CFMap M,N;
2184  int best_level= myCompress (A, B, M, N, false);
2185
2186  if (best_level == 0)
2187    return B.genOne();
2188
2189  A= M(A);
2190  B= M(B);
2191
2192  Variable x= Variable (1);
2193
2194  //univariate case
2195  if (A.isUnivariate() && B.isUnivariate())
2196    return N (gcd (A, B));
2197
2198  CanonicalForm skel= M(skeleton);
2199  CanonicalForm cA, cB;    // content of A and B
2200  CanonicalForm ppA, ppB;    // primitive part of A and B
2201  CanonicalForm gcdcAcB;
2202  cA = uni_content (A);
2203  cB = uni_content (B);
2204  gcdcAcB= gcd (cA, cB);
2205  ppA= A/cA;
2206  ppB= B/cB;
2207
2208  CanonicalForm lcA, lcB;  // leading coefficients of A and B
2209  CanonicalForm gcdlcAlcB;
2210  lcA= uni_lcoeff (ppA);
2211  lcB= uni_lcoeff (ppB);
2212
2213  if (fdivides (lcA, lcB))
2214  {
2215    if (fdivides (A, B))
2216      return F/Lc(F);
2217  }
2218  if (fdivides (lcB, lcA))
2219  {
2220    if (fdivides (B, A))
2221      return G/Lc(G);
2222  }
2223
2224  gcdlcAlcB= gcd (lcA, lcB);
2225  int skelSize= size (skel, skel.mvar());
2226
2227  int j= 0;
2228  int biggestSize= 0;
2229
2230  for (CFIterator i= skel; i.hasTerms(); i++, j++)
2231    biggestSize= tmax (biggestSize, size (i.coeff()));
2232
2233  CanonicalForm g, Aeval, Beval;
2234
2235  CFList evalPoints;
2236  bool evalFail= false;
2237  CFList list;
2238  bool GF= false;
2239  CanonicalForm LCA= LC (A);
2240  CanonicalForm tmp;
2241  CFArray gcds= CFArray (biggestSize);
2242  CFList * pEvalPoints= new CFList [biggestSize];
2243  Variable V_buf= alpha, V_buf4= alpha;
2244  CFList source, dest;
2245  CanonicalForm prim_elem, im_prim_elem;
2246  CanonicalForm prim_elem_alpha, im_prim_elem_alpha;
2247  for (int i= 0; i < biggestSize; i++)
2248  {
2249    if (i == 0)
2250      evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf, evalFail,
2251                                    list);
2252    else
2253    {
2254      mult (evalPoints, pEvalPoints [0]);
2255      eval (A, B, Aeval, Beval, evalPoints);
2256    }
2257
2258    if (evalFail)
2259    {
2260      if (V_buf.level() != 1)
2261      {
2262        do
2263        {
2264          Variable V_buf3= V_buf;
2265          V_buf= chooseExtension (V_buf);
2266          source= CFList();
2267          dest= CFList();
2268
2269          bool prim_fail= false;
2270          Variable V_buf2;
2271          prim_elem= primitiveElement (V_buf4, V_buf2, prim_fail);
2272          if (V_buf4 == alpha && alpha.level() != 1)
2273            prim_elem_alpha= prim_elem;
2274
2275          ASSERT (!prim_fail, "failure in integer factorizer");
2276          if (prim_fail)
2277            ; //ERROR
2278          else
2279            im_prim_elem= mapPrimElem (prim_elem, V_buf4, V_buf);
2280
2281          DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf));
2282          DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf2));
2283
2284          if (V_buf4 == alpha && alpha.level() != 1)
2285            im_prim_elem_alpha= im_prim_elem;
2286          else if (alpha.level() != 1)
2287            im_prim_elem_alpha= mapUp (im_prim_elem_alpha, V_buf4, V_buf,
2288                                       prim_elem, im_prim_elem, source, dest);
2289
2290          for (CFListIterator j= list; j.hasItem(); j++)
2291            j.getItem()= mapUp (j.getItem(), V_buf4, V_buf, prim_elem,
2292                                im_prim_elem, source, dest);
2293          for (int k= 0; k < i; k++)
2294          {
2295            for (CFListIterator j= pEvalPoints[k]; j.hasItem(); j++)
2296              j.getItem()= mapUp (j.getItem(), V_buf4, V_buf, prim_elem,
2297                                  im_prim_elem, source, dest);
2298            gcds[k]= mapUp (gcds[k], V_buf4, V_buf, prim_elem, im_prim_elem,
2299                            source, dest);
2300          }
2301
2302          if (alpha.level() != 1)
2303          {
2304            A= mapUp (A, V_buf4, V_buf, prim_elem, im_prim_elem, source,dest);
2305            B= mapUp (B, V_buf4, V_buf, prim_elem, im_prim_elem, source,dest);
2306          }
2307          V_buf4= V_buf;
2308          evalFail= false;
2309          evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf,
2310                                        evalFail, list);
2311        } while (evalFail);
2312      }
2313      else
2314      {
2315        CanonicalForm mipo;
2316        int deg= 2;
2317        bool initialized= false;
2318        do
2319        {
2320          mipo= randomIrredpoly (deg, x);
2321          if (initialized)
2322            setMipo (V_buf, mipo);
2323          else
2324            V_buf= rootOf (mipo);
2325          evalFail= false;
2326          initialized= true;
2327          evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf,
2328                                        evalFail, list);
2329          deg++;
2330        } while (evalFail);
2331        V_buf4= V_buf;
2332      }
2333    }
2334
2335    g= gcd (Aeval, Beval);
2336    g /= Lc (g);
2337
2338    if (degree (g) != degree (skel) || (skelSize != size (g)))
2339    {
2340      delete[] pEvalPoints;
2341      fail= true;
2342      if (alpha.level() != 1 && V_buf != alpha)
2343        prune1 (alpha);
2344      return 0;
2345    }
2346    CFIterator l= skel;
2347    for (CFIterator k= g; k.hasTerms(); k++, l++)
2348    {
2349      if (k.exp() != l.exp())
2350      {
2351        delete[] pEvalPoints;
2352        fail= true;
2353        if (alpha.level() != 1 && V_buf != alpha)
2354          prune1 (alpha);
2355        return 0;
2356      }
2357    }
2358    pEvalPoints[i]= evalPoints;
2359    gcds[i]= g;
2360
2361    tmp= 0;
2362    int j= 0;
2363    for (CFListIterator k= evalPoints; k.hasItem(); k++, j++)
2364      tmp += k.getItem()*power (x, j);
2365    list.append (tmp);
2366  }
2367
2368  if (Monoms.size() == 0)
2369    Monoms= getMonoms (skel);
2370
2371  coeffMonoms= new CFArray [skelSize];
2372  j= 0;
2373  for (CFIterator i= skel; i.hasTerms(); i++, j++)
2374    coeffMonoms[j]= getMonoms (i.coeff());
2375
2376  CFArray* pL= new CFArray [skelSize];
2377  CFArray* pM= new CFArray [skelSize];
2378  for (int i= 0; i < biggestSize; i++)
2379  {
2380    CFIterator l= gcds [i];
2381    evalPoints= pEvalPoints [i];
2382    for (int k= 0; k < skelSize; k++, l++)
2383    {
2384      if (i == 0)
2385        pL[k]= CFArray (biggestSize);
2386      pL[k] [i]= l.coeff();
2387
2388      if (i == 0)
2389        pM[k]= evaluate (coeffMonoms [k], evalPoints);
2390    }
2391  }
2392
2393  CFArray solution;
2394  CanonicalForm result= 0;
2395  int ind= 0;
2396  CFArray bufArray;
2397  CFMatrix Mat;
2398  for (int k= 0; k < skelSize; k++)
2399  {
2400    if (biggestSize != coeffMonoms[k].size())
2401    {
2402      bufArray= CFArray (coeffMonoms[k].size());
2403      for (int i= 0; i < coeffMonoms[k].size(); i++)
2404        bufArray [i]= pL[k] [i];
2405      solution= solveGeneralVandermonde (pM [k], bufArray);
2406    }
2407    else
2408      solution= solveGeneralVandermonde (pM [k], pL[k]);
2409
2410    if (solution.size() == 0)
2411    {
2412      delete[] pEvalPoints;
2413      delete[] pM;
2414      delete[] pL;
2415      delete[] coeffMonoms;
2416      fail= true;
2417      if (alpha.level() != 1 && V_buf != alpha)
2418        prune1 (alpha);
2419      return 0;
2420    }
2421    for (int l= 0; l < solution.size(); l++)
2422      result += solution[l]*Monoms [ind + l];
2423    ind += solution.size();
2424  }
2425
2426  delete[] pEvalPoints;
2427  delete[] pM;
2428  delete[] pL;
2429  delete[] coeffMonoms;
2430
2431  if (alpha.level() != 1 && V_buf != alpha)
2432  {
2433    CFList u, v;
2434    result= mapDown (result, prim_elem_alpha, im_prim_elem_alpha, alpha, u, v);
2435    prune1 (alpha);
2436  }
2437
2438  result= N(result);
2439  if (fdivides (result, F) && fdivides (result, G))
2440    return result;
2441  else
2442  {
2443    fail= true;
2444    return 0;
2445  }
2446}
2447#endif
2448
2449#ifdef HAVE_NTL // primitiveElement, FindRoot
2450CanonicalForm
2451nonMonicSparseInterpol (const CanonicalForm& F, const CanonicalForm& G,
2452                        const CanonicalForm& skeleton, const Variable& alpha,
2453                        bool& fail, CFArray*& coeffMonoms, CFArray& Monoms
2454                       )
2455{
2456  CanonicalForm A= F;
2457  CanonicalForm B= G;
2458  if (F.isZero() && degree(G) > 0) return G/Lc(G);
2459  else if (G.isZero() && degree (F) > 0) return F/Lc(F);
2460  else if (F.isZero() && G.isZero()) return F.genOne();
2461  if (F.inBaseDomain() || G.inBaseDomain()) return F.genOne();
2462  if (F.isUnivariate() && fdivides(F, G)) return F/Lc(F);
2463  if (G.isUnivariate() && fdivides(G, F)) return G/Lc(G);
2464  if (F == G) return F/Lc(F);
2465
2466  ASSERT (degree (A, 1) == 0, "expected degree (F, 1) == 0");
2467  ASSERT (degree (B, 1) == 0, "expected degree (G, 1) == 0");
2468
2469  CFMap M,N;
2470  int best_level= myCompress (A, B, M, N, false);
2471
2472  if (best_level == 0)
2473    return B.genOne();
2474
2475  A= M(A);
2476  B= M(B);
2477
2478  Variable x= Variable (1);
2479
2480  //univariate case
2481  if (A.isUnivariate() && B.isUnivariate())
2482    return N (gcd (A, B));
2483
2484  CanonicalForm skel= M(skeleton);
2485
2486  CanonicalForm cA, cB;    // content of A and B
2487  CanonicalForm ppA, ppB;    // primitive part of A and B
2488  CanonicalForm gcdcAcB;
2489  cA = uni_content (A);
2490  cB = uni_content (B);
2491  gcdcAcB= gcd (cA, cB);
2492  ppA= A/cA;
2493  ppB= B/cB;
2494
2495  CanonicalForm lcA, lcB;  // leading coefficients of A and B
2496  CanonicalForm gcdlcAlcB;
2497  lcA= uni_lcoeff (ppA);
2498  lcB= uni_lcoeff (ppB);
2499
2500  if (fdivides (lcA, lcB))
2501  {
2502    if (fdivides (A, B))
2503      return F/Lc(F);
2504  }
2505  if (fdivides (lcB, lcA))
2506  {
2507    if (fdivides (B, A))
2508      return G/Lc(G);
2509  }
2510
2511  gcdlcAlcB= gcd (lcA, lcB);
2512  int skelSize= size (skel, skel.mvar());
2513
2514  int j= 0;
2515  int biggestSize= 0;
2516  int bufSize;
2517  int numberUni= 0;
2518  for (CFIterator i= skel; i.hasTerms(); i++, j++)
2519  {
2520    bufSize= size (i.coeff());
2521    biggestSize= tmax (biggestSize, bufSize);
2522    numberUni += bufSize;
2523  }
2524  numberUni--;
2525  numberUni= (int) ceil ( (double) numberUni / (skelSize - 1));
2526  biggestSize= tmax (biggestSize , numberUni);
2527
2528  numberUni= biggestSize + size (LC(skel)) - 1;
2529  int biggestSize2= tmax (numberUni, biggestSize);
2530
2531  CanonicalForm g, Aeval, Beval;
2532
2533  CFList evalPoints;
2534  CFArray coeffEval;
2535  bool evalFail= false;
2536  CFList list;
2537  bool GF= false;
2538  CanonicalForm LCA= LC (A);
2539  CanonicalForm tmp;
2540  CFArray gcds= CFArray (biggestSize);
2541  CFList * pEvalPoints= new CFList [biggestSize];
2542  Variable V_buf= alpha, V_buf4= alpha;
2543  CFList source, dest;
2544  CanonicalForm prim_elem, im_prim_elem;
2545  CanonicalForm prim_elem_alpha, im_prim_elem_alpha;
2546  for (int i= 0; i < biggestSize; i++)
2547  {
2548    if (i == 0)
2549    {
2550      if (getCharacteristic() > 3)
2551        evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf,
2552                                    evalFail, list);
2553      else
2554        evalFail= true;
2555
2556      if (evalFail)
2557      {
2558        if (V_buf.level() != 1)
2559        {
2560          do
2561          {
2562            Variable V_buf3= V_buf;
2563            V_buf= chooseExtension (V_buf);
2564            source= CFList();
2565            dest= CFList();
2566
2567            bool prim_fail= false;
2568            Variable V_buf2;
2569            prim_elem= primitiveElement (V_buf4, V_buf2, prim_fail);
2570            if (V_buf4 == alpha && alpha.level() != 1)
2571              prim_elem_alpha= prim_elem;
2572
2573            ASSERT (!prim_fail, "failure in integer factorizer");
2574            if (prim_fail)
2575              ; //ERROR
2576            else
2577              im_prim_elem= mapPrimElem (prim_elem, V_buf4, V_buf);
2578
2579            DEBOUTLN (cerr, "getMipo (V_buf)= " << getMipo (V_buf));
2580            DEBOUTLN (cerr, "getMipo (V_buf2)= " << getMipo (V_buf2));
2581
2582            if (V_buf4 == alpha && alpha.level() != 1)
2583              im_prim_elem_alpha= im_prim_elem;
2584            else if (alpha.level() != 1)
2585              im_prim_elem_alpha= mapUp (im_prim_elem_alpha, V_buf4, V_buf,
2586                                         prim_elem, im_prim_elem, source, dest);
2587
2588            for (CFListIterator i= list; i.hasItem(); i++)
2589              i.getItem()= mapUp (i.getItem(), V_buf4, V_buf, prim_elem,
2590                                im_prim_elem, source, dest);
2591            if (alpha.level() != 1)
2592            {
2593              A= mapUp (A, V_buf4, V_buf, prim_elem, im_prim_elem, source,dest);
2594              B= mapUp (B, V_buf4, V_buf, prim_elem, im_prim_elem, source,dest);
2595            }
2596            evalFail= false;
2597            V_buf4= V_buf;
2598            evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf,
2599                                          evalFail, list);
2600          } while (evalFail);
2601        }
2602        else
2603        {
2604          CanonicalForm mipo;
2605          int deg= 2;
2606          bool initialized= false;
2607          do
2608          {
2609            mipo= randomIrredpoly (deg, x);
2610            if (initialized)
2611              setMipo (V_buf, mipo);
2612            else
2613              V_buf= rootOf (mipo);
2614            evalFail= false;
2615            initialized= true;
2616            evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf,
2617                                          evalFail, list);
2618            deg++;
2619          } while (evalFail);
2620          V_buf4= V_buf;
2621        }
2622      }
2623    }
2624    else
2625    {
2626      mult (evalPoints, pEvalPoints[0]);
2627      eval (A, B, Aeval, Beval, evalPoints);
2628    }
2629
2630    g= gcd (Aeval, Beval);
2631    g /= Lc (g);
2632
2633    if (degree (g) != degree (skel) || (skelSize != size (g)))
2634    {
2635      delete[] pEvalPoints;
2636      fail= true;
2637      if (alpha.level() != 1 && V_buf != alpha)
2638        prune1 (alpha);
2639      return 0;
2640    }
2641    CFIterator l= skel;
2642    for (CFIterator k= g; k.hasTerms(); k++, l++)
2643    {
2644      if (k.exp() != l.exp())
2645      {
2646        delete[] pEvalPoints;
2647        fail= true;
2648        if (alpha.level() != 1 && V_buf != alpha)
2649          prune1 (alpha);
2650        return 0;
2651      }
2652    }
2653    pEvalPoints[i]= evalPoints;
2654    gcds[i]= g;
2655
2656    tmp= 0;
2657    int j= 0;
2658    for (CFListIterator k= evalPoints; k.hasItem(); k++, j++)
2659      tmp += k.getItem()*power (x, j);
2660    list.append (tmp);
2661  }
2662
2663  if (Monoms.size() == 0)
2664    Monoms= getMonoms (skel);
2665
2666  coeffMonoms= new CFArray [skelSize];
2667
2668  j= 0;
2669  for (CFIterator i= skel; i.hasTerms(); i++, j++)
2670    coeffMonoms[j]= getMonoms (i.coeff());
2671
2672  int minimalColumnsIndex;
2673  if (skelSize > 1)
2674    minimalColumnsIndex= 1;
2675  else
2676    minimalColumnsIndex= 0;
2677  int minimalColumns=-1;
2678
2679  CFArray* pM= new CFArray [skelSize];
2680  CFMatrix Mat;
2681  // find the Matrix with minimal number of columns
2682  for (int i= 0; i < skelSize; i++)
2683  {
2684    pM[i]= CFArray (coeffMonoms[i].size());
2685    if (i == 1)
2686      minimalColumns= coeffMonoms[i].size();
2687    if (i > 1)
2688    {
2689      if (minimalColumns > coeffMonoms[i].size())
2690      {
2691        minimalColumns= coeffMonoms[i].size();
2692        minimalColumnsIndex= i;
2693      }
2694    }
2695  }
2696  CFMatrix* pMat= new CFMatrix [2];
2697  pMat[0]= CFMatrix (biggestSize, coeffMonoms[0].size());
2698  pMat[1]= CFMatrix (biggestSize, coeffMonoms[minimalColumnsIndex].size());
2699  CFArray* pL= new CFArray [skelSize];
2700  for (int i= 0; i < biggestSize; i++)
2701  {
2702    CFIterator l= gcds [i];
2703    evalPoints= pEvalPoints [i];
2704    for (int k= 0; k < skelSize; k++, l++)
2705    {
2706      if (i == 0)
2707        pL[k]= CFArray (biggestSize);
2708      pL[k] [i]= l.coeff();
2709
2710      if (i == 0 && k != 0 && k != minimalColumnsIndex)
2711        pM[k]= evaluate (coeffMonoms[k], evalPoints);
2712      else if (i == 0 && (k == 0 || k == minimalColumnsIndex))
2713        coeffEval= evaluate (coeffMonoms[k], evalPoints);
2714      if (i > 0 && (k == 0 || k == minimalColumnsIndex))
2715        coeffEval= evaluate (coeffMonoms[k], evalPoints);
2716
2717      if (k == 0)
2718      {
2719        if (pMat[k].rows() >= i + 1)
2720        {
2721          for (int ind= 1; ind < coeffEval.size() + 1; ind++)
2722            pMat[k] (i + 1, ind)= coeffEval[ind - 1];
2723        }
2724      }
2725      if (k == minimalColumnsIndex)
2726      {
2727        if (pMat[1].rows() >= i + 1)
2728        {
2729          for (int ind= 1; ind < coeffEval.size() + 1; ind++)
2730            pMat[1] (i + 1, ind)= coeffEval[ind - 1];
2731        }
2732      }
2733    }
2734  }
2735
2736  CFArray solution;
2737  CanonicalForm result= 0;
2738  int ind= 1;
2739  int matRows, matColumns;
2740  matRows= pMat[1].rows();
2741  matColumns= pMat[0].columns() - 1;
2742  matColumns += pMat[1].columns();
2743
2744  Mat= CFMatrix (matRows, matColumns);
2745  for (int i= 1; i <= matRows; i++)
2746    for (int j= 1; j <= pMat[1].columns(); j++)
2747      Mat (i, j)= pMat[1] (i, j);
2748
2749  for (int j= pMat[1].columns() + 1; j <= matColumns;
2750       j++, ind++)
2751  {
2752    for (int i= 1; i <= matRows; i++)
2753      Mat (i, j)= (-pMat [0] (i, ind + 1))*pL[minimalColumnsIndex] [i - 1];
2754  }
2755
2756  CFArray firstColumn= CFArray (pMat[0].rows());
2757  for (int i= 0; i < pMat[0].rows(); i++)
2758    firstColumn [i]= pMat[0] (i + 1, 1);
2759
2760  CFArray bufArray= pL[minimalColumnsIndex];
2761
2762  for (int i= 0; i < biggestSize; i++)
2763    pL[minimalColumnsIndex] [i] *= firstColumn[i];
2764
2765  CFMatrix bufMat= pMat[1];
2766  pMat[1]= Mat;
2767
2768  if (V_buf.level() != 1)
2769    solution= solveSystemFq (pMat[1],
2770                             pL[minimalColumnsIndex], V_buf);
2771  else
2772    solution= solveSystemFp (pMat[1],
2773                             pL[minimalColumnsIndex]);
2774
2775  if (solution.size() == 0)
2776  { //Javadi's method failed, try deKleine, Monagan, Wittkopf instead
2777    CFMatrix bufMat0= pMat[0];
2778    delete [] pMat;
2779    pMat= new CFMatrix [skelSize];
2780    pL[minimalColumnsIndex]= bufArray;
2781    CFList* bufpEvalPoints= NULL;
2782    CFArray bufGcds;
2783    if (biggestSize != biggestSize2)
2784    {
2785      bufpEvalPoints= pEvalPoints;
2786      pEvalPoints= new CFList [biggestSize2];
2787      bufGcds= gcds;
2788      gcds= CFArray (biggestSize2);
2789      for (int i= 0; i < biggestSize; i++)
2790      {
2791        pEvalPoints[i]= bufpEvalPoints [i];
2792        gcds[i]= bufGcds[i];
2793      }
2794      for (int i= 0; i < biggestSize2 - biggestSize; i++)
2795      {
2796        mult (evalPoints, pEvalPoints[0]);
2797        eval (A, B, Aeval, Beval, evalPoints);
2798        g= gcd (Aeval, Beval);
2799        g /= Lc (g);
2800        if (degree (g) != degree (skel) || (skelSize != size (g)))
2801        {
2802          delete[] pEvalPoints;
2803          delete[] pMat;
2804          delete[] pL;
2805          delete[] coeffMonoms;
2806          delete[] pM;
2807          if (bufpEvalPoints != NULL)
2808            delete [] bufpEvalPoints;
2809          fail= true;
2810          if (alpha.level() != 1 && V_buf != alpha)
2811            prune1 (alpha);
2812          return 0;
2813        }
2814        CFIterator l= skel;
2815        for (CFIterator k= g; k.hasTerms(); k++, l++)
2816        {
2817          if (k.exp() != l.exp())
2818          {
2819            delete[] pEvalPoints;
2820            delete[] pMat;
2821            delete[] pL;
2822            delete[] coeffMonoms;
2823            delete[] pM;
2824            if (bufpEvalPoints != NULL)
2825              delete [] bufpEvalPoints;
2826            fail= true;
2827            if (alpha.level() != 1 && V_buf != alpha)
2828              prune1 (alpha);
2829            return 0;
2830          }
2831        }
2832        pEvalPoints[i + biggestSize]= evalPoints;
2833        gcds[i + biggestSize]= g;
2834      }
2835    }
2836    for (int i= 0; i < biggestSize; i++)
2837    {
2838      CFIterator l= gcds [i];
2839      evalPoints= pEvalPoints [i];
2840      for (int k= 1; k < skelSize; k++, l++)
2841      {
2842        if (i == 0)
2843          pMat[k]= CFMatrix (biggestSize2,coeffMonoms[k].size()+biggestSize2-1);
2844        if (k == minimalColumnsIndex)
2845          continue;
2846        coeffEval= evaluate (coeffMonoms[k], evalPoints);
2847        if (pMat[k].rows() >= i + 1)
2848        {
2849          for (int ind= 1; ind < coeffEval.size() + 1; ind++)
2850            pMat[k] (i + 1, ind)= coeffEval[ind - 1];
2851        }
2852      }
2853    }
2854    Mat= bufMat0;
2855    pMat[0]= CFMatrix (biggestSize2, coeffMonoms[0].size() + biggestSize2 - 1);
2856    for (int j= 1; j <= Mat.rows(); j++)
2857      for (int k= 1; k <= Mat.columns(); k++)
2858        pMat [0] (j,k)= Mat (j,k);
2859    Mat= bufMat;
2860    for (int j= 1; j <= Mat.rows(); j++)
2861      for (int k= 1; k <= Mat.columns(); k++)
2862        pMat [minimalColumnsIndex] (j,k)= Mat (j,k);
2863    // write old matrix entries into new matrices
2864    for (int i= 0; i < skelSize; i++)
2865    {
2866      bufArray= pL[i];
2867      pL[i]= CFArray (biggestSize2);
2868      for (int j= 0; j < bufArray.size(); j++)
2869        pL[i] [j]= bufArray [j];
2870    }
2871    //write old vector entries into new and add new entries to old matrices
2872    for (int i= 0; i < biggestSize2 - biggestSize; i++)
2873    {
2874      CFIterator l= gcds [i + biggestSize];
2875      evalPoints= pEvalPoints [i + biggestSize];
2876      for (int k= 0; k < skelSize; k++, l++)
2877      {
2878        pL[k] [i + biggestSize]= l.coeff();
2879        coeffEval= evaluate (coeffMonoms[k], evalPoints);
2880        if (pMat[k].rows() >= i + biggestSize + 1)
2881        {
2882          for (int ind= 1; ind < coeffEval.size() + 1; ind++)
2883            pMat[k] (i + biggestSize + 1, ind)= coeffEval[ind - 1];
2884        }
2885      }
2886    }
2887    // begin new
2888    for (int i= 0; i < skelSize; i++)
2889    {
2890      if (pL[i].size() > 1)
2891      {
2892        for (int j= 2; j <= pMat[i].rows(); j++)
2893          pMat[i] (j, coeffMonoms[i].size() + j - 1)=
2894              -pL[i] [j - 1];
2895      }
2896    }
2897
2898    matColumns= biggestSize2 - 1;
2899    matRows= 0;
2900    for (int i= 0; i < skelSize; i++)
2901    {
2902      if (V_buf.level() == 1)
2903        (void) gaussianElimFp (pMat[i], pL[i]);
2904      else
2905        (void) gaussianElimFq (pMat[i], pL[i], V_buf);
2906
2907      if (pMat[i] (coeffMonoms[i].size(), coeffMonoms[i].size()) == 0)
2908      {
2909        delete[] pEvalPoints;
2910        delete[] pMat;
2911        delete[] pL;
2912        delete[] coeffMonoms;
2913        delete[] pM;
2914        if (bufpEvalPoints != NULL)
2915          delete [] bufpEvalPoints;
2916        fail= true;
2917        if (alpha.level() != 1 && V_buf != alpha)
2918          prune1 (alpha);
2919        return 0;
2920      }
2921      matRows += pMat[i].rows() - coeffMonoms[i].size();
2922    }
2923
2924    CFMatrix bufMat;
2925    Mat= CFMatrix (matRows, matColumns);
2926    ind= 0;
2927    bufArray= CFArray (matRows);
2928    CFArray bufArray2;
2929    for (int i= 0; i < skelSize; i++)
2930    {
2931      if (coeffMonoms[i].size() + 1 >= pMat[i].rows() || coeffMonoms[i].size() + 1 >= pMat[i].columns())
2932      {
2933        delete[] pEvalPoints;
2934        delete[] pMat;
2935        delete[] pL;
2936        delete[] coeffMonoms;
2937        delete[] pM;
2938        if (bufpEvalPoints != NULL)
2939          delete [] bufpEvalPoints;
2940        fail= true;
2941        if (alpha.level() != 1 && V_buf != alpha)
2942          prune1 (alpha);
2943        return 0;
2944      }
2945      bufMat= pMat[i] (coeffMonoms[i].size() + 1, pMat[i].rows(),
2946                       coeffMonoms[i].size() + 1, pMat[i].columns());
2947
2948      for (int j= 1; j <= bufMat.rows(); j++)
2949        for (int k= 1; k <= bufMat.columns(); k++)
2950          Mat (j + ind, k)= bufMat(j, k);
2951      bufArray2= coeffMonoms[i].size();
2952      for (int j= 1; j <= pMat[i].rows(); j++)
2953      {
2954        if (j > coeffMonoms[i].size())
2955          bufArray [j-coeffMonoms[i].size() + ind - 1]= pL[i] [j - 1];
2956        else
2957          bufArray2 [j - 1]= pL[i] [j - 1];
2958      }
2959      pL[i]= bufArray2;
2960      ind += bufMat.rows();
2961      pMat[i]= pMat[i] (1, coeffMonoms[i].size(), 1, pMat[i].columns());
2962    }
2963
2964    if (V_buf.level() != 1)
2965      solution= solveSystemFq (Mat, bufArray, V_buf);
2966    else
2967      solution= solveSystemFp (Mat, bufArray);
2968
2969    if (solution.size() == 0)
2970    {
2971      delete[] pEvalPoints;
2972      delete[] pMat;
2973      delete[] pL;
2974      delete[] coeffMonoms;
2975      delete[] pM;
2976      if (bufpEvalPoints != NULL)
2977        delete [] bufpEvalPoints;
2978      fail= true;
2979      if (alpha.level() != 1 && V_buf != alpha)
2980        prune1 (alpha);
2981      return 0;
2982    }
2983
2984    ind= 0;
2985    result= 0;
2986    CFArray bufSolution;
2987    for (int i= 0; i < skelSize; i++)
2988    {
2989      bufSolution= readOffSolution (pMat[i], pL[i], solution);
2990      for (int i= 0; i < bufSolution.size(); i++, ind++)
2991        result += Monoms [ind]*bufSolution[i];
2992    }
2993    if (alpha.level() != 1 && V_buf != alpha)
2994    {
2995      CFList u, v;
2996      result= mapDown (result,prim_elem_alpha, im_prim_elem_alpha, alpha, u, v);
2997      prune1 (alpha);
2998    }
2999    result= N(result);
3000    delete[] pEvalPoints;
3001    delete[] pMat;
3002    delete[] pL;
3003    delete[] coeffMonoms;
3004    delete[] pM;
3005
3006    if (bufpEvalPoints != NULL)
3007      delete [] bufpEvalPoints;
3008    if (fdivides (result, F) && fdivides (result, G))
3009      return result;
3010    else
3011    {
3012      fail= true;
3013      return 0;
3014    }
3015  } // end of deKleine, Monagan & Wittkopf
3016
3017  result += Monoms[0];
3018  int ind2= 0, ind3= 2;
3019  ind= 0;
3020  for (int l= 0; l < minimalColumnsIndex; l++)
3021    ind += coeffMonoms[l].size();
3022  for (int l= solution.size() - pMat[0].columns() + 1; l < solution.size();
3023       l++, ind2++, ind3++)
3024  {
3025    result += solution[l]*Monoms [1 + ind2];
3026    for (int i= 0; i < pMat[0].rows(); i++)
3027      firstColumn[i] += solution[l]*pMat[0] (i + 1, ind3);
3028  }
3029  for (int l= 0; l < solution.size() + 1 - pMat[0].columns(); l++)
3030    result += solution[l]*Monoms [ind + l];
3031  ind= coeffMonoms[0].size();
3032  for (int k= 1; k < skelSize; k++)
3033  {
3034    if (k == minimalColumnsIndex)
3035    {
3036      ind += coeffMonoms[k].size();
3037      continue;
3038    }
3039    if (k != minimalColumnsIndex)
3040    {
3041      for (int i= 0; i < biggestSize; i++)
3042        pL[k] [i] *= firstColumn [i];
3043    }
3044
3045    if (biggestSize != coeffMonoms[k].size() && k != minimalColumnsIndex)
3046    {
3047      bufArray= CFArray (coeffMonoms[k].size());
3048      for (int i= 0; i < bufArray.size(); i++)
3049        bufArray [i]= pL[k] [i];
3050      solution= solveGeneralVandermonde (pM[k], bufArray);
3051    }
3052    else
3053      solution= solveGeneralVandermonde (pM[k], pL[k]);
3054
3055    if (solution.size() == 0)
3056    {
3057      delete[] pEvalPoints;
3058      delete[] pMat;
3059      delete[] pL;
3060      delete[] coeffMonoms;
3061      delete[] pM;
3062      fail= true;
3063      if (alpha.level() != 1 && V_buf != alpha)
3064        prune1 (alpha);
3065      return 0;
3066    }
3067    if (k != minimalColumnsIndex)
3068    {
3069      for (int l= 0; l < solution.size(); l++)
3070        result += solution[l]*Monoms [ind + l];
3071      ind += solution.size();
3072    }
3073  }
3074
3075  delete[] pEvalPoints;
3076  delete[] pMat;
3077  delete[] pL;
3078  delete[] pM;
3079  delete[] coeffMonoms;
3080
3081  if (alpha.level() != 1 && V_buf != alpha)
3082  {
3083    CFList u, v;
3084    result= mapDown (result, prim_elem, im_prim_elem, alpha, u, v);
3085    prune1 (alpha);
3086  }
3087  result= N(result);
3088
3089  if (fdivides (result, F) && fdivides (result, G))
3090    return result;
3091  else
3092  {
3093    fail= true;
3094    return 0;
3095  }
3096}
3097#endif
3098
3099#ifdef HAVE_NTL // primitiveElement,FindRoot
3100CanonicalForm sparseGCDFq (const CanonicalForm& F, const CanonicalForm& G,
3101                           const Variable & alpha, CFList& l, bool& topLevel)
3102{
3103  CanonicalForm A= F;
3104  CanonicalForm B= G;
3105  if (F.isZero() && degree(G) > 0) return G/Lc(G);
3106  else if (G.isZero() && degree (F) > 0) return F/Lc(F);
3107  else if (F.isZero() && G.isZero()) return F.genOne();
3108  if (F.inBaseDomain() || G.inBaseDomain()) return F.genOne();
3109  if (F.isUnivariate() && fdivides(F, G)) return F/Lc(F);
3110  if (G.isUnivariate() && fdivides(G, F)) return G/Lc(G);
3111  if (F == G) return F/Lc(F);
3112
3113  CFMap M,N;
3114  int best_level= myCompress (A, B, M, N, topLevel);
3115
3116  if (best_level == 0) return B.genOne();
3117
3118  A= M(A);
3119  B= M(B);
3120
3121  Variable x= Variable (1);
3122
3123  //univariate case
3124  if (A.isUnivariate() && B.isUnivariate())
3125    return N (gcd (A, B));
3126
3127  CanonicalForm cA, cB;    // content of A and B
3128  CanonicalForm ppA, ppB;    // primitive part of A and B
3129  CanonicalForm gcdcAcB;
3130
3131  cA = uni_content (A);
3132  cB = uni_content (B);
3133  gcdcAcB= gcd (cA, cB);
3134  ppA= A/cA;
3135  ppB= B/cB;
3136
3137  CanonicalForm lcA, lcB;  // leading coefficients of A and B
3138  CanonicalForm gcdlcAlcB;
3139  lcA= uni_lcoeff (ppA);
3140  lcB= uni_lcoeff (ppB);
3141
3142  if (fdivides (lcA, lcB))
3143  {
3144    if (fdivides (A, B))
3145      return F/Lc(F);
3146  }
3147  if (fdivides (lcB, lcA))
3148  {
3149    if (fdivides (B, A))
3150      return G/Lc(G);
3151  }
3152
3153  gcdlcAlcB= gcd (lcA, lcB);
3154
3155  int d= totaldegree (ppA, Variable (2), Variable (ppA.level()));
3156  int d0;
3157
3158  if (d == 0)
3159    return N(gcdcAcB);
3160  d0= totaldegree (ppB, Variable (2), Variable (ppB.level()));
3161
3162  if (d0 < d)
3163    d= d0;
3164
3165  if (d == 0)
3166    return N(gcdcAcB);
3167
3168  CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH, skeleton;
3169  CanonicalForm newtonPoly= 1;
3170  m= gcdlcAlcB;
3171  G_m= 0;
3172  H= 0;
3173  bool fail= false;
3174  topLevel= false;
3175  bool inextension= false;
3176  Variable V_buf= alpha, V_buf4= alpha;
3177  CanonicalForm prim_elem, im_prim_elem;
3178  CanonicalForm prim_elem_alpha, im_prim_elem_alpha;
3179  CFList source, dest;
3180  do // first do
3181  {
3182    random_element= randomElement (m, V_buf, l, fail);
3183    if (random_element == 0 && !fail)
3184    {
3185      if (!find (l, random_element))
3186        l.append (random_element);
3187      continue;
3188    }
3189    if (fail)
3190    {
3191      source= CFList();
3192      dest= CFList();
3193
3194      Variable V_buf3= V_buf;
3195      V_buf= chooseExtension (V_buf);
3196      bool prim_fail= false;
3197      Variable V_buf2;
3198      prim_elem= primitiveElement (V_buf4, V_buf2, prim_fail);
3199      if (V_buf4 == alpha)
3200        prim_elem_alpha= prim_elem;
3201
3202      if (V_buf3 != V_buf4)
3203      {
3204        m= mapDown (m, prim_elem, im_prim_elem, V_buf4, source, dest);
3205        G_m= mapDown (m, prim_elem, im_prim_elem, V_buf4, source, dest);
3206        newtonPoly= mapDown (newtonPoly, prim_elem, im_prim_elem, V_buf4,
3207                             source, dest);
3208        ppA= mapDown (ppA, prim_elem, im_prim_elem, V_buf4, source, dest);
3209        ppB= mapDown (ppB, prim_elem, im_prim_elem, V_buf4, source, dest);
3210        gcdlcAlcB= mapDown (gcdlcAlcB, prim_elem, im_prim_elem, V_buf4, source,
3211                            dest);
3212        for (CFListIterator i= l; i.hasItem(); i++)
3213          i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, V_buf4,
3214                                source, dest);
3215      }
3216
3217      ASSERT (!prim_fail, "failure in integer factorizer");
3218      if (prim_fail)
3219        ; //ERROR
3220      else
3221        im_prim_elem= mapPrimElem (prim_elem, V_buf4, V_buf);
3222
3223      if (V_buf4 == alpha)
3224        im_prim_elem_alpha= im_prim_elem;
3225      else
3226        im_prim_elem_alpha= mapUp (im_prim_elem_alpha, V_buf4, V_buf, prim_elem,
3227                                   im_prim_elem, source, dest);
3228
3229      DEBOUTLN (cerr, "getMipo (V_buf4)= " << getMipo (V_buf4));
3230      DEBOUTLN (cerr, "getMipo (V_buf2)= " << getMipo (V_buf2));
3231      inextension= true;
3232      for (CFListIterator i= l; i.hasItem(); i++)
3233        i.getItem()= mapUp (i.getItem(), V_buf4, V_buf, prim_elem,
3234                             im_prim_elem, source, dest);
3235      m= mapUp (m, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
3236      G_m= mapUp (G_m, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
3237      newtonPoly= mapUp (newtonPoly, V_buf4, V_buf, prim_elem, im_prim_elem,
3238                          source, dest);
3239      ppA= mapUp (ppA, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
3240      ppB= mapUp (ppB, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
3241      gcdlcAlcB= mapUp (gcdlcAlcB, V_buf4, V_buf, prim_elem, im_prim_elem,
3242                         source, dest);
3243
3244      fail= false;
3245      random_element= randomElement (m, V_buf, l, fail );
3246      DEBOUTLN (cerr, "fail= " << fail);
3247      CFList list;
3248      TIMING_START (gcd_recursion);
3249      G_random_element=
3250      sparseGCDFq (ppA (random_element, x), ppB (random_element, x), V_buf,
3251                        list, topLevel);
3252      TIMING_END_AND_PRINT (gcd_recursion,
3253                            "time for recursive call: ");
3254      DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3255      V_buf4= V_buf;
3256    }
3257    else
3258    {
3259      CFList list;
3260      TIMING_START (gcd_recursion);
3261      G_random_element=
3262      sparseGCDFq (ppA(random_element, x), ppB(random_element, x), V_buf,
3263                        list, topLevel);
3264      TIMING_END_AND_PRINT (gcd_recursion,
3265                            "time for recursive call: ");
3266      DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3267    }
3268
3269    if (!G_random_element.inCoeffDomain())
3270      d0= totaldegree (G_random_element, Variable(2),
3271                       Variable (G_random_element.level()));
3272    else
3273      d0= 0;
3274
3275    if (d0 == 0)
3276    {
3277      if (inextension)
3278        prune1 (alpha);
3279      return N(gcdcAcB);
3280    }
3281    if (d0 >  d)
3282    {
3283      if (!find (l, random_element))
3284        l.append (random_element);
3285      continue;
3286    }
3287
3288    G_random_element=
3289    (gcdlcAlcB(random_element, x)/uni_lcoeff (G_random_element))
3290    * G_random_element;
3291
3292    skeleton= G_random_element;
3293    if (!G_random_element.inCoeffDomain())
3294      d0= totaldegree (G_random_element, Variable(2),
3295                       Variable (G_random_element.level()));
3296    else
3297      d0= 0;
3298
3299    if (d0 <  d)
3300    {
3301      m= gcdlcAlcB;
3302      newtonPoly= 1;
3303      G_m= 0;
3304      d= d0;
3305    }
3306
3307    H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
3308    if (uni_lcoeff (H) == gcdlcAlcB)
3309    {
3310      cH= uni_content (H);
3311      ppH= H/cH;
3312      if (inextension)
3313      {
3314        CFList u, v;
3315        //maybe it's better to test if ppH is an element of F(\alpha) before
3316        //mapping down
3317        if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
3318        {
3319          DEBOUTLN (cerr, "ppH before mapDown= " << ppH);
3320          ppH= mapDown (ppH, prim_elem_alpha, im_prim_elem_alpha, alpha, u, v);
3321          ppH /= Lc(ppH);
3322          DEBOUTLN (cerr, "ppH after mapDown= " << ppH);
3323          prune1 (alpha);
3324          return N(gcdcAcB*ppH);
3325        }
3326      }
3327      else if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
3328        return N(gcdcAcB*ppH);
3329    }
3330    G_m= H;
3331    newtonPoly= newtonPoly*(x - random_element);
3332    m= m*(x - random_element);
3333    if (!find (l, random_element))
3334      l.append (random_element);
3335
3336    if (getCharacteristic () > 3 && size (skeleton) < 100)
3337    {
3338      CFArray Monoms;
3339      CFArray *coeffMonoms;
3340      do //second do
3341      {
3342        random_element= randomElement (m, V_buf, l, fail);
3343        if (random_element == 0 && !fail)
3344        {
3345          if (!find (l, random_element))
3346            l.append (random_element);
3347          continue;
3348        }
3349        if (fail)
3350        {
3351          source= CFList();
3352          dest= CFList();
3353
3354          Variable V_buf3= V_buf;
3355          V_buf= chooseExtension (V_buf);
3356          bool prim_fail= false;
3357          Variable V_buf2;
3358          prim_elem= primitiveElement (V_buf4, V_buf2, prim_fail);
3359          if (V_buf4 == alpha)
3360            prim_elem_alpha= prim_elem;
3361
3362          if (V_buf3 != V_buf4)
3363          {
3364            m= mapDown (m, prim_elem, im_prim_elem, V_buf4, source, dest);
3365            G_m= mapDown (m, prim_elem, im_prim_elem, V_buf4, source, dest);
3366            newtonPoly= mapDown (newtonPoly, prim_elem, im_prim_elem, V_buf4,
3367                                 source, dest);
3368            ppA= mapDown (ppA, prim_elem, im_prim_elem, V_buf4, source, dest);
3369            ppB= mapDown (ppB, prim_elem, im_prim_elem, V_buf4, source, dest);
3370            gcdlcAlcB= mapDown (gcdlcAlcB, prim_elem, im_prim_elem, V_buf4,
3371                                source, dest);
3372            for (CFListIterator i= l; i.hasItem(); i++)
3373              i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, V_buf4,
3374                                    source, dest);
3375          }
3376
3377          ASSERT (!prim_fail, "failure in integer factorizer");
3378          if (prim_fail)
3379            ; //ERROR
3380          else
3381            im_prim_elem= mapPrimElem (prim_elem, V_buf4, V_buf);
3382
3383          if (V_buf4 == alpha)
3384            im_prim_elem_alpha= im_prim_elem;
3385          else
3386            im_prim_elem_alpha= mapUp (im_prim_elem_alpha, V_buf4, V_buf,
3387                                       prim_elem, im_prim_elem, source, dest);
3388
3389          DEBOUTLN (cerr, "getMipo (V_buf4)= " << getMipo (V_buf4));
3390          DEBOUTLN (cerr, "getMipo (V_buf2)= " << getMipo (V_buf2));
3391          inextension= true;
3392          for (CFListIterator i= l; i.hasItem(); i++)
3393            i.getItem()= mapUp (i.getItem(), V_buf4, V_buf, prim_elem,
3394                                im_prim_elem, source, dest);
3395          m= mapUp (m, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
3396          G_m= mapUp (G_m, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
3397          newtonPoly= mapUp (newtonPoly, V_buf4, V_buf, prim_elem, im_prim_elem,
3398                              source, dest);
3399          ppA= mapUp (ppA, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
3400          ppB= mapUp (ppB, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
3401
3402          gcdlcAlcB= mapUp (gcdlcAlcB, V_buf4, V_buf, prim_elem, im_prim_elem,
3403                            source, dest);
3404
3405          fail= false;
3406          random_element= randomElement (m, V_buf, l, fail);
3407          DEBOUTLN (cerr, "fail= " << fail);
3408          CFList list;
3409          TIMING_START (gcd_recursion);
3410
3411          V_buf4= V_buf;
3412
3413          //sparseInterpolation
3414          bool sparseFail= false;
3415          if (LC (skeleton).inCoeffDomain())
3416            G_random_element=
3417            monicSparseInterpol (ppA (random_element, x), ppB(random_element,x),
3418                                skeleton,V_buf, sparseFail, coeffMonoms,Monoms);
3419          else
3420            G_random_element=
3421            nonMonicSparseInterpol (ppA(random_element,x),ppB(random_element,x),
3422                                    skeleton, V_buf, sparseFail, coeffMonoms,
3423                                    Monoms);
3424          if (sparseFail)
3425            break;
3426
3427          TIMING_END_AND_PRINT (gcd_recursion,
3428                                "time for recursive call: ");
3429          DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3430        }
3431        else
3432        {
3433          CFList list;
3434          TIMING_START (gcd_recursion);
3435          bool sparseFail= false;
3436          if (LC (skeleton).inCoeffDomain())
3437            G_random_element=
3438            monicSparseInterpol (ppA (random_element, x),ppB(random_element, x),
3439                                skeleton,V_buf, sparseFail,coeffMonoms, Monoms);
3440          else
3441            G_random_element=
3442            nonMonicSparseInterpol (ppA(random_element,x),ppB(random_element,x),
3443                                    skeleton, V_buf, sparseFail, coeffMonoms,
3444                                    Monoms);
3445          if (sparseFail)
3446            break;
3447
3448          TIMING_END_AND_PRINT (gcd_recursion,
3449                                "time for recursive call: ");
3450          DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3451        }
3452
3453        if (!G_random_element.inCoeffDomain())
3454          d0= totaldegree (G_random_element, Variable(2),
3455                           Variable (G_random_element.level()));
3456        else
3457          d0= 0;
3458
3459        if (d0 == 0)
3460        {
3461          if (inextension)
3462            prune1 (alpha);
3463          return N(gcdcAcB);
3464        }
3465        if (d0 >  d)
3466        {
3467          if (!find (l, random_element))
3468            l.append (random_element);
3469          continue;
3470        }
3471
3472        G_random_element=
3473        (gcdlcAlcB(random_element, x)/uni_lcoeff (G_random_element))
3474        * G_random_element;
3475
3476        if (!G_random_element.inCoeffDomain())
3477          d0= totaldegree (G_random_element, Variable(2),
3478                          Variable (G_random_element.level()));
3479        else
3480          d0= 0;
3481
3482        if (d0 <  d)
3483        {
3484          m= gcdlcAlcB;
3485          newtonPoly= 1;
3486          G_m= 0;
3487          d= d0;
3488        }
3489
3490        TIMING_START (newton_interpolation);
3491        H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
3492        TIMING_END_AND_PRINT (newton_interpolation,
3493                              "time for newton interpolation: ");
3494
3495        //termination test
3496        if (uni_lcoeff (H) == gcdlcAlcB)
3497        {
3498          cH= uni_content (H);
3499          ppH= H/cH;
3500          if (inextension)
3501          {
3502            CFList u, v;
3503            //maybe it's better to test if ppH is an element of F(\alpha) before
3504            //mapping down
3505            if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
3506            {
3507              DEBOUTLN (cerr, "ppH before mapDown= " << ppH);
3508              ppH= mapDown (ppH, prim_elem_alpha, im_prim_elem_alpha, alpha, u, v);
3509              ppH /= Lc(ppH);
3510              DEBOUTLN (cerr, "ppH after mapDown= " << ppH);
3511              prune1 (alpha);
3512              return N(gcdcAcB*ppH);
3513            }
3514          }
3515          else if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
3516          {
3517            return N(gcdcAcB*ppH);
3518          }
3519        }
3520
3521        G_m= H;
3522        newtonPoly= newtonPoly*(x - random_element);
3523        m= m*(x - random_element);
3524        if (!find (l, random_element))
3525          l.append (random_element);
3526
3527      } while (1);
3528    }
3529  } while (1);
3530}
3531#endif
3532
3533#ifdef HAVE_NTL // primitiveElement,FindRoot
3534CanonicalForm sparseGCDFp (const CanonicalForm& F, const CanonicalForm& G,
3535                           bool& topLevel, CFList& l)
3536{
3537  CanonicalForm A= F;
3538  CanonicalForm B= G;
3539  if (F.isZero() && degree(G) > 0) return G/Lc(G);
3540  else if (G.isZero() && degree (F) > 0) return F/Lc(F);
3541  else if (F.isZero() && G.isZero()) return F.genOne();
3542  if (F.inBaseDomain() || G.inBaseDomain()) return F.genOne();
3543  if (F.isUnivariate() && fdivides(F, G)) return F/Lc(F);
3544  if (G.isUnivariate() && fdivides(G, F)) return G/Lc(G);
3545  if (F == G) return F/Lc(F);
3546
3547  CFMap M,N;
3548  int best_level= myCompress (A, B, M, N, topLevel);
3549
3550  if (best_level == 0) return B.genOne();
3551
3552  A= M(A);
3553  B= M(B);
3554
3555  Variable x= Variable (1);
3556
3557  //univariate case
3558  if (A.isUnivariate() && B.isUnivariate())
3559    return N (gcd (A, B));
3560
3561  CanonicalForm cA, cB;    // content of A and B
3562  CanonicalForm ppA, ppB;    // primitive part of A and B
3563  CanonicalForm gcdcAcB;
3564
3565  cA = uni_content (A);
3566  cB = uni_content (B);
3567  gcdcAcB= gcd (cA, cB);
3568  ppA= A/cA;
3569  ppB= B/cB;
3570
3571  CanonicalForm lcA, lcB;  // leading coefficients of A and B
3572  CanonicalForm gcdlcAlcB;
3573  lcA= uni_lcoeff (ppA);
3574  lcB= uni_lcoeff (ppB);
3575
3576  if (fdivides (lcA, lcB))
3577  {
3578    if (fdivides (A, B))
3579      return F/Lc(F);
3580  }
3581  if (fdivides (lcB, lcA))
3582  {
3583    if (fdivides (B, A))
3584      return G/Lc(G);
3585  }
3586
3587  gcdlcAlcB= gcd (lcA, lcB);
3588
3589  int d= totaldegree (ppA, Variable (2), Variable (ppA.level()));
3590  int d0;
3591
3592  if (d == 0)
3593    return N(gcdcAcB);
3594
3595  d0= totaldegree (ppB, Variable (2), Variable (ppB.level()));
3596
3597  if (d0 < d)
3598    d= d0;
3599
3600  if (d == 0)
3601    return N(gcdcAcB);
3602
3603  CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH, skeleton;
3604  CanonicalForm newtonPoly= 1;
3605  m= gcdlcAlcB;
3606  G_m= 0;
3607  H= 0;
3608  bool fail= false;
3609  topLevel= false;
3610  bool inextension= false;
3611  Variable V_buf, alpha, cleanUp;
3612  CanonicalForm prim_elem, im_prim_elem;
3613  CFList source, dest;
3614  do //first do
3615  {
3616    if (inextension)
3617      random_element= randomElement (m, V_buf, l, fail);
3618    else
3619      random_element= FpRandomElement (m, l, fail);
3620    if (random_element == 0 && !fail)
3621    {
3622      if (!find (l, random_element))
3623        l.append (random_element);
3624      continue;
3625    }
3626
3627    if (!fail && !inextension)
3628    {
3629      CFList list;
3630      TIMING_START (gcd_recursion);
3631      G_random_element=
3632      sparseGCDFp (ppA (random_element,x), ppB (random_element,x), topLevel,
3633                   list);
3634      TIMING_END_AND_PRINT (gcd_recursion,
3635                            "time for recursive call: ");
3636      DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3637    }
3638    else if (!fail && inextension)
3639    {
3640      CFList list;
3641      TIMING_START (gcd_recursion);
3642      G_random_element=
3643      sparseGCDFq (ppA (random_element, x), ppB (random_element, x), alpha,
3644                        list, topLevel);
3645      TIMING_END_AND_PRINT (gcd_recursion,
3646                            "time for recursive call: ");
3647      DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3648    }
3649    else if (fail && !inextension)
3650    {
3651      source= CFList();
3652      dest= CFList();
3653      CFList list;
3654      CanonicalForm mipo;
3655      int deg= 2;
3656      bool initialized= false;
3657      do
3658      {
3659        mipo= randomIrredpoly (deg, x);
3660        if (initialized)
3661          setMipo (alpha, mipo);
3662        else
3663          alpha= rootOf (mipo);
3664        inextension= true;
3665        fail= false;
3666        initialized= true;
3667        random_element= randomElement (m, alpha, l, fail);
3668        deg++;
3669      } while (fail);
3670      cleanUp= alpha;
3671      V_buf= alpha;
3672      list= CFList();
3673      TIMING_START (gcd_recursion);
3674      G_random_element=
3675      sparseGCDFq (ppA (random_element, x), ppB (random_element, x), alpha,
3676                        list, topLevel);
3677      TIMING_END_AND_PRINT (gcd_recursion,
3678                            "time for recursive call: ");
3679      DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3680    }
3681    else if (fail && inextension)
3682    {
3683      source= CFList();
3684      dest= CFList();
3685
3686      Variable V_buf3= V_buf;
3687      V_buf= chooseExtension (V_buf);
3688      bool prim_fail= false;
3689      Variable V_buf2;
3690      CanonicalForm prim_elem, im_prim_elem;
3691      prim_elem= primitiveElement (alpha, V_buf2, prim_fail);
3692
3693      if (V_buf3 != alpha)
3694      {
3695        m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
3696        G_m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
3697        newtonPoly= mapDown (newtonPoly, prim_elem, im_prim_elem, alpha, source,
3698                             dest);
3699        ppA= mapDown (ppA, prim_elem, im_prim_elem, alpha, source, dest);
3700        ppB= mapDown (ppB, prim_elem, im_prim_elem, alpha, source, dest);
3701        gcdlcAlcB= mapDown (gcdlcAlcB, prim_elem, im_prim_elem, alpha, source,
3702                            dest);
3703        for (CFListIterator i= l; i.hasItem(); i++)
3704          i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, alpha,
3705                                source, dest);
3706      }
3707
3708      ASSERT (!prim_fail, "failure in integer factorizer");
3709      if (prim_fail)
3710        ; //ERROR
3711      else
3712        im_prim_elem= mapPrimElem (prim_elem, alpha, V_buf);
3713
3714      DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (alpha));
3715      DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf2));
3716
3717      for (CFListIterator i= l; i.hasItem(); i++)
3718        i.getItem()= mapUp (i.getItem(), alpha, V_buf, prim_elem,
3719                             im_prim_elem, source, dest);
3720      m= mapUp (m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3721      G_m= mapUp (G_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3722      newtonPoly= mapUp (newtonPoly, alpha, V_buf, prim_elem, im_prim_elem,
3723                          source, dest);
3724      ppA= mapUp (ppA, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3725      ppB= mapUp (ppB, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3726      gcdlcAlcB= mapUp (gcdlcAlcB, alpha, V_buf, prim_elem, im_prim_elem,
3727                         source, dest);
3728      fail= false;
3729      random_element= randomElement (m, V_buf, l, fail );
3730      DEBOUTLN (cerr, "fail= " << fail);
3731      CFList list;
3732      TIMING_START (gcd_recursion);
3733      G_random_element=
3734      sparseGCDFq (ppA (random_element, x), ppB (random_element, x), V_buf,
3735                        list, topLevel);
3736      TIMING_END_AND_PRINT (gcd_recursion,
3737                            "time for recursive call: ");
3738      DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3739      alpha= V_buf;
3740    }
3741
3742    if (!G_random_element.inCoeffDomain())
3743      d0= totaldegree (G_random_element, Variable(2),
3744                       Variable (G_random_element.level()));
3745    else
3746      d0= 0;
3747
3748    if (d0 == 0)
3749    {
3750      if (inextension)
3751        prune (cleanUp);
3752      return N(gcdcAcB);
3753    }
3754    if (d0 >  d)
3755    {
3756      if (!find (l, random_element))
3757        l.append (random_element);
3758      continue;
3759    }
3760
3761    G_random_element=
3762    (gcdlcAlcB(random_element, x)/uni_lcoeff (G_random_element))
3763    * G_random_element;
3764
3765    skeleton= G_random_element;
3766
3767    if (!G_random_element.inCoeffDomain())
3768      d0= totaldegree (G_random_element, Variable(2),
3769                       Variable (G_random_element.level()));
3770    else
3771      d0= 0;
3772
3773    if (d0 <  d)
3774    {
3775      m= gcdlcAlcB;
3776      newtonPoly= 1;
3777      G_m= 0;
3778      d= d0;
3779    }
3780
3781    H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
3782
3783    if (uni_lcoeff (H) == gcdlcAlcB)
3784    {
3785      cH= uni_content (H);
3786      ppH= H/cH;
3787      ppH /= Lc (ppH);
3788      DEBOUTLN (cerr, "ppH= " << ppH);
3789
3790      if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
3791      {
3792        if (inextension)
3793          prune (cleanUp);
3794        return N(gcdcAcB*ppH);
3795      }
3796    }
3797    G_m= H;
3798    newtonPoly= newtonPoly*(x - random_element);
3799    m= m*(x - random_element);
3800    if (!find (l, random_element))
3801      l.append (random_element);
3802
3803    if ((getCharacteristic() > 3 && size (skeleton) < 200))
3804    {
3805      CFArray Monoms;
3806      CFArray* coeffMonoms;
3807
3808      do //second do
3809      {
3810        if (inextension)
3811          random_element= randomElement (m, alpha, l, fail);
3812        else
3813          random_element= FpRandomElement (m, l, fail);
3814        if (random_element == 0 && !fail)
3815        {
3816          if (!find (l, random_element))
3817            l.append (random_element);
3818          continue;
3819        }
3820
3821        bool sparseFail= false;
3822        if (!fail && !inextension)
3823        {
3824          CFList list;
3825          TIMING_START (gcd_recursion);
3826          if (LC (skeleton).inCoeffDomain())
3827            G_random_element=
3828            monicSparseInterpol(ppA(random_element, x), ppB (random_element, x),
3829                                skeleton, x, sparseFail, coeffMonoms,
3830                                Monoms);
3831          else
3832            G_random_element=
3833            nonMonicSparseInterpol(ppA(random_element,x), ppB(random_element,x),
3834                                    skeleton, x, sparseFail,
3835                                    coeffMonoms, Monoms);
3836          TIMING_END_AND_PRINT (gcd_recursion,
3837                                "time for recursive call: ");
3838          DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3839        }
3840        else if (!fail && inextension)
3841        {
3842          CFList list;
3843          TIMING_START (gcd_recursion);
3844          if (LC (skeleton).inCoeffDomain())
3845            G_random_element=
3846            monicSparseInterpol(ppA (random_element,x), ppB (random_element, x),
3847                                skeleton, alpha, sparseFail, coeffMonoms,
3848                                Monoms);
3849          else
3850            G_random_element=
3851            nonMonicSparseInterpol(ppA(random_element,x), ppB(random_element,x),
3852                                   skeleton, alpha, sparseFail, coeffMonoms,
3853                                   Monoms);
3854          TIMING_END_AND_PRINT (gcd_recursion,
3855                                "time for recursive call: ");
3856          DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3857        }
3858        else if (fail && !inextension)
3859        {
3860          source= CFList();
3861          dest= CFList();
3862          CFList list;
3863          CanonicalForm mipo;
3864          int deg= 2;
3865          bool initialized= false;
3866          do
3867          {
3868            mipo= randomIrredpoly (deg, x);
3869            if (initialized)
3870              setMipo (alpha, mipo);
3871            else
3872              alpha= rootOf (mipo);
3873            inextension= true;
3874            fail= false;
3875            initialized= true;
3876            random_element= randomElement (m, alpha, l, fail);
3877            deg++;
3878          } while (fail);
3879          cleanUp= alpha;
3880          V_buf= alpha;
3881          list= CFList();
3882          TIMING_START (gcd_recursion);
3883          if (LC (skeleton).inCoeffDomain())
3884            G_random_element=
3885            monicSparseInterpol (ppA (random_element,x), ppB (random_element,x),
3886                                skeleton, alpha, sparseFail, coeffMonoms,
3887                                Monoms);
3888          else
3889            G_random_element=
3890            nonMonicSparseInterpol(ppA(random_element,x), ppB(random_element,x),
3891                                   skeleton, alpha, sparseFail, coeffMonoms,
3892                                   Monoms);
3893          TIMING_END_AND_PRINT (gcd_recursion,
3894                                "time for recursive call: ");
3895          DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3896        }
3897        else if (fail && inextension)
3898        {
3899          source= CFList();
3900          dest= CFList();
3901
3902          Variable V_buf3= V_buf;
3903          V_buf= chooseExtension (V_buf);
3904          bool prim_fail= false;
3905          Variable V_buf2;
3906          CanonicalForm prim_elem, im_prim_elem;
3907          prim_elem= primitiveElement (alpha, V_buf2, prim_fail);
3908
3909          if (V_buf3 != alpha)
3910          {
3911            m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
3912            G_m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
3913            newtonPoly= mapDown (newtonPoly, prim_elem, im_prim_elem, alpha,
3914                                 source, dest);
3915            ppA= mapDown (ppA, prim_elem, im_prim_elem, alpha, source, dest);
3916            ppB= mapDown (ppB, prim_elem, im_prim_elem, alpha, source, dest);
3917            gcdlcAlcB= mapDown (gcdlcAlcB, prim_elem, im_prim_elem, alpha,
3918                                source, dest);
3919            for (CFListIterator i= l; i.hasItem(); i++)
3920              i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, alpha,
3921                                    source, dest);
3922          }
3923
3924          ASSERT (!prim_fail, "failure in integer factorizer");
3925          if (prim_fail)
3926            ; //ERROR
3927          else
3928            im_prim_elem= mapPrimElem (prim_elem, alpha, V_buf);
3929
3930          DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (alpha));
3931          DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf2));
3932
3933          for (CFListIterator i= l; i.hasItem(); i++)
3934            i.getItem()= mapUp (i.getItem(), alpha, V_buf, prim_elem,
3935                                im_prim_elem, source, dest);
3936          m= mapUp (m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3937          G_m= mapUp (G_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3938          newtonPoly= mapUp (newtonPoly, alpha, V_buf, prim_elem, im_prim_elem,
3939                              source, dest);
3940          ppA= mapUp (ppA, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3941          ppB= mapUp (ppB, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3942          gcdlcAlcB= mapUp (gcdlcAlcB, alpha, V_buf, prim_elem, im_prim_elem,
3943                            source, dest);
3944          fail= false;
3945          random_element= randomElement (m, V_buf, l, fail );
3946          DEBOUTLN (cerr, "fail= " << fail);
3947          CFList list;
3948          TIMING_START (gcd_recursion);
3949          if (LC (skeleton).inCoeffDomain())
3950            G_random_element=
3951            monicSparseInterpol (ppA (random_element, x), ppB (random_element, x),
3952                                skeleton, V_buf, sparseFail, coeffMonoms,
3953                                Monoms);
3954          else
3955            G_random_element=
3956            nonMonicSparseInterpol (ppA(random_element,x), ppB(random_element,x),
3957                                    skeleton, V_buf, sparseFail,
3958                                    coeffMonoms, Monoms);
3959          TIMING_END_AND_PRINT (gcd_recursion,
3960                                "time for recursive call: ");
3961          DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3962          alpha= V_buf;
3963        }
3964
3965        if (sparseFail)
3966          break;
3967
3968        if (!G_random_element.inCoeffDomain())
3969          d0= totaldegree (G_random_element, Variable(2),
3970                           Variable (G_random_element.level()));
3971        else
3972          d0= 0;
3973
3974        if (d0 == 0)
3975        {
3976          if (inextension)
3977            prune (cleanUp);
3978          return N(gcdcAcB);
3979        }
3980        if (d0 >  d)
3981        {
3982          if (!find (l, random_element))
3983            l.append (random_element);
3984          continue;
3985        }
3986
3987        G_random_element=
3988        (gcdlcAlcB(random_element, x)/uni_lcoeff (G_random_element))
3989        * G_random_element;
3990
3991        if (!G_random_element.inCoeffDomain())
3992          d0= totaldegree (G_random_element, Variable(2),
3993                           Variable (G_random_element.level()));
3994        else
3995          d0= 0;
3996
3997        if (d0 <  d)
3998        {
3999          m= gcdlcAlcB;
4000          newtonPoly= 1;
4001          G_m= 0;
4002          d= d0;
4003        }
4004
4005        TIMING_START (newton_interpolation);
4006        H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
4007        TIMING_END_AND_PRINT (newton_interpolation,
4008                              "time for newton interpolation: ");
4009
4010        //termination test
4011        if (uni_lcoeff (H) == gcdlcAlcB)
4012        {
4013          cH= uni_content (H);
4014          ppH= H/cH;
4015          ppH /= Lc (ppH);
4016          DEBOUTLN (cerr, "ppH= " << ppH);
4017          if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
4018          {
4019            if (inextension)
4020              prune (cleanUp);
4021            return N(gcdcAcB*ppH);
4022          }
4023        }
4024
4025        G_m= H;
4026        newtonPoly= newtonPoly*(x - random_element);
4027        m= m*(x - random_element);
4028        if (!find (l, random_element))
4029          l.append (random_element);
4030
4031      } while (1); //end of second do
4032    }
4033    else
4034    {
4035      if (inextension)
4036        prune (cleanUp);
4037      return N(gcdcAcB*modGCDFp (ppA, ppB));
4038    }
4039  } while (1); //end of first do
4040}
4041#endif
4042
4043TIMING_DEFINE_PRINT(modZ_termination)
4044TIMING_DEFINE_PRINT(modZ_recursion)
4045
4046/// modular gcd algorithm, see Keith, Czapor, Geddes "Algorithms for Computer
4047/// Algebra", Algorithm 7.1
4048CanonicalForm modGCDZ ( const CanonicalForm & FF, const CanonicalForm & GG )
4049{
4050  CanonicalForm f, g, cl, q(0), Dp, newD, D, newq, newqh;
4051  int p, i, dp_deg, d_deg=-1;
4052
4053  CanonicalForm cd ( bCommonDen( FF ));
4054  f=cd*FF;
4055  Variable x= Variable (1);
4056  CanonicalForm cf, cg;
4057  cf= icontent (f);
4058  f /= cf;
4059  //cd = bCommonDen( f ); f *=cd;
4060  //f /=vcontent(f,Variable(1));
4061
4062  cd = bCommonDen( GG );
4063  g=cd*GG;
4064  cg= icontent (g);
4065  g /= cg;
4066  //cd = bCommonDen( g ); g *=cd;
4067  //g /=vcontent(g,Variable(1));
4068
4069  CanonicalForm Dn, test= 0;
4070  CanonicalForm lcf, lcg;
4071  lcf= f.lc();
4072  lcg= g.lc();
4073  cl =  gcd (f.lc(),g.lc());
4074  CanonicalForm gcdcfcg= gcd (cf, cg);
4075  CanonicalForm fp, gp;
4076  CanonicalForm b= 1;
4077  int minCommonDeg= 0;
4078  for (i= tmax (f.level(), g.level()); i > 0; i--)
4079  {
4080    if (degree (f, i) <= 0 || degree (g, i) <= 0)
4081      continue;
4082    else
4083    {
4084      minCommonDeg= tmin (degree (g, i), degree (f, i));
4085      break;
4086    }
4087  }
4088  if (i == 0)
4089    return gcdcfcg;
4090  for (; i > 0; i--)
4091  {
4092    if (degree (f, i) <= 0 || degree (g, i) <= 0)
4093      continue;
4094    else
4095      minCommonDeg= tmin (minCommonDeg, tmin (degree (g, i), degree (f, i)));
4096  }
4097  b= 2*tmin (maxNorm (f), maxNorm (g))*abs (cl)*
4098     power (CanonicalForm (2), minCommonDeg);
4099  bool equal= false;
4100  i = cf_getNumBigPrimes() - 1;
4101
4102  CanonicalForm cof, cog, cofp, cogp, newCof, newCog, cofn, cogn, cDn;
4103  int maxNumVars= tmax (getNumVars (f), getNumVars (g));
4104  //Off (SW_RATIONAL);
4105  while ( true )
4106  {
4107    p = cf_getBigPrime( i );
4108    i--;
4109    while ( i >= 0 && mod( cl*(lc(f)/cl)*(lc(g)/cl), p ) == 0 )
4110    {
4111      p = cf_getBigPrime( i );
4112      i--;
4113    }
4114    //printf("try p=%d\n",p);
4115    setCharacteristic( p );
4116    fp= mapinto (f);
4117    gp= mapinto (g);
4118    TIMING_START (modZ_recursion)
4119#ifdef HAVE_NTL
4120    if (size (fp)/maxNumVars > 500 && size (gp)/maxNumVars > 500)
4121      Dp = modGCDFp (fp, gp, cofp, cogp);
4122    else
4123    {
4124      Dp= gcd_poly (fp, gp);
4125      cofp= fp/Dp;
4126      cogp= gp/Dp;
4127    }
4128#else
4129    Dp= gcd_poly (fp, gp);
4130    cofp= fp/Dp;
4131    cogp= gp/Dp;
4132#endif
4133    TIMING_END_AND_PRINT (modZ_recursion,
4134                          "time for gcd mod p in modular gcd: ");
4135    Dp /=Dp.lc();
4136    Dp *= mapinto (cl);
4137    cofp /= lc (cofp);
4138    cofp *= mapinto (lcf);
4139    cogp /= lc (cogp);
4140    cogp *= mapinto (lcg);
4141    setCharacteristic( 0 );
4142    dp_deg=totaldegree(Dp);
4143    if ( dp_deg == 0 )
4144    {
4145      //printf(" -> 1\n");
4146      return CanonicalForm(gcdcfcg);
4147    }
4148    if ( q.isZero() )
4149    {
4150      D = mapinto( Dp );
4151      cof= mapinto (cofp);
4152      cog= mapinto (cogp);
4153      d_deg=dp_deg;
4154      q = p;
4155      Dn= balance_p (D, p);
4156      cofn= balance_p (cof, p);
4157      cogn= balance_p (cog, p);
4158    }
4159    else
4160    {
4161      if ( dp_deg == d_deg )
4162      {
4163        chineseRemainder( D, q, mapinto( Dp ), p, newD, newq );
4164        chineseRemainder( cof, q, mapinto (cofp), p, newCof, newq);
4165        chineseRemainder( cog, q, mapinto (cogp), p, newCog, newq);
4166        cof= newCof;
4167        cog= newCog;
4168        newqh= newq/2;
4169        Dn= balance_p (newD, newq, newqh);
4170        cofn= balance_p (newCof, newq, newqh);
4171        cogn= balance_p (newCog, newq, newqh);
4172        if (test != Dn) //balance_p (newD, newq))
4173          test= balance_p (newD, newq);
4174        else
4175          equal= true;
4176        q = newq;
4177        D = newD;
4178      }
4179      else if ( dp_deg < d_deg )
4180      {
4181        //printf(" were all bad, try more\n");
4182        // all previous p's are bad primes
4183        q = p;
4184        D = mapinto( Dp );
4185        cof= mapinto (cof);
4186        cog= mapinto (cog);
4187        d_deg=dp_deg;
4188        test= 0;
4189        equal= false;
4190        Dn= balance_p (D, p);
4191        cofn= balance_p (cof, p);
4192        cogn= balance_p (cog, p);
4193      }
4194      else
4195      {
4196        //printf(" was bad, try more\n");
4197      }
4198      //else dp_deg > d_deg: bad prime
4199    }
4200    if ( i >= 0 )
4201    {
4202      cDn= icontent (Dn);
4203      Dn /= cDn;
4204      cofn /= cl/cDn;
4205      //cofn /= icontent (cofn);
4206      cogn /= cl/cDn;
4207      //cogn /= icontent (cogn);
4208      TIMING_START (modZ_termination);
4209      if ((terminationTest (f,g, cofn, cogn, Dn)) ||
4210          ((equal || q > b) && fdivides (Dn, f) && fdivides (Dn, g)))
4211      {
4212        TIMING_END_AND_PRINT (modZ_termination,
4213                            "time for successful termination in modular gcd: ");
4214        //printf(" -> success\n");
4215        return Dn*gcdcfcg;
4216      }
4217      TIMING_END_AND_PRINT (modZ_termination,
4218                          "time for unsuccessful termination in modular gcd: ");
4219      equal= false;
4220      //else: try more primes
4221    }
4222    else
4223    { // try other method
4224      //printf("try other gcd\n");
4225      Off(SW_USE_CHINREM_GCD);
4226      D=gcd_poly( f, g );
4227      On(SW_USE_CHINREM_GCD);
4228      return D*gcdcfcg;
4229    }
4230  }
4231}
4232
4233
4234#endif
Note: See TracBrowser for help on using the repository browser.