source: git/factory/cfModGcd.cc @ e64ec88

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