source: git/factory/cfModGcd.cc @ ea1479

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