source: git/factory/cfModGcd.cc @ 8d1432e

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