source: git/factory/cf_gcd_smallp.cc @ 102a88e

spielwiese
Last change on this file since 102a88e was 52cccfd, checked in by Martin Lee <martinlee84@…>, 11 years ago
fix: int overflow
  • Property mode set to 100644
File size: 133.2 KB
Line 
1// -*- c++ -*-
2//*****************************************************************************
3/** @file cf_gcd_smallp.cc
4 *
5 * @author Martin Lee
6 * @date 22.10.2009
7 *
8 * This file implements the GCD of two polynomials over \f$ F_{p} \f$ ,
9 * \f$ F_{p}(\alpha ) \f$ or GF based on Alg. 7.2. as described in "Algorithms
10 * for Computer Algebra" by Geddes, Czapor, Labahnn
11 *
12 * @par Copyright:
13 *   (c) by The SINGULAR Team, see LICENSE file
14 *
15**/
16//*****************************************************************************
17
18#include "config.h"
19
20#include "cf_assert.h"
21#include "debug.h"
22#include "timing.h"
23
24#include "canonicalform.h"
25#include "algext.h"
26#include "cf_map.h"
27#include "cf_util.h"
28#include "templates/ftmpl_functions.h"
29#include "cf_random.h"
30#include "cf_reval.h"
31#include "facHensel.h"
32#include "cf_iter.h"
33#include "cfNewtonPolygon.h"
34#include "cf_algorithm.h"
35#include "cf_primes.h"
36
37// iinline helper functions:
38#include "cf_map_ext.h"
39
40#ifdef HAVE_NTL
41#include <NTLconvert.h>
42#endif
43
44#ifdef HAVE_FLINT
45#include "FLINTconvert.h"
46#endif
47
48#include "cf_gcd_smallp.h"
49
50TIMING_DEFINE_PRINT(gcd_recursion)
51TIMING_DEFINE_PRINT(newton_interpolation)
52TIMING_DEFINE_PRINT(termination_test)
53TIMING_DEFINE_PRINT(ez_p_compress)
54TIMING_DEFINE_PRINT(ez_p_hensel_lift)
55TIMING_DEFINE_PRINT(ez_p_eval)
56TIMING_DEFINE_PRINT(ez_p_content)
57
58bool
59terminationTest (const CanonicalForm& F, const CanonicalForm& G,
60                 const CanonicalForm& coF, const CanonicalForm& coG,
61                 const CanonicalForm& cand)
62{
63  CanonicalForm LCCand= abs (LC (cand));
64  if (LCCand*abs (LC (coF)) == abs (LC (F)))
65  {
66    if (LCCand*abs (LC (coG)) == abs (LC (G)))
67    {
68      if (abs (cand)*abs (coF) == abs (F))
69      {
70        if (abs (cand)*abs (coG) == abs (G))
71          return true;
72      }
73      return false;
74    }
75    return false;
76  }
77  return false;
78}
79
80#ifdef HAVE_NTL
81
82static const double log2exp= 1.442695041;
83
84/// compressing two polynomials F and G, M is used for compressing,
85/// N to reverse the compression
86int myCompress (const CanonicalForm& F, const CanonicalForm& G, CFMap & M,
87                CFMap & N, bool topLevel)
88{
89  int n= tmax (F.level(), G.level());
90  int * degsf= new int [n + 1];
91  int * degsg= new int [n + 1];
92
93  for (int i = 0; i <= n; i++)
94    degsf[i]= degsg[i]= 0;
95
96  degsf= degrees (F, degsf);
97  degsg= degrees (G, degsg);
98
99  int both_non_zero= 0;
100  int f_zero= 0;
101  int g_zero= 0;
102  int both_zero= 0;
103
104  if (topLevel)
105  {
106    for (int i= 1; i <= n; i++)
107    {
108      if (degsf[i] != 0 && degsg[i] != 0)
109      {
110        both_non_zero++;
111        continue;
112      }
113      if (degsf[i] == 0 && degsg[i] != 0 && i <= G.level())
114      {
115        f_zero++;
116        continue;
117      }
118      if (degsg[i] == 0 && degsf[i] && i <= F.level())
119      {
120        g_zero++;
121        continue;
122      }
123    }
124
125    if (both_non_zero == 0)
126    {
127      delete [] degsf;
128      delete [] degsg;
129      return 0;
130    }
131
132    // map Variables which do not occur in both polynomials to higher levels
133    int k= 1;
134    int l= 1;
135    for (int i= 1; i <= n; i++)
136    {
137      if (degsf[i] != 0 && degsg[i] == 0 && i <= F.level())
138      {
139        if (k + both_non_zero != i)
140        {
141          M.newpair (Variable (i), Variable (k + both_non_zero));
142          N.newpair (Variable (k + both_non_zero), Variable (i));
143        }
144        k++;
145      }
146      if (degsf[i] == 0 && degsg[i] != 0 && i <= G.level())
147      {
148        if (l + g_zero + both_non_zero != i)
149        {
150          M.newpair (Variable (i), Variable (l + g_zero + both_non_zero));
151          N.newpair (Variable (l + g_zero + both_non_zero), Variable (i));
152        }
153        l++;
154      }
155    }
156
157    // sort Variables x_{i} in increasing order of
158    // min(deg_{x_{i}}(f),deg_{x_{i}}(g))
159    int m= tmax (F.level(), G.level());
160    int min_max_deg;
161    k= both_non_zero;
162    l= 0;
163    int i= 1;
164    while (k > 0)
165    {
166      if (degsf [i] != 0 && degsg [i] != 0)
167        min_max_deg= tmax (degsf[i], degsg[i]);
168      else
169        min_max_deg= 0;
170      while (min_max_deg == 0)
171      {
172        i++;
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      }
178      for (int j= i + 1; j <=  m; j++)
179      {
180        if (degsf[j] != 0 && degsg [j] != 0 &&
181            tmax (degsf[j],degsg[j]) <= min_max_deg)
182        {
183          min_max_deg= tmax (degsf[j], degsg[j]);
184          l= j;
185        }
186      }
187      if (l != 0)
188      {
189        if (l != k)
190        {
191          M.newpair (Variable (l), Variable(k));
192          N.newpair (Variable (k), Variable(l));
193          degsf[l]= 0;
194          degsg[l]= 0;
195          l= 0;
196        }
197        else
198        {
199          degsf[l]= 0;
200          degsg[l]= 0;
201          l= 0;
202        }
203      }
204      else if (l == 0)
205      {
206        if (i != k)
207        {
208          M.newpair (Variable (i), Variable (k));
209          N.newpair (Variable (k), Variable (i));
210          degsf[i]= 0;
211          degsg[i]= 0;
212        }
213        else
214        {
215          degsf[i]= 0;
216          degsg[i]= 0;
217        }
218        i++;
219      }
220      k--;
221    }
222  }
223  else
224  {
225    //arrange Variables such that no gaps occur
226    for (int i= 1; i <= n; i++)
227    {
228      if (degsf[i] == 0 && degsg[i] == 0)
229      {
230        both_zero++;
231        continue;
232      }
233      else
234      {
235        if (both_zero != 0)
236        {
237          M.newpair (Variable (i), Variable (i - both_zero));
238          N.newpair (Variable (i - both_zero), Variable (i));
239        }
240      }
241    }
242  }
243
244  delete [] degsf;
245  delete [] degsg;
246
247  return 1;
248}
249
250static inline CanonicalForm
251uni_content (const CanonicalForm & F);
252
253CanonicalForm
254uni_content (const CanonicalForm& F, const Variable& x)
255{
256  if (F.inCoeffDomain())
257    return F.genOne();
258  if (F.level() == x.level() && F.isUnivariate())
259    return F;
260  if (F.level() != x.level() && F.isUnivariate())
261    return F.genOne();
262
263  if (x.level() != 1)
264  {
265    CanonicalForm f= swapvar (F, x, Variable (1));
266    CanonicalForm result= uni_content (f);
267    return swapvar (result, x, Variable (1));
268  }
269  else
270    return uni_content (F);
271}
272
273/// compute the content of F, where F is considered as an element of
274/// \f$ R[x_{1}][x_{2},\ldots ,x_{n}] \f$
275static inline CanonicalForm
276uni_content (const CanonicalForm & F)
277{
278  if (F.inBaseDomain())
279    return F.genOne();
280  if (F.level() == 1 && F.isUnivariate())
281    return F;
282  if (F.level() != 1 && F.isUnivariate())
283    return F.genOne();
284  if (degree (F,1) == 0) return F.genOne();
285
286  int l= F.level();
287  if (l == 2)
288    return content(F);
289  else
290  {
291    CanonicalForm pol, c = 0;
292    CFIterator i = F;
293    for (; i.hasTerms(); i++)
294    {
295      pol= i.coeff();
296      pol= uni_content (pol);
297      c= gcd (c, pol);
298      if (c.isOne())
299        return c;
300    }
301    return c;
302  }
303}
304
305CanonicalForm
306extractContents (const CanonicalForm& F, const CanonicalForm& G,
307                 CanonicalForm& contentF, CanonicalForm& contentG,
308                 CanonicalForm& ppF, CanonicalForm& ppG, const int d)
309{
310  CanonicalForm uniContentF, uniContentG, gcdcFcG;
311  contentF= 1;
312  contentG= 1;
313  ppF= F;
314  ppG= G;
315  CanonicalForm result= 1;
316  for (int i= 1; i <= d; i++)
317  {
318    uniContentF= uni_content (F, Variable (i));
319    uniContentG= uni_content (G, Variable (i));
320    gcdcFcG= gcd (uniContentF, uniContentG);
321    contentF *= uniContentF;
322    contentG *= uniContentG;
323    ppF /= uniContentF;
324    ppG /= uniContentG;
325    result *= gcdcFcG;
326  }
327  return result;
328}
329
330/// compute the leading coefficient of F, where F is considered as an element
331/// of \f$ R[x_{1}][x_{2},\ldots ,x_{n}] \f$, order on Mon (x_{2},\ldots ,x_{n})
332/// is dp.
333static inline
334CanonicalForm uni_lcoeff (const CanonicalForm& F)
335{
336  if (F.level() > 1)
337  {
338    Variable x= Variable (2);
339    int deg= totaldegree (F, x, F.mvar());
340    for (CFIterator i= F; i.hasTerms(); i++)
341    {
342      if (i.exp() + totaldegree (i.coeff(), x, i.coeff().mvar()) == deg)
343        return uni_lcoeff (i.coeff());
344    }
345  }
346  return F;
347}
348
349/// Newton interpolation - Incremental algorithm.
350/// Given a list of values alpha_i and a list of polynomials u_i, 1 <= i <= n,
351/// computes the interpolation polynomial assuming that
352/// the polynomials in u are the results of evaluating the variabe x
353/// of the unknown polynomial at the alpha values.
354/// This incremental version receives only the values of alpha_n and u_n and
355/// the previous interpolation polynomial for points 1 <= i <= n-1, and computes
356/// the polynomial interpolating in all the points.
357/// newtonPoly must be equal to (x - alpha_1) * ... * (x - alpha_{n-1})
358static inline CanonicalForm
359newtonInterp(const CanonicalForm alpha, const CanonicalForm u,
360             const CanonicalForm newtonPoly, const CanonicalForm oldInterPoly,
361             const Variable & x)
362{
363  CanonicalForm interPoly;
364
365  interPoly= oldInterPoly + ((u - oldInterPoly(alpha, x))/newtonPoly(alpha, x))
366             *newtonPoly;
367  return interPoly;
368}
369
370/// compute a random element a of \f$ F_{p}(\alpha ) \f$ ,
371/// s.t. F(a) \f$ \neq 0 \f$ , F is a univariate polynomial, returns
372/// fail if there are no field elements left which have not been used before
373static inline CanonicalForm
374randomElement (const CanonicalForm & F, const Variable & alpha, CFList & list,
375               bool & fail)
376{
377  fail= false;
378  Variable x= F.mvar();
379  AlgExtRandomF genAlgExt (alpha);
380  FFRandom genFF;
381  CanonicalForm random, mipo;
382  mipo= getMipo (alpha);
383  int p= getCharacteristic ();
384  int d= degree (mipo);
385  double bound= pow ((double) p, (double) d);
386  do
387  {
388    if (list.length() == bound)
389    {
390      fail= true;
391      break;
392    }
393    if (list.length() < p)
394    {
395      random= genFF.generate();
396      while (find (list, random))
397        random= genFF.generate();
398    }
399    else
400    {
401      random= genAlgExt.generate();
402      while (find (list, random))
403        random= genAlgExt.generate();
404    }
405    if (F (random, x) == 0)
406    {
407      list.append (random);
408      continue;
409    }
410  } while (find (list, random));
411  return random;
412}
413
414static inline
415Variable chooseExtension (const Variable & alpha)
416{
417  if (fac_NTL_char != getCharacteristic())
418  {
419    fac_NTL_char= getCharacteristic();
420    zz_p::init (getCharacteristic());
421  }
422  zz_pX NTLIrredpoly;
423  int i, m;
424  // extension of F_p needed
425  if (alpha.level() == 1)
426  {
427    i= 1;
428    m= 2;
429  } //extension of F_p(alpha)
430  if (alpha.level() != 1)
431  {
432    i= 4;
433    m= degree (getMipo (alpha));
434  }
435  BuildIrred (NTLIrredpoly, i*m);
436  CanonicalForm newMipo= convertNTLzzpX2CF (NTLIrredpoly, Variable (1));
437  return rootOf (newMipo);
438}
439
440CanonicalForm
441GCD_Fp_extension (const CanonicalForm& F, const CanonicalForm& G,
442                  CanonicalForm& coF, CanonicalForm& coG,
443                  Variable & alpha, CFList& l, bool& topLevel);
444
445CanonicalForm
446GCD_Fp_extension (const CanonicalForm& F, const CanonicalForm& G,
447                  Variable & alpha, CFList& l, bool& topLevel)
448{
449  CanonicalForm dummy1, dummy2;
450  CanonicalForm result= GCD_Fp_extension (F, G, dummy1, dummy2, alpha, l,
451                                          topLevel);
452  return result;
453}
454
455/// GCD of F and G over \f$ F_{p}(\alpha ) \f$ ,
456/// l and topLevel are only used internally, output is monic
457/// based on Alg. 7.2. as described in "Algorithms for
458/// Computer Algebra" by Geddes, Czapor, Labahn
459CanonicalForm
460GCD_Fp_extension (const CanonicalForm& F, const CanonicalForm& G,
461                  CanonicalForm& coF, CanonicalForm& coG,
462                  Variable & alpha, CFList& l, bool& topLevel)
463{
464  CanonicalForm A= F;
465  CanonicalForm B= G;
466  if (F.isZero() && degree(G) > 0)
467  {
468    coF= 0;
469    coG= Lc (G);
470    return G/Lc(G);
471  }
472  else if (G.isZero() && degree (F) > 0)
473  {
474    coF= Lc (F);
475    coG= 0;
476    return F/Lc(F);
477  }
478  else if (F.isZero() && G.isZero())
479  {
480    coF= coG= 0;
481    return F.genOne();
482  }
483  if (F.inBaseDomain() || G.inBaseDomain())
484  {
485    coF= F;
486    coG= G;
487    return F.genOne();
488  }
489  if (F.isUnivariate() && fdivides(F, G, coG))
490  {
491    coF= Lc (F);
492    return F/Lc(F);
493  }
494  if (G.isUnivariate() && fdivides(G, F, coF))
495  {
496    coG= Lc (G);
497    return G/Lc(G);
498  }
499  if (F == G)
500  {
501    coF= coG= Lc (F);
502    return F/Lc(F);
503  }
504
505  CFMap M,N;
506  int best_level= myCompress (A, B, M, N, topLevel);
507
508  if (best_level == 0)
509  {
510    coF= F;
511    coG= G;
512    return B.genOne();
513  }
514
515  A= M(A);
516  B= M(B);
517
518  Variable x= Variable(1);
519
520  //univariate case
521  if (A.isUnivariate() && B.isUnivariate())
522  {
523    CanonicalForm result= gcd (A, B);
524    coF= N (A/result);
525    coG= N (B/result);
526    return N (result);
527  }
528
529  CanonicalForm cA, cB;    // content of A and B
530  CanonicalForm ppA, ppB;    // primitive part of A and B
531  CanonicalForm gcdcAcB;
532
533  cA = uni_content (A);
534  cB = uni_content (B);
535  gcdcAcB= gcd (cA, cB);
536  ppA= A/cA;
537  ppB= B/cB;
538
539  int sizeNewtonPolyg;
540  int ** newtonPolyg= NULL;
541  mat_ZZ MM;
542  vec_ZZ V;
543  bool compressConvexDense= false; //(ppA.level() == 2 && ppB.level() == 2);
544  if (compressConvexDense)
545  {
546    CanonicalForm bufcA= cA;
547    CanonicalForm bufcB= cB;
548    cA= content (ppA, 1);
549    cB= content (ppB, 1);
550    ppA /= cA;
551    ppB /= cB;
552    gcdcAcB *= gcd (cA, cB);
553    cA *= bufcA;
554    cB *= bufcB;
555    if (ppA.isUnivariate() || ppB.isUnivariate())
556    {
557      if (ppA.level() == ppB.level())
558      {
559        CanonicalForm result= gcd (ppA, ppB);
560        coF= N ((ppA/result)*(cA/gcdcAcB));
561        coG= N ((ppB/result)*(cB/gcdcAcB));
562        return N (result*gcdcAcB);
563      }
564      else
565      {
566        coF= N (ppA*(cA/gcdcAcB));
567        coG= N (ppB*(cB/gcdcAcB));
568        return N (gcdcAcB);
569      }
570    }
571
572    newtonPolyg= newtonPolygon (ppA,ppB, sizeNewtonPolyg);
573    convexDense (newtonPolyg, sizeNewtonPolyg, MM, V);
574
575    for (int i= 0; i < sizeNewtonPolyg; i++)
576      delete [] newtonPolyg[i];
577    delete [] newtonPolyg;
578
579    ppA= compress (ppA, MM, V, false);
580    ppB= compress (ppB, MM, V, false);
581    MM= inv (MM);
582
583    if (ppA.isUnivariate() && ppB.isUnivariate())
584    {
585      if (ppA.level() == ppB.level())
586      {
587        CanonicalForm result= gcd (ppA, ppB);
588        coF= N (decompress ((ppA/result), MM, V)*(cA/gcdcAcB));
589        coG= N (decompress ((ppB/result), MM, V)*(cB/gcdcAcB));
590        return N (decompress (result, MM, V)*gcdcAcB);
591      }
592      else
593      {
594        coF= N (decompress (ppA, MM, V));
595        coG= N (decompress (ppB, MM, V));
596        return N (gcdcAcB);
597      }
598    }
599  }
600
601  CanonicalForm lcA, lcB;  // leading coefficients of A and B
602  CanonicalForm gcdlcAlcB;
603
604  lcA= uni_lcoeff (ppA);
605  lcB= uni_lcoeff (ppB);
606
607  /*if (fdivides (lcA, lcB))
608  {
609    if (fdivides (A, B))
610      return F/Lc(F);
611  }
612  if (fdivides (lcB, lcA))
613  {
614    if (fdivides (B, A))
615      return G/Lc(G);
616  }*/
617
618  gcdlcAlcB= gcd (lcA, lcB);
619
620  int d= totaldegree (ppA, Variable(2), Variable (ppA.level()));
621
622  if (d == 0)
623  {
624    coF= N (ppA*(cA/gcdcAcB));
625    coG= N (ppB*(cB/gcdcAcB));
626    return N(gcdcAcB);
627  }
628
629  int d0= totaldegree (ppB, Variable(2), Variable (ppB.level()));
630  if (d0 < d)
631    d= d0;
632  if (d == 0)
633  {
634    coF= N (ppA*(cA/gcdcAcB));
635    coG= N (ppB*(cB/gcdcAcB));
636    return N(gcdcAcB);
637  }
638
639  CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH;
640  CanonicalForm newtonPoly, coF_random_element, coG_random_element, coF_m,
641                coG_m, ppCoF, ppCoG;
642
643  newtonPoly= 1;
644  m= gcdlcAlcB;
645  G_m= 0;
646  coF= 0;
647  coG= 0;
648  H= 0;
649  bool fail= false;
650  topLevel= false;
651  bool inextension= false;
652  Variable V_buf= alpha;
653  CanonicalForm prim_elem, im_prim_elem;
654  CFList source, dest;
655  int bound1= degree (ppA, 1);
656  int bound2= degree (ppB, 1);
657  do
658  {
659    random_element= randomElement (m*lcA*lcB, V_buf, l, fail);
660    if (fail)
661    {
662      source= CFList();
663      dest= CFList();
664
665      Variable V_buf3= V_buf;
666      V_buf= chooseExtension (V_buf);
667      bool prim_fail= false;
668      Variable V_buf2;
669      prim_elem= primitiveElement (alpha, V_buf2, prim_fail);
670
671      if (V_buf3 != alpha)
672      {
673        m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
674        G_m= mapDown (G_m, prim_elem, im_prim_elem, alpha, source, dest);
675        coF_m= mapDown (coF_m, prim_elem, im_prim_elem, alpha, source, dest);
676        coG_m= mapDown (coG_m, prim_elem, im_prim_elem, alpha, source, dest);
677        newtonPoly= mapDown (newtonPoly, prim_elem, im_prim_elem, alpha,
678                             source, dest);
679        ppA= mapDown (ppA, prim_elem, im_prim_elem, alpha, source, dest);
680        ppB= mapDown (ppB, prim_elem, im_prim_elem, alpha, source, dest);
681        gcdlcAlcB= mapDown (gcdlcAlcB, prim_elem, im_prim_elem, alpha,
682                            source, dest);
683        lcA= mapDown (lcA, prim_elem, im_prim_elem, alpha, source, dest);
684        lcB= mapDown (lcB, prim_elem, im_prim_elem, alpha, source, dest);
685        for (CFListIterator i= l; i.hasItem(); i++)
686          i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, alpha,
687                                source, dest);
688      }
689
690      ASSERT (!prim_fail, "failure in integer factorizer");
691      if (prim_fail)
692        ; //ERROR
693      else
694        im_prim_elem= mapPrimElem (prim_elem, alpha, V_buf);
695
696      DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (alpha));
697      DEBOUTLN (cerr, "getMipo (V_buf2)= " << getMipo (V_buf2));
698      inextension= true;
699      for (CFListIterator i= l; i.hasItem(); i++)
700        i.getItem()= mapUp (i.getItem(), alpha, V_buf, prim_elem,
701                             im_prim_elem, source, dest);
702      m= mapUp (m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
703      G_m= mapUp (G_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
704      coF_m= mapUp (coF_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
705      coG_m= mapUp (coG_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
706      newtonPoly= mapUp (newtonPoly, alpha, V_buf, prim_elem, im_prim_elem,
707                          source, dest);
708      ppA= mapUp (ppA, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
709      ppB= mapUp (ppB, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
710      gcdlcAlcB= mapUp (gcdlcAlcB, alpha, V_buf, prim_elem, im_prim_elem,
711                         source, dest);
712      lcA= mapUp (lcA, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
713      lcB= mapUp (lcB, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
714
715      fail= false;
716      random_element= randomElement (m*lcA*lcB, V_buf, l, fail );
717      DEBOUTLN (cerr, "fail= " << fail);
718      CFList list;
719      TIMING_START (gcd_recursion);
720      G_random_element=
721      GCD_Fp_extension (ppA (random_element, x), ppB (random_element, x),
722                        coF_random_element, coG_random_element, V_buf,
723                        list, topLevel);
724      TIMING_END_AND_PRINT (gcd_recursion,
725                            "time for recursive call: ");
726      DEBOUTLN (cerr, "G_random_element= " << G_random_element);
727    }
728    else
729    {
730      CFList list;
731      TIMING_START (gcd_recursion);
732      G_random_element=
733      GCD_Fp_extension (ppA(random_element, x), ppB(random_element, x),
734                        coF_random_element, coG_random_element, V_buf,
735                        list, topLevel);
736      TIMING_END_AND_PRINT (gcd_recursion,
737                            "time for recursive call: ");
738      DEBOUTLN (cerr, "G_random_element= " << G_random_element);
739    }
740
741    if (!G_random_element.inCoeffDomain())
742      d0= totaldegree (G_random_element, Variable(2),
743                       Variable (G_random_element.level()));
744    else
745      d0= 0;
746
747    if (d0 == 0)
748    {
749      coF= N (ppA*(cA/gcdcAcB));
750      coG= N (ppB*(cB/gcdcAcB));
751      return N(gcdcAcB);
752    }
753    if (d0 >  d)
754    {
755      if (!find (l, random_element))
756        l.append (random_element);
757      continue;
758    }
759
760    G_random_element=
761    (gcdlcAlcB(random_element, x)/uni_lcoeff (G_random_element))
762    * G_random_element;
763    coF_random_element= (lcA(random_element,x)/uni_lcoeff(coF_random_element))
764                        *coF_random_element;
765    coG_random_element= (lcB(random_element,x)/uni_lcoeff(coG_random_element))
766                        *coG_random_element;
767
768    if (!G_random_element.inCoeffDomain())
769      d0= totaldegree (G_random_element, Variable(2),
770                       Variable (G_random_element.level()));
771    else
772      d0= 0;
773
774    if (d0 <  d)
775    {
776      m= gcdlcAlcB;
777      newtonPoly= 1;
778      G_m= 0;
779      d= d0;
780      coF_m= 0;
781      coG_m= 0;
782    }
783
784    TIMING_START (newton_interpolation);
785    H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
786    coF= newtonInterp (random_element, coF_random_element, newtonPoly, coF_m,x);
787    coG= newtonInterp (random_element, coG_random_element, newtonPoly, coG_m,x);
788    TIMING_END_AND_PRINT (newton_interpolation,
789                          "time for newton interpolation: ");
790
791    //termination test
792    if ((uni_lcoeff (H) == gcdlcAlcB) || (G_m == H))
793    {
794      TIMING_START (termination_test);
795      if (gcdlcAlcB.isOne())
796        cH= 1;
797      else
798        cH= uni_content (H);
799      ppH= H/cH;
800      ppH /= Lc (ppH);
801      CanonicalForm lcppH= gcdlcAlcB/cH;
802      CanonicalForm ccoF= lcppH/Lc (lcppH);
803      CanonicalForm ccoG= lcppH/Lc (lcppH);
804      ppCoF= coF/ccoF;
805      ppCoG= coG/ccoG;
806      if (inextension)
807      {
808        if (((degree (ppCoF,1)+degree (ppH,1) == bound1) &&
809             (degree (ppCoG,1)+degree (ppH,1) == bound2) &&
810             terminationTest (ppA, ppB, ppCoF, ppCoG, ppH)) || 
811             (fdivides (ppH, ppA, ppCoF) && fdivides (ppH, ppB, ppCoG)))
812        {
813          CFList u, v;
814          DEBOUTLN (cerr, "ppH before mapDown= " << ppH);
815          ppH= mapDown (ppH, prim_elem, im_prim_elem, alpha, u, v);
816          ppCoF= mapDown (ppCoF, prim_elem, im_prim_elem, alpha, u, v);
817          ppCoG= mapDown (ppCoG, prim_elem, im_prim_elem, alpha, u, v);
818          DEBOUTLN (cerr, "ppH after mapDown= " << ppH);
819          if (compressConvexDense)
820          {
821            ppH= decompress (ppH, MM, V);
822            ppCoF= decompress (ppCoF, MM, V);
823            ppCoG= decompress (ppCoG, MM, V);
824          }
825          coF= N ((cA/gcdcAcB)*ppCoF);
826          coG= N ((cB/gcdcAcB)*ppCoG);
827          TIMING_END_AND_PRINT (termination_test,
828                                "time for successful termination test Fq: ");
829          return N(gcdcAcB*ppH);
830        }
831      }
832      else if (((degree (ppCoF,1)+degree (ppH,1) == bound1) &&
833                (degree (ppCoG,1)+degree (ppH,1) == bound2) &&
834                terminationTest (ppA, ppB, ppCoF, ppCoG, ppH)) ||
835                (fdivides (ppH, ppA, ppCoF) && fdivides (ppH, ppB, ppCoG)))
836      {
837        if (compressConvexDense)
838        {
839          ppH= decompress (ppH, MM, V);
840          ppCoF= decompress (ppCoF, MM, V);
841          ppCoG= decompress (ppCoG, MM, V);
842        }
843        coF= N ((cA/gcdcAcB)*ppCoF);
844        coG= N ((cB/gcdcAcB)*ppCoG);
845        TIMING_END_AND_PRINT (termination_test,
846                              "time for successful termination test Fq: ");
847        return N(gcdcAcB*ppH);
848      }
849      TIMING_END_AND_PRINT (termination_test,
850                            "time for unsuccessful termination test Fq: ");
851    }
852
853    G_m= H;
854    coF_m= coF;
855    coG_m= coG;
856    newtonPoly= newtonPoly*(x - random_element);
857    m= m*(x - random_element);
858    if (!find (l, random_element))
859      l.append (random_element);
860  } while (1);
861}
862
863/// compute a random element a of GF, s.t. F(a) \f$ \neq 0 \f$ , F is a
864/// univariate polynomial, returns fail if there are no field elements left
865/// which have not been used before
866static inline
867CanonicalForm
868GFRandomElement (const CanonicalForm& F, CFList& list, bool& fail)
869{
870  fail= false;
871  Variable x= F.mvar();
872  GFRandom genGF;
873  CanonicalForm random;
874  int p= getCharacteristic();
875  int d= getGFDegree();
876  int bound= ipower (p, d);
877  do
878  {
879    if (list.length() == bound)
880    {
881      fail= true;
882      break;
883    }
884    if (list.length() < 1)
885      random= 0;
886    else
887    {
888      random= genGF.generate();
889      while (find (list, random))
890        random= genGF.generate();
891    }
892    if (F (random, x) == 0)
893    {
894      list.append (random);
895      continue;
896    }
897  } while (find (list, random));
898  return random;
899}
900
901CanonicalForm
902GCD_GF (const CanonicalForm& F, const CanonicalForm& G,
903        CanonicalForm& coF, CanonicalForm& coG,
904        CFList& l, bool& topLevel);
905
906CanonicalForm
907GCD_GF (const CanonicalForm& F, const CanonicalForm& G, CFList& l,
908        bool& topLevel)
909{
910  CanonicalForm dummy1, dummy2;
911  CanonicalForm result= GCD_GF (F, G, dummy1, dummy2, l, topLevel);
912  return result;
913}
914
915/// GCD of F and G over GF, based on Alg. 7.2. as described in "Algorithms for
916/// Computer Algebra" by Geddes, Czapor, Labahn
917/// Usually this algorithm will be faster than GCD_Fp_extension since GF has
918/// faster field arithmetics, however it might fail if the input is large since
919/// the size of the base field is bounded by 2^16, output is monic
920CanonicalForm
921GCD_GF (const CanonicalForm& F, const CanonicalForm& G,
922        CanonicalForm& coF, CanonicalForm& coG,
923        CFList& l, bool& topLevel)
924{
925  CanonicalForm A= F;
926  CanonicalForm B= G;
927  if (F.isZero() && degree(G) > 0)
928  {
929    coF= 0;
930    coG= Lc (G);
931    return G/Lc(G);
932  }
933  else if (G.isZero() && degree (F) > 0)
934  {
935    coF= Lc (F);
936    coG= 0;
937    return F/Lc(F);
938  }
939  else if (F.isZero() && G.isZero())
940  {
941    coF= coG= 0;
942    return F.genOne();
943  }
944  if (F.inBaseDomain() || G.inBaseDomain())
945  {
946    coF= F;
947    coG= G;
948    return F.genOne();
949  }
950  if (F.isUnivariate() && fdivides(F, G, coG))
951  {
952    coF= Lc (F);
953    return F/Lc(F);
954  }
955  if (G.isUnivariate() && fdivides(G, F, coF))
956  {
957    coG= Lc (G);
958    return G/Lc(G);
959  }
960  if (F == G)
961  {
962    coF= coG= Lc (F);
963    return F/Lc(F);
964  }
965
966  CFMap M,N;
967  int best_level= myCompress (A, B, M, N, topLevel);
968
969  if (best_level == 0)
970  {
971    coF= F;
972    coG= G;
973    return B.genOne();
974  }
975
976  A= M(A);
977  B= M(B);
978
979  Variable x= Variable(1);
980
981  //univariate case
982  if (A.isUnivariate() && B.isUnivariate())
983  {
984    CanonicalForm result= gcd (A, B);
985    coF= N (A/result);
986    coG= N (B/result);
987    return N (result);
988  }
989
990  CanonicalForm cA, cB;    // content of A and B
991  CanonicalForm ppA, ppB;    // primitive part of A and B
992  CanonicalForm gcdcAcB;
993
994  cA = uni_content (A);
995  cB = uni_content (B);
996  gcdcAcB= gcd (cA, cB);
997  ppA= A/cA;
998  ppB= B/cB;
999
1000  int sizeNewtonPolyg;
1001  int ** newtonPolyg= NULL;
1002  mat_ZZ MM;
1003  vec_ZZ V;
1004  bool compressConvexDense= false; //(ppA.level() == 2 && ppB.level() == 2);
1005  if (compressConvexDense)
1006  {
1007    CanonicalForm bufcA= cA;
1008    CanonicalForm bufcB= cB;
1009    cA= content (ppA, 1);
1010    cB= content (ppB, 1);
1011    ppA /= cA;
1012    ppB /= cB;
1013    gcdcAcB *= gcd (cA, cB);
1014    cA *= bufcA;
1015    cB *= bufcB;
1016    if (ppA.isUnivariate() || ppB.isUnivariate())
1017    {
1018      if (ppA.level() == ppB.level())
1019      {
1020        CanonicalForm result= gcd (ppA, ppB);
1021        coF= N ((ppA/result)*(cA/gcdcAcB));
1022        coG= N ((ppB/result)*(cB/gcdcAcB));
1023        return N (result*gcdcAcB);
1024      }
1025      else
1026      {
1027        coF= N (ppA*(cA/gcdcAcB));
1028        coG= N (ppB*(cB/gcdcAcB));
1029        return N (gcdcAcB);
1030      }
1031    }
1032
1033    newtonPolyg= newtonPolygon (ppA,ppB, sizeNewtonPolyg);
1034    convexDense (newtonPolyg, sizeNewtonPolyg, MM, V);
1035
1036    for (int i= 0; i < sizeNewtonPolyg; i++)
1037      delete [] newtonPolyg[i];
1038    delete [] newtonPolyg;
1039
1040    ppA= compress (ppA, MM, V, false);
1041    ppB= compress (ppB, MM, V, false);
1042    MM= inv (MM);
1043
1044    if (ppA.isUnivariate() && ppB.isUnivariate())
1045    {
1046      if (ppA.level() == ppB.level())
1047      {
1048        CanonicalForm result= gcd (ppA, ppB);
1049        coF= N (decompress ((ppA/result), MM, V)*(cA/gcdcAcB));
1050        coG= N (decompress ((ppB/result), MM, V)*(cB/gcdcAcB));
1051        return N (decompress (result, MM, V)*gcdcAcB);
1052      }
1053      else
1054      {
1055        coF= N (decompress (ppA, MM, V));
1056        coG= N (decompress (ppB, MM, V));
1057        return N (gcdcAcB);
1058      }
1059    }
1060  }
1061
1062  CanonicalForm lcA, lcB;  // leading coefficients of A and B
1063  CanonicalForm gcdlcAlcB;
1064
1065  lcA= uni_lcoeff (ppA);
1066  lcB= uni_lcoeff (ppB);
1067
1068  /*if (fdivides (lcA, lcB))
1069  {
1070    if (fdivides (ppA, ppB, coG))
1071    {
1072      coF= 1;
1073      if (compressConvexDense)
1074        coG= decompress (coG, MM, V);
1075      coG= N (coG*(cB/gcdcAcB));
1076      return F;
1077    }
1078  }
1079  if (fdivides (lcB, lcA))
1080  {
1081    if (fdivides (ppB, ppA, coF))
1082    {
1083      coG= 1;
1084      if (compressConvexDense)
1085        coF= decompress (coF, MM, V);
1086      coF= N (coF*(cA/gcdcAcB));
1087      return G;
1088    }
1089  }*/
1090
1091  gcdlcAlcB= gcd (lcA, lcB);
1092
1093  int d= totaldegree (ppA, Variable(2), Variable (ppA.level()));
1094  if (d == 0)
1095  {
1096    coF= N (ppA*(cA/gcdcAcB));
1097    coG= N (ppB*(cB/gcdcAcB));
1098    return N(gcdcAcB);
1099  }
1100  int d0= totaldegree (ppB, Variable(2), Variable (ppB.level()));
1101  if (d0 < d)
1102    d= d0;
1103  if (d == 0)
1104  {
1105    coF= N (ppA*(cA/gcdcAcB));
1106    coG= N (ppB*(cB/gcdcAcB));
1107    return N(gcdcAcB);
1108  }
1109
1110  CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH;
1111  CanonicalForm newtonPoly, coF_random_element, coG_random_element, coF_m,
1112                coG_m, ppCoF, ppCoG;
1113
1114  newtonPoly= 1;
1115  m= gcdlcAlcB;
1116  G_m= 0;
1117  coF= 0;
1118  coG= 0;
1119  H= 0;
1120  bool fail= false;
1121  //topLevel= false;
1122  bool inextension= false;
1123  int p=-1;
1124  int k= getGFDegree();
1125  int kk;
1126  int expon;
1127  char gf_name_buf= gf_name;
1128  int bound1= degree (ppA, 1);
1129  int bound2= degree (ppB, 1);
1130  do
1131  {
1132    random_element= GFRandomElement (m*lcA*lcB, l, fail);
1133    if (fail)
1134    {
1135      p= getCharacteristic();
1136      expon= 2;
1137      kk= getGFDegree();
1138      if (ipower (p, kk*expon) < (1 << 16))
1139        setCharacteristic (p, kk*(int)expon, 'b');
1140      else
1141      {
1142        expon= (int) floor((log ((double)((1<<16) - 1)))/(log((double)p)*kk));
1143        ASSERT (expon >= 2, "not enough points in GCD_GF");
1144        setCharacteristic (p, (int)(kk*expon), 'b');
1145      }
1146      inextension= true;
1147      fail= false;
1148      for (CFListIterator i= l; i.hasItem(); i++)
1149        i.getItem()= GFMapUp (i.getItem(), kk);
1150      m= GFMapUp (m, kk);
1151      G_m= GFMapUp (G_m, kk);
1152      newtonPoly= GFMapUp (newtonPoly, kk);
1153      coF_m= GFMapUp (coF_m, kk);
1154      coG_m= GFMapUp (coG_m, kk);
1155      ppA= GFMapUp (ppA, kk);
1156      ppB= GFMapUp (ppB, kk);
1157      gcdlcAlcB= GFMapUp (gcdlcAlcB, kk);
1158      lcA= GFMapUp (lcA, kk);
1159      lcB= GFMapUp (lcB, kk);
1160      random_element= GFRandomElement (m*lcA*lcB, l, fail);
1161      DEBOUTLN (cerr, "fail= " << fail);
1162      CFList list;
1163      TIMING_START (gcd_recursion);
1164      G_random_element= GCD_GF (ppA(random_element, x), ppB(random_element, x),
1165                                coF_random_element, coG_random_element,
1166                                list, topLevel);
1167      TIMING_END_AND_PRINT (gcd_recursion,
1168                            "time for recursive call: ");
1169      DEBOUTLN (cerr, "G_random_element= " << G_random_element);
1170    }
1171    else
1172    {
1173      CFList list;
1174      TIMING_START (gcd_recursion);
1175      G_random_element= GCD_GF (ppA(random_element, x), ppB(random_element, x),
1176                                coF_random_element, coG_random_element,
1177                                list, topLevel);
1178      TIMING_END_AND_PRINT (gcd_recursion,
1179                            "time for recursive call: ");
1180      DEBOUTLN (cerr, "G_random_element= " << G_random_element);
1181    }
1182
1183    if (!G_random_element.inCoeffDomain())
1184      d0= totaldegree (G_random_element, Variable(2),
1185                       Variable (G_random_element.level()));
1186    else
1187      d0= 0;
1188
1189    if (d0 == 0)
1190    {
1191      if (inextension)
1192      {
1193        ppA= GFMapDown (ppA, k);
1194        ppB= GFMapDown (ppB, k);
1195        setCharacteristic (p, k, gf_name_buf);
1196      }
1197      coF= N (ppA*(cA/gcdcAcB));
1198      coG= N (ppB*(cB/gcdcAcB));
1199      return N(gcdcAcB);
1200    }
1201    if (d0 >  d)
1202    {
1203      if (!find (l, random_element))
1204        l.append (random_element);
1205      continue;
1206    }
1207
1208    G_random_element=
1209    (gcdlcAlcB(random_element, x)/uni_lcoeff(G_random_element)) *
1210     G_random_element;
1211
1212    coF_random_element= (lcA(random_element,x)/uni_lcoeff(coF_random_element))
1213                        *coF_random_element;
1214    coG_random_element= (lcB(random_element,x)/uni_lcoeff(coG_random_element))
1215                        *coG_random_element;
1216
1217    if (!G_random_element.inCoeffDomain())
1218      d0= totaldegree (G_random_element, Variable(2),
1219                       Variable (G_random_element.level()));
1220    else
1221      d0= 0;
1222
1223    if (d0 < d)
1224    {
1225      m= gcdlcAlcB;
1226      newtonPoly= 1;
1227      G_m= 0;
1228      d= d0;
1229      coF_m= 0;
1230      coG_m= 0;
1231    }
1232
1233    TIMING_START (newton_interpolation);
1234    H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
1235    coF= newtonInterp (random_element, coF_random_element, newtonPoly, coF_m,x);
1236    coG= newtonInterp (random_element, coG_random_element, newtonPoly, coG_m,x);
1237    TIMING_END_AND_PRINT (newton_interpolation,
1238                          "time for newton interpolation: ");
1239
1240    //termination test
1241    if ((uni_lcoeff (H) == gcdlcAlcB) || (G_m == H))
1242    {
1243      TIMING_START (termination_test);
1244      if (gcdlcAlcB.isOne())
1245        cH= 1;
1246      else
1247        cH= uni_content (H);
1248      ppH= H/cH;
1249      ppH /= Lc (ppH);
1250      CanonicalForm lcppH= gcdlcAlcB/cH;
1251      CanonicalForm ccoF= lcppH/Lc (lcppH);
1252      CanonicalForm ccoG= lcppH/Lc (lcppH);
1253      ppCoF= coF/ccoF;
1254      ppCoG= coG/ccoG;
1255      if (inextension)
1256      {
1257        if (((degree (ppCoF,1)+degree (ppH,1) == bound1) &&
1258             (degree (ppCoG,1)+degree (ppH,1) == bound2) &&
1259             terminationTest (ppA, ppB, ppCoF, ppCoG, ppH)) ||
1260             (fdivides (ppH, ppA, ppCoF) && fdivides (ppH, ppB, ppCoG)))
1261        {
1262          DEBOUTLN (cerr, "ppH before mapDown= " << ppH);
1263          ppH= GFMapDown (ppH, k);
1264          ppCoF= GFMapDown (ppCoF, k);
1265          ppCoG= GFMapDown (ppCoG, k);
1266          DEBOUTLN (cerr, "ppH after mapDown= " << ppH);
1267          if (compressConvexDense)
1268          {
1269            ppH= decompress (ppH, MM, V);
1270            ppCoF= decompress (ppCoF, MM, V);
1271            ppCoG= decompress (ppCoG, MM, V);
1272          }
1273          coF= N ((cA/gcdcAcB)*ppCoF);
1274          coG= N ((cB/gcdcAcB)*ppCoG);
1275          setCharacteristic (p, k, gf_name_buf);
1276          TIMING_END_AND_PRINT (termination_test,
1277                                "time for successful termination GF: ");
1278          return N(gcdcAcB*ppH);
1279        }
1280      }
1281      else
1282      {
1283      if (((degree (ppCoF,1)+degree (ppH,1) == bound1) &&
1284           (degree (ppCoG,1)+degree (ppH,1) == bound2) &&
1285           terminationTest (ppA, ppB, ppCoF, ppCoG, ppH)) ||
1286           (fdivides (ppH, ppA, ppCoF) && fdivides (ppH, ppB, ppCoG)))
1287        {
1288          if (compressConvexDense)
1289          {
1290            ppH= decompress (ppH, MM, V);
1291            ppCoF= decompress (ppCoF, MM, V);
1292            ppCoG= decompress (ppCoG, MM, V);
1293          }
1294          coF= N ((cA/gcdcAcB)*ppCoF);
1295          coG= N ((cB/gcdcAcB)*ppCoG);
1296          TIMING_END_AND_PRINT (termination_test,
1297                                "time for successful termination GF: ");
1298          return N(gcdcAcB*ppH);
1299        }
1300      }
1301      TIMING_END_AND_PRINT (termination_test,
1302                            "time for unsuccessful termination GF: ");
1303    }
1304
1305    G_m= H;
1306    coF_m= coF;
1307    coG_m= coG;
1308    newtonPoly= newtonPoly*(x - random_element);
1309    m= m*(x - random_element);
1310    if (!find (l, random_element))
1311      l.append (random_element);
1312  } while (1);
1313}
1314
1315/// F is assumed to be an univariate polynomial in x,
1316/// computes a random monic irreducible univariate polynomial of random
1317/// degree < i in x which does not divide F
1318CanonicalForm
1319randomIrredpoly (int i, const Variable & x)
1320{
1321  int p= getCharacteristic();
1322  if (fac_NTL_char != p)
1323  {
1324    fac_NTL_char= p;
1325    zz_p::init (p);
1326  }
1327  zz_pX NTLirredpoly;
1328  CanonicalForm CFirredpoly;
1329  BuildIrred (NTLirredpoly, i + 1);
1330  CFirredpoly= convertNTLzzpX2CF (NTLirredpoly, x);
1331  return CFirredpoly;
1332}
1333
1334static inline
1335CanonicalForm
1336FpRandomElement (const CanonicalForm& F, CFList& list, bool& fail)
1337{
1338  fail= false;
1339  Variable x= F.mvar();
1340  FFRandom genFF;
1341  CanonicalForm random;
1342  int p= getCharacteristic();
1343  int bound= p;
1344  do
1345  {
1346    if (list.length() == bound)
1347    {
1348      fail= true;
1349      break;
1350    }
1351    if (list.length() < 1)
1352      random= 0;
1353    else
1354    {
1355      random= genFF.generate();
1356      while (find (list, random))
1357        random= genFF.generate();
1358    }
1359    if (F (random, x) == 0)
1360    {
1361      list.append (random);
1362      continue;
1363    }
1364  } while (find (list, random));
1365  return random;
1366}
1367
1368CanonicalForm
1369GCD_small_p (const CanonicalForm& F, const CanonicalForm&  G,
1370             CanonicalForm& coF, CanonicalForm& coG,
1371             bool& topLevel, CFList& l);
1372
1373CanonicalForm
1374GCD_small_p (const CanonicalForm& F, const CanonicalForm& G,
1375             bool& topLevel, CFList& l)
1376{
1377  CanonicalForm dummy1, dummy2;
1378  CanonicalForm result= GCD_small_p (F, G, dummy1, dummy2, topLevel, l);
1379  return result;
1380}
1381
1382CanonicalForm
1383GCD_small_p (const CanonicalForm& F, const CanonicalForm&  G,
1384             CanonicalForm& coF, CanonicalForm& coG,
1385             bool& topLevel, CFList& l)
1386{
1387  CanonicalForm A= F;
1388  CanonicalForm B= G;
1389  if (F.isZero() && degree(G) > 0)
1390  {
1391    coF= 0;
1392    coG= Lc (G);
1393    return G/Lc(G);
1394  }
1395  else if (G.isZero() && degree (F) > 0)
1396  {
1397    coF= Lc (F);
1398    coG= 0;
1399    return F/Lc(F);
1400  }
1401  else if (F.isZero() && G.isZero())
1402  {
1403    coF= coG= 0;
1404    return F.genOne();
1405  }
1406  if (F.inBaseDomain() || G.inBaseDomain())
1407  {
1408    coF= F;
1409    coG= G;
1410    return F.genOne();
1411  }
1412  if (F.isUnivariate() && fdivides(F, G, coG))
1413  {
1414    coF= Lc (F);
1415    return F/Lc(F);
1416  }
1417  if (G.isUnivariate() && fdivides(G, F, coF))
1418  {
1419    coG= Lc (G);
1420    return G/Lc(G);
1421  }
1422  if (F == G)
1423  {
1424    coF= coG= Lc (F);
1425    return F/Lc(F);
1426  }
1427  CFMap M,N;
1428  int best_level= myCompress (A, B, M, N, topLevel);
1429
1430  if (best_level == 0)
1431  {
1432    coF= F;
1433    coG= G;
1434    return B.genOne();
1435  }
1436
1437  A= M(A);
1438  B= M(B);
1439
1440  Variable x= Variable (1);
1441
1442  //univariate case
1443  if (A.isUnivariate() && B.isUnivariate())
1444  {
1445    CanonicalForm result= gcd (A, B);
1446    coF= N (A/result);
1447    coG= N (B/result);
1448    return N (result);
1449  }
1450
1451  CanonicalForm cA, cB;    // content of A and B
1452  CanonicalForm ppA, ppB;    // primitive part of A and B
1453  CanonicalForm gcdcAcB;
1454
1455  cA = uni_content (A);
1456  cB = uni_content (B);
1457  gcdcAcB= gcd (cA, cB);
1458  ppA= A/cA;
1459  ppB= B/cB;
1460
1461  int sizeNewtonPolyg;
1462  int ** newtonPolyg= NULL;
1463  mat_ZZ MM;
1464  vec_ZZ V;
1465  bool compressConvexDense= false; //(ppA.level() == 2 && ppB.level() == 2);
1466  if (compressConvexDense)
1467  {
1468    CanonicalForm bufcA= cA;
1469    CanonicalForm bufcB= cB;
1470    cA= content (ppA, 1);
1471    cB= content (ppB, 1);
1472    ppA /= cA;
1473    ppB /= cB;
1474    gcdcAcB *= gcd (cA, cB);
1475    cA *= bufcA;
1476    cB *= bufcB;
1477    if (ppA.isUnivariate() || ppB.isUnivariate())
1478    {
1479      if (ppA.level() == ppB.level())
1480      {
1481        CanonicalForm result= gcd (ppA, ppB);
1482        coF= N ((ppA/result)*(cA/gcdcAcB));
1483        coG= N ((ppB/result)*(cB/gcdcAcB));
1484        return N (result*gcdcAcB);
1485      }
1486      else
1487      {
1488        coF= N (ppA*(cA/gcdcAcB));
1489        coG= N (ppB*(cB/gcdcAcB));
1490        return N (gcdcAcB);
1491      }
1492    }
1493
1494    newtonPolyg= newtonPolygon (ppA,ppB, sizeNewtonPolyg);
1495    convexDense (newtonPolyg, sizeNewtonPolyg, MM, V);
1496
1497    for (int i= 0; i < sizeNewtonPolyg; i++)
1498      delete [] newtonPolyg[i];
1499    delete [] newtonPolyg;
1500
1501    ppA= compress (ppA, MM, V, false);
1502    ppB= compress (ppB, MM, V, false);
1503    MM= inv (MM);
1504
1505    if (ppA.isUnivariate() && ppB.isUnivariate())
1506    {
1507      if (ppA.level() == ppB.level())
1508      {
1509        CanonicalForm result= gcd (ppA, ppB);
1510        coF= N (decompress ((ppA/result), MM, V)*(cA/gcdcAcB));
1511        coG= N (decompress ((ppB/result), MM, V)*(cB/gcdcAcB));
1512        return N (decompress (result, MM, V)*gcdcAcB);
1513      }
1514      else
1515      {
1516        coF= N (decompress (ppA, MM, V));
1517        coG= N (decompress (ppB, MM, V));
1518        return N (gcdcAcB);
1519      }
1520    }
1521  }
1522
1523  CanonicalForm lcA, lcB;  // leading coefficients of A and B
1524  CanonicalForm gcdlcAlcB;
1525  lcA= uni_lcoeff (ppA);
1526  lcB= uni_lcoeff (ppB);
1527
1528  /*if (fdivides (lcA, lcB))
1529  {
1530    if (fdivides (A, B))
1531      return F/Lc(F);
1532  }
1533  if (fdivides (lcB, lcA))
1534  {
1535    if (fdivides (B, A))
1536      return G/Lc(G);
1537  }*/
1538
1539  gcdlcAlcB= gcd (lcA, lcB);
1540
1541  int d= totaldegree (ppA, Variable (2), Variable (ppA.level()));
1542  int d0;
1543
1544  if (d == 0)
1545  {
1546    coF= N (ppA*(cA/gcdcAcB));
1547    coG= N (ppB*(cB/gcdcAcB));
1548    return N(gcdcAcB);
1549  }
1550
1551  d0= totaldegree (ppB, Variable (2), Variable (ppB.level()));
1552
1553  if (d0 < d)
1554    d= d0;
1555
1556  if (d == 0)
1557  {
1558    coF= N (ppA*(cA/gcdcAcB));
1559    coG= N (ppB*(cB/gcdcAcB));
1560    return N(gcdcAcB);
1561  }
1562
1563  CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH;
1564  CanonicalForm newtonPoly, coF_random_element, coG_random_element,
1565                coF_m, coG_m, ppCoF, ppCoG;
1566
1567  newtonPoly= 1;
1568  m= gcdlcAlcB;
1569  H= 0;
1570  coF= 0;
1571  coG= 0;
1572  G_m= 0;
1573  Variable alpha, V_buf;
1574  bool fail= false;
1575  bool inextension= false;
1576  topLevel= false;
1577  CFList source, dest;
1578  int bound1= degree (ppA, 1);
1579  int bound2= degree (ppB, 1);
1580  do
1581  {
1582    if (inextension)
1583      random_element= randomElement (m*lcA*lcB, V_buf, l, fail);
1584    else
1585      random_element= FpRandomElement (m*lcA*lcB, l, fail);
1586
1587    if (!fail && !inextension)
1588    {
1589      CFList list;
1590      TIMING_START (gcd_recursion);
1591      G_random_element=
1592      GCD_small_p (ppA (random_element,x), ppB (random_element,x),
1593                   coF_random_element, coG_random_element, topLevel,
1594                   list);
1595      TIMING_END_AND_PRINT (gcd_recursion,
1596                            "time for recursive call: ");
1597      DEBOUTLN (cerr, "G_random_element= " << G_random_element);
1598    }
1599    else if (!fail && inextension)
1600    {
1601      CFList list;
1602      TIMING_START (gcd_recursion);
1603      G_random_element=
1604      GCD_Fp_extension (ppA (random_element, x), ppB (random_element, x),
1605                        coF_random_element, coG_random_element, V_buf,
1606                        list, topLevel);
1607      TIMING_END_AND_PRINT (gcd_recursion,
1608                            "time for recursive call: ");
1609      DEBOUTLN (cerr, "G_random_element= " << G_random_element);
1610    }
1611    else if (fail && !inextension)
1612    {
1613      source= CFList();
1614      dest= CFList();
1615      CFList list;
1616      CanonicalForm mipo;
1617      int deg= 2;
1618      do {
1619        mipo= randomIrredpoly (deg, x);
1620        alpha= rootOf (mipo);
1621        inextension= true;
1622        fail= false;
1623        random_element= randomElement (m*lcA*lcB, alpha, l, fail);
1624        deg++;
1625      } while (fail);
1626      list= CFList();
1627      V_buf= alpha;
1628      TIMING_START (gcd_recursion);
1629      G_random_element=
1630      GCD_Fp_extension (ppA (random_element, x), ppB (random_element, x),
1631                        coF_random_element, coG_random_element, alpha,
1632                        list, topLevel);
1633      TIMING_END_AND_PRINT (gcd_recursion,
1634                            "time for recursive call: ");
1635      DEBOUTLN (cerr, "G_random_element= " << G_random_element);
1636    }
1637    else if (fail && inextension)
1638    {
1639      source= CFList();
1640      dest= CFList();
1641
1642      Variable V_buf3= V_buf;
1643      V_buf= chooseExtension (V_buf);
1644      bool prim_fail= false;
1645      Variable V_buf2;
1646      CanonicalForm prim_elem, im_prim_elem;
1647      prim_elem= primitiveElement (alpha, V_buf2, prim_fail);
1648
1649      if (V_buf3 != alpha)
1650      {
1651        m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
1652        G_m= mapDown (G_m, prim_elem, im_prim_elem, alpha, source, dest);
1653        coF_m= mapDown (coF_m, prim_elem, im_prim_elem, alpha, source, dest);
1654        coG_m= mapDown (coG_m, prim_elem, im_prim_elem, alpha, source, dest);
1655        newtonPoly= mapDown (newtonPoly, prim_elem, im_prim_elem, alpha,
1656                             source, dest);
1657        ppA= mapDown (ppA, prim_elem, im_prim_elem, alpha, source, dest);
1658        ppB= mapDown (ppB, prim_elem, im_prim_elem, alpha, source, dest);
1659        gcdlcAlcB= mapDown (gcdlcAlcB, prim_elem, im_prim_elem, alpha, source,
1660                            dest);
1661        lcA= mapDown (lcA, prim_elem, im_prim_elem, alpha, source, dest);
1662        lcB= mapDown (lcB, prim_elem, im_prim_elem, alpha, source, dest);
1663        for (CFListIterator i= l; i.hasItem(); i++)
1664          i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, alpha,
1665                                source, dest);
1666      }
1667
1668      ASSERT (!prim_fail, "failure in integer factorizer");
1669      if (prim_fail)
1670        ; //ERROR
1671      else
1672        im_prim_elem= mapPrimElem (prim_elem, alpha, V_buf);
1673
1674      DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (alpha));
1675      DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf2));
1676
1677      for (CFListIterator i= l; i.hasItem(); i++)
1678        i.getItem()= mapUp (i.getItem(), alpha, V_buf, prim_elem,
1679                             im_prim_elem, source, dest);
1680      m= mapUp (m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
1681      G_m= mapUp (G_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
1682      coF_m= mapUp (coF_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
1683      coG_m= mapUp (coG_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
1684      newtonPoly= mapUp (newtonPoly, alpha, V_buf, prim_elem, im_prim_elem,
1685                          source, dest);
1686      ppA= mapUp (ppA, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
1687      ppB= mapUp (ppB, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
1688      gcdlcAlcB= mapUp (gcdlcAlcB, alpha, V_buf, prim_elem, im_prim_elem,
1689                         source, dest);
1690      lcA= mapUp (lcA, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
1691      lcB= mapUp (lcB, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
1692      fail= false;
1693      random_element= randomElement (m*lcA*lcB, V_buf, l, fail );
1694      DEBOUTLN (cerr, "fail= " << fail);
1695      CFList list;
1696      TIMING_START (gcd_recursion);
1697      G_random_element=
1698      GCD_Fp_extension (ppA (random_element, x), ppB (random_element, x),
1699                        coF_random_element, coG_random_element, V_buf,
1700                        list, topLevel);
1701      TIMING_END_AND_PRINT (gcd_recursion,
1702                            "time for recursive call: ");
1703      DEBOUTLN (cerr, "G_random_element= " << G_random_element);
1704    }
1705
1706    if (!G_random_element.inCoeffDomain())
1707      d0= totaldegree (G_random_element, Variable(2),
1708                       Variable (G_random_element.level()));
1709    else
1710      d0= 0;
1711
1712    if (d0 == 0)
1713    {
1714      coF= N (ppA*(cA/gcdcAcB));
1715      coG= N (ppB*(cB/gcdcAcB));
1716      return N(gcdcAcB);
1717    }
1718
1719    if (d0 >  d)
1720    {
1721      if (!find (l, random_element))
1722        l.append (random_element);
1723      continue;
1724    }
1725
1726    G_random_element= (gcdlcAlcB(random_element,x)/uni_lcoeff(G_random_element))
1727                       *G_random_element;
1728
1729    coF_random_element= (lcA(random_element,x)/uni_lcoeff(coF_random_element))
1730                        *coF_random_element;
1731    coG_random_element= (lcB(random_element,x)/uni_lcoeff(coG_random_element))
1732                        *coG_random_element;
1733
1734    if (!G_random_element.inCoeffDomain())
1735      d0= totaldegree (G_random_element, Variable(2),
1736                       Variable (G_random_element.level()));
1737    else
1738      d0= 0;
1739
1740    if (d0 <  d)
1741    {
1742      m= gcdlcAlcB;
1743      newtonPoly= 1;
1744      G_m= 0;
1745      d= d0;
1746      coF_m= 0;
1747      coG_m= 0;
1748    }
1749
1750    TIMING_START (newton_interpolation);
1751    H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
1752    coF= newtonInterp (random_element, coF_random_element, newtonPoly, coF_m,x);
1753    coG= newtonInterp (random_element, coG_random_element, newtonPoly, coG_m,x);
1754    TIMING_END_AND_PRINT (newton_interpolation,
1755                          "time for newton_interpolation: ");
1756
1757    //termination test
1758    if ((uni_lcoeff (H) == gcdlcAlcB) || (G_m == H))
1759    {
1760      TIMING_START (termination_test);
1761      if (gcdlcAlcB.isOne())
1762        cH= 1;
1763      else
1764        cH= uni_content (H);
1765      ppH= H/cH;
1766      ppH /= Lc (ppH);
1767      CanonicalForm lcppH= gcdlcAlcB/cH;
1768      CanonicalForm ccoF= lcppH/Lc (lcppH);
1769      CanonicalForm ccoG= lcppH/Lc (lcppH);
1770      ppCoF= coF/ccoF;
1771      ppCoG= coG/ccoG;
1772      DEBOUTLN (cerr, "ppH= " << ppH);
1773      if (((degree (ppCoF,1)+degree (ppH,1) == bound1) &&
1774           (degree (ppCoG,1)+degree (ppH,1) == bound2) &&
1775           terminationTest (ppA, ppB, ppCoF, ppCoG, ppH)) ||
1776           (fdivides (ppH, ppA, ppCoF) && fdivides (ppH, ppB, ppCoG)))
1777      {
1778        if (compressConvexDense)
1779        {
1780          ppH= decompress (ppH, MM, V);
1781          ppCoF= decompress (ppCoF, MM, V);
1782          ppCoG= decompress (ppCoG, MM, V);
1783        }
1784        coF= N ((cA/gcdcAcB)*ppCoF);
1785        coG= N ((cB/gcdcAcB)*ppCoG);
1786        TIMING_END_AND_PRINT (termination_test,
1787                              "time for successful termination Fp: ");
1788        return N(gcdcAcB*ppH);
1789      }
1790      TIMING_END_AND_PRINT (termination_test,
1791                            "time for unsuccessful termination Fp: ");
1792    }
1793
1794    G_m= H;
1795    coF_m= coF;
1796    coG_m= coG;
1797    newtonPoly= newtonPoly*(x - random_element);
1798    m= m*(x - random_element);
1799    if (!find (l, random_element))
1800      l.append (random_element);
1801  } while (1);
1802}
1803
1804CFArray
1805solveVandermonde (const CFArray& M, const CFArray& A)
1806{
1807  int r= M.size();
1808  ASSERT (A.size() == r, "vector does not have right size");
1809
1810  if (r == 1)
1811  {
1812    CFArray result= CFArray (1);
1813    result [0]= A [0] / M [0];
1814    return result;
1815  }
1816  // check solvability
1817  bool notDistinct= false;
1818  for (int i= 0; i < r - 1; i++)
1819  {
1820    for (int j= i + 1; j < r; j++)
1821    {
1822      if (M [i] == M [j])
1823      {
1824        notDistinct= true;
1825        break;
1826      }
1827    }
1828  }
1829  if (notDistinct)
1830    return CFArray();
1831
1832  CanonicalForm master= 1;
1833  Variable x= Variable (1);
1834  for (int i= 0; i < r; i++)
1835    master *= x - M [i];
1836  CFList Pj;
1837  CanonicalForm tmp;
1838  for (int i= 0; i < r; i++)
1839  {
1840    tmp= master/(x - M [i]);
1841    tmp /= tmp (M [i], 1);
1842    Pj.append (tmp);
1843  }
1844  CFArray result= CFArray (r);
1845
1846  CFListIterator j= Pj;
1847  for (int i= 1; i <= r; i++, j++)
1848  {
1849    tmp= 0;
1850    for (int l= 0; l < A.size(); l++)
1851      tmp += A[l]*j.getItem()[l];
1852    result[i - 1]= tmp;
1853  }
1854  return result;
1855}
1856
1857CFArray
1858solveGeneralVandermonde (const CFArray& M, const CFArray& A)
1859{
1860  int r= M.size();
1861  ASSERT (A.size() == r, "vector does not have right size");
1862  if (r == 1)
1863  {
1864    CFArray result= CFArray (1);
1865    result [0]= A[0] / M [0];
1866    return result;
1867  }
1868  // check solvability
1869  bool notDistinct= false;
1870  for (int i= 0; i < r - 1; i++)
1871  {
1872    for (int j= i + 1; j < r; j++)
1873    {
1874      if (M [i] == M [j])
1875      {
1876        notDistinct= true;
1877        break;
1878      }
1879    }
1880  }
1881  if (notDistinct)
1882    return CFArray();
1883
1884  CanonicalForm master= 1;
1885  Variable x= Variable (1);
1886  for (int i= 0; i < r; i++)
1887    master *= x - M [i];
1888  master *= x;
1889  CFList Pj;
1890  CanonicalForm tmp;
1891  for (int i= 0; i < r; i++)
1892  {
1893    tmp= master/(x - M [i]);
1894    tmp /= tmp (M [i], 1);
1895    Pj.append (tmp);
1896  }
1897
1898  CFArray result= CFArray (r);
1899
1900  CFListIterator j= Pj;
1901  for (int i= 1; i <= r; i++, j++)
1902  {
1903    tmp= 0;
1904
1905    for (int l= 1; l <= A.size(); l++)
1906      tmp += A[l - 1]*j.getItem()[l];
1907    result[i - 1]= tmp;
1908  }
1909  return result;
1910}
1911
1912/// M in row echolon form, rk rank of M
1913CFArray
1914readOffSolution (const CFMatrix& M, const long rk)
1915{
1916  CFArray result= CFArray (rk);
1917  CanonicalForm tmp1, tmp2, tmp3;
1918  for (int i= rk; i >= 1; i--)
1919  {
1920    tmp3= 0;
1921    tmp1= M (i, M.columns());
1922    for (int j= M.columns() - 1; j >= 1; j--)
1923    {
1924      tmp2= M (i, j);
1925      if (j == i)
1926        break;
1927      else
1928        tmp3 += tmp2*result[j - 1];
1929    }
1930    result[i - 1]= (tmp1 - tmp3)/tmp2;
1931  }
1932  return result;
1933}
1934
1935CFArray
1936readOffSolution (const CFMatrix& M, const CFArray& L, const CFArray& partialSol)
1937{
1938  CFArray result= CFArray (M.rows());
1939  CanonicalForm tmp1, tmp2, tmp3;
1940  int k;
1941  for (int i= M.rows(); i >= 1; i--)
1942  {
1943    tmp3= 0;
1944    tmp1= L[i - 1];
1945    k= 0;
1946    for (int j= M.columns(); j >= 1; j--, k++)
1947    {
1948      tmp2= M (i, j);
1949      if (j == i)
1950        break;
1951      else
1952      {
1953        if (k > partialSol.size() - 1)
1954          tmp3 += tmp2*result[j - 1];
1955        else
1956          tmp3 += tmp2*partialSol[partialSol.size() - k - 1];
1957      }
1958    }
1959    result [i - 1]= (tmp1 - tmp3)/tmp2;
1960  }
1961  return result;
1962}
1963
1964long
1965gaussianElimFp (CFMatrix& M, CFArray& L)
1966{
1967  ASSERT (L.size() <= M.rows(), "dimension exceeded");
1968  CFMatrix *N;
1969  N= new CFMatrix (M.rows(), M.columns() + 1);
1970
1971  for (int i= 1; i <= M.rows(); i++)
1972    for (int j= 1; j <= M.columns(); j++)
1973      (*N) (i, j)= M (i, j);
1974
1975  int j= 1;
1976  for (int i= 0; i < L.size(); i++, j++)
1977    (*N) (j, M.columns() + 1)= L[i];
1978#ifdef HAVE_FLINT
1979  nmod_mat_t FLINTN;
1980  convertFacCFMatrix2nmod_mat_t (FLINTN, *N);
1981  long rk= nmod_mat_rref (FLINTN);
1982
1983  N= convertNmod_mat_t2FacCFMatrix (FLINTN);
1984  nmod_mat_clear (FLINTN);
1985#else
1986  int p= getCharacteristic ();
1987  if (fac_NTL_char != p)
1988  {
1989    fac_NTL_char= p;
1990    zz_p::init (p);
1991  }
1992  mat_zz_p *NTLN= convertFacCFMatrix2NTLmat_zz_p(*N);
1993  long rk= gauss (*NTLN);
1994
1995  N= convertNTLmat_zz_p2FacCFMatrix (*NTLN);
1996#endif
1997
1998  L= CFArray (M.rows());
1999  for (int i= 0; i < M.rows(); i++)
2000    L[i]= (*N) (i + 1, M.columns() + 1);
2001  M= (*N) (1, M.rows(), 1, M.columns());
2002  delete N;
2003  return rk;
2004}
2005
2006long
2007gaussianElimFq (CFMatrix& M, CFArray& L, const Variable& alpha)
2008{
2009  ASSERT (L.size() <= M.rows(), "dimension exceeded");
2010  CFMatrix *N;
2011  N= new CFMatrix (M.rows(), M.columns() + 1);
2012
2013  for (int i= 1; i <= M.rows(); i++)
2014    for (int j= 1; j <= M.columns(); j++)
2015      (*N) (i, j)= M (i, j);
2016
2017  int j= 1;
2018  for (int i= 0; i < L.size(); i++, j++)
2019    (*N) (j, M.columns() + 1)= L[i];
2020  int p= getCharacteristic ();
2021  if (fac_NTL_char != p)
2022  {
2023    fac_NTL_char= p;
2024    zz_p::init (p);
2025  }
2026  zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
2027  zz_pE::init (NTLMipo);
2028  mat_zz_pE *NTLN= convertFacCFMatrix2NTLmat_zz_pE(*N);
2029  long rk= gauss (*NTLN);
2030
2031  N= convertNTLmat_zz_pE2FacCFMatrix (*NTLN, alpha);
2032
2033  M= (*N) (1, M.rows(), 1, M.columns());
2034  L= CFArray (M.rows());
2035  for (int i= 0; i < M.rows(); i++)
2036    L[i]= (*N) (i + 1, M.columns() + 1);
2037
2038  delete N;
2039  return rk;
2040}
2041
2042CFArray
2043solveSystemFp (const CFMatrix& M, const CFArray& L)
2044{
2045  ASSERT (L.size() <= M.rows(), "dimension exceeded");
2046  CFMatrix *N;
2047  N= new CFMatrix (M.rows(), M.columns() + 1);
2048
2049  for (int i= 1; i <= M.rows(); i++)
2050    for (int j= 1; j <= M.columns(); j++)
2051      (*N) (i, j)= M (i, j);
2052
2053  int j= 1;
2054  for (int i= 0; i < L.size(); i++, j++)
2055    (*N) (j, M.columns() + 1)= L[i];
2056
2057#ifdef HAVE_FLINT
2058  nmod_mat_t FLINTN;
2059  convertFacCFMatrix2nmod_mat_t (FLINTN, *N);
2060  long rk= nmod_mat_rref (FLINTN);
2061#else
2062  int p= getCharacteristic ();
2063  if (fac_NTL_char != p)
2064  {
2065    fac_NTL_char= p;
2066    zz_p::init (p);
2067  }
2068  mat_zz_p *NTLN= convertFacCFMatrix2NTLmat_zz_p(*N);
2069  long rk= gauss (*NTLN);
2070#endif
2071  if (rk != M.columns())
2072  {
2073#ifdef HAVE_FLINT
2074    nmod_mat_clear (FLINTN);
2075#endif
2076    delete N;
2077    return CFArray();
2078  }
2079#ifdef HAVE_FLINT
2080  N= convertNmod_mat_t2FacCFMatrix (FLINTN);
2081  nmod_mat_clear (FLINTN);
2082#else
2083  N= convertNTLmat_zz_p2FacCFMatrix (*NTLN);
2084#endif
2085  CFArray A= readOffSolution (*N, rk);
2086
2087  delete N;
2088  return A;
2089}
2090
2091CFArray
2092solveSystemFq (const CFMatrix& M, const CFArray& L, const Variable& alpha)
2093{
2094  ASSERT (L.size() <= M.rows(), "dimension exceeded");
2095  CFMatrix *N;
2096  N= new CFMatrix (M.rows(), M.columns() + 1);
2097
2098  for (int i= 1; i <= M.rows(); i++)
2099    for (int j= 1; j <= M.columns(); j++)
2100      (*N) (i, j)= M (i, j);
2101  int j= 1;
2102  for (int i= 0; i < L.size(); i++, j++)
2103    (*N) (j, M.columns() + 1)= L[i];
2104  int p= getCharacteristic ();
2105  if (fac_NTL_char != p)
2106  {
2107    fac_NTL_char= p;
2108    zz_p::init (p);
2109  }
2110  zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
2111  zz_pE::init (NTLMipo);
2112  mat_zz_pE *NTLN= convertFacCFMatrix2NTLmat_zz_pE(*N);
2113  long rk= gauss (*NTLN);
2114  if (rk != M.columns())
2115  {
2116    delete N;
2117    return CFArray();
2118  }
2119  N= convertNTLmat_zz_pE2FacCFMatrix (*NTLN, alpha);
2120
2121  CFArray A= readOffSolution (*N, rk);
2122
2123  delete N;
2124  return A;
2125}
2126#endif
2127
2128CFArray
2129getMonoms (const CanonicalForm& F)
2130{
2131  if (F.inCoeffDomain())
2132  {
2133    CFArray result= CFArray (1);
2134    result [0]= 1;
2135    return result;
2136  }
2137  if (F.isUnivariate())
2138  {
2139    CFArray result= CFArray (size(F));
2140    int j= 0;
2141    for (CFIterator i= F; i.hasTerms(); i++, j++)
2142      result[j]= power (F.mvar(), i.exp());
2143    return result;
2144  }
2145  int numMon= size (F);
2146  CFArray result= CFArray (numMon);
2147  int j= 0;
2148  CFArray recResult;
2149  Variable x= F.mvar();
2150  CanonicalForm powX;
2151  for (CFIterator i= F; i.hasTerms(); i++)
2152  {
2153    powX= power (x, i.exp());
2154    recResult= getMonoms (i.coeff());
2155    for (int k= 0; k < recResult.size(); k++)
2156      result[j+k]= powX*recResult[k];
2157    j += recResult.size();
2158  }
2159  return result;
2160}
2161
2162#ifdef HAVE_NTL
2163CFArray
2164evaluateMonom (const CanonicalForm& F, const CFList& evalPoints)
2165{
2166  if (F.inCoeffDomain())
2167  {
2168    CFArray result= CFArray (1);
2169    result [0]= F;
2170    return result;
2171  }
2172  if (F.isUnivariate())
2173  {
2174    ASSERT (evalPoints.length() == 1,
2175            "expected an eval point with only one component");
2176    CFArray result= CFArray (size(F));
2177    int j= 0;
2178    CanonicalForm evalPoint= evalPoints.getLast();
2179    for (CFIterator i= F; i.hasTerms(); i++, j++)
2180      result[j]= power (evalPoint, i.exp());
2181    return result;
2182  }
2183  int numMon= size (F);
2184  CFArray result= CFArray (numMon);
2185  int j= 0;
2186  CanonicalForm evalPoint= evalPoints.getLast();
2187  CFList buf= evalPoints;
2188  buf.removeLast();
2189  CFArray recResult;
2190  CanonicalForm powEvalPoint;
2191  for (CFIterator i= F; i.hasTerms(); i++)
2192  {
2193    powEvalPoint= power (evalPoint, i.exp());
2194    recResult= evaluateMonom (i.coeff(), buf);
2195    for (int k= 0; k < recResult.size(); k++)
2196      result[j+k]= powEvalPoint*recResult[k];
2197    j += recResult.size();
2198  }
2199  return result;
2200}
2201
2202CFArray
2203evaluate (const CFArray& A, const CFList& evalPoints)
2204{
2205  CFArray result= A.size();
2206  CanonicalForm tmp;
2207  int k;
2208  for (int i= 0; i < A.size(); i++)
2209  {
2210    tmp= A[i];
2211    k= 1;
2212    for (CFListIterator j= evalPoints; j.hasItem(); j++, k++)
2213      tmp= tmp (j.getItem(), k);
2214    result[i]= tmp;
2215  }
2216  return result;
2217}
2218
2219CFList
2220evaluationPoints (const CanonicalForm& F, const CanonicalForm& G,
2221                  CanonicalForm& Feval, CanonicalForm& Geval,
2222                  const CanonicalForm& LCF, const bool& GF,
2223                  const Variable& alpha, bool& fail, CFList& list
2224                 )
2225{
2226  int k= tmax (F.level(), G.level()) - 1;
2227  Variable x= Variable (1);
2228  CFList result;
2229  FFRandom genFF;
2230  GFRandom genGF;
2231  int p= getCharacteristic ();
2232  double bound;
2233  if (alpha != Variable (1))
2234  {
2235    bound= pow ((double) p, (double) degree (getMipo(alpha)));
2236    bound= pow (bound, (double) k);
2237  }
2238  else if (GF)
2239  {
2240    bound= pow ((double) p, (double) getGFDegree());
2241    bound= pow ((double) bound, (double) k);
2242  }
2243  else
2244    bound= pow ((double) p, (double) k);
2245
2246  CanonicalForm random;
2247  int j;
2248  bool zeroOneOccured= false;
2249  bool allEqual= false;
2250  CanonicalForm buf;
2251  do
2252  {
2253    random= 0;
2254    // possible overflow if list.length() does not fit into a int
2255    if (list.length() >= bound)
2256    {
2257      fail= true;
2258      break;
2259    }
2260    for (int i= 0; i < k; i++)
2261    {
2262      if (GF)
2263      {
2264        result.append (genGF.generate());
2265        random += result.getLast()*power (x, i);
2266      }
2267      else if (alpha.level() != 1)
2268      {
2269        AlgExtRandomF genAlgExt (alpha);
2270        result.append (genAlgExt.generate());
2271        random += result.getLast()*power (x, i);
2272      }
2273      else
2274      {
2275        result.append (genFF.generate());
2276        random += result.getLast()*power (x, i);
2277      }
2278      if (result.getLast().isOne() || result.getLast().isZero())
2279        zeroOneOccured= true;
2280    }
2281    if (find (list, random))
2282    {
2283      zeroOneOccured= false;
2284      allEqual= false;
2285      result= CFList();
2286      continue;
2287    }
2288    if (zeroOneOccured)
2289    {
2290      list.append (random);
2291      zeroOneOccured= false;
2292      allEqual= false;
2293      result= CFList();
2294      continue;
2295    }
2296    // no zero at this point
2297    if (k > 1)
2298    {
2299      allEqual= true;
2300      CFIterator iter= random;
2301      buf= iter.coeff();
2302      iter++;
2303      for (; iter.hasTerms(); iter++)
2304        if (buf != iter.coeff())
2305          allEqual= false;
2306    }
2307    if (allEqual)
2308    {
2309      list.append (random);
2310      allEqual= false;
2311      zeroOneOccured= false;
2312      result= CFList();
2313      continue;
2314    }
2315
2316    Feval= F;
2317    Geval= G;
2318    CanonicalForm LCeval= LCF;
2319    j= 1;
2320    for (CFListIterator i= result; i.hasItem(); i++, j++)
2321    {
2322      Feval= Feval (i.getItem(), j);
2323      Geval= Geval (i.getItem(), j);
2324      LCeval= LCeval (i.getItem(), j);
2325    }
2326
2327    if (LCeval.isZero())
2328    {
2329      if (!find (list, random))
2330        list.append (random);
2331      zeroOneOccured= false;
2332      allEqual= false;
2333      result= CFList();
2334      continue;
2335    }
2336
2337    if (list.length() >= bound)
2338    {
2339      fail= true;
2340      break;
2341    }
2342  } while (find (list, random));
2343
2344  return result;
2345}
2346
2347/// multiply two lists componentwise
2348void mult (CFList& L1, const CFList& L2)
2349{
2350  ASSERT (L1.length() == L2.length(), "lists of the same size expected");
2351
2352  CFListIterator j= L2;
2353  for (CFListIterator i= L1; i.hasItem(); i++, j++)
2354    i.getItem() *= j.getItem();
2355}
2356
2357void eval (const CanonicalForm& A, const CanonicalForm& B, CanonicalForm& Aeval,
2358           CanonicalForm& Beval, const CFList& L)
2359{
2360  Aeval= A;
2361  Beval= B;
2362  int j= 1;
2363  for (CFListIterator i= L; i.hasItem(); i++, j++)
2364  {
2365    Aeval= Aeval (i.getItem(), j);
2366    Beval= Beval (i.getItem(), j);
2367  }
2368}
2369
2370CanonicalForm
2371monicSparseInterpol (const CanonicalForm& F, const CanonicalForm& G,
2372                     const CanonicalForm& skeleton, const Variable& alpha,
2373                     bool& fail, CFArray*& coeffMonoms, CFArray& Monoms
2374                    )
2375{
2376  CanonicalForm A= F;
2377  CanonicalForm B= G;
2378  if (F.isZero() && degree(G) > 0) return G/Lc(G);
2379  else if (G.isZero() && degree (F) > 0) return F/Lc(F);
2380  else if (F.isZero() && G.isZero()) return F.genOne();
2381  if (F.inBaseDomain() || G.inBaseDomain()) return F.genOne();
2382  if (F.isUnivariate() && fdivides(F, G)) return F/Lc(F);
2383  if (G.isUnivariate() && fdivides(G, F)) return G/Lc(G);
2384  if (F == G) return F/Lc(F);
2385
2386  CFMap M,N;
2387  int best_level= myCompress (A, B, M, N, false);
2388
2389  if (best_level == 0)
2390    return B.genOne();
2391
2392  A= M(A);
2393  B= M(B);
2394
2395  Variable x= Variable (1);
2396  ASSERT (degree (A, x) == 0, "expected degree (F, 1) == 0");
2397  ASSERT (degree (B, x) == 0, "expected degree (G, 1) == 0");
2398
2399  //univariate case
2400  if (A.isUnivariate() && B.isUnivariate())
2401    return N (gcd (A, B));
2402
2403  CanonicalForm skel= M(skeleton);
2404  CanonicalForm cA, cB;    // content of A and B
2405  CanonicalForm ppA, ppB;    // primitive part of A and B
2406  CanonicalForm gcdcAcB;
2407  cA = uni_content (A);
2408  cB = uni_content (B);
2409  gcdcAcB= gcd (cA, cB);
2410  ppA= A/cA;
2411  ppB= B/cB;
2412
2413  CanonicalForm lcA, lcB;  // leading coefficients of A and B
2414  CanonicalForm gcdlcAlcB;
2415  lcA= uni_lcoeff (ppA);
2416  lcB= uni_lcoeff (ppB);
2417
2418  if (fdivides (lcA, lcB))
2419  {
2420    if (fdivides (A, B))
2421      return F/Lc(F);
2422  }
2423  if (fdivides (lcB, lcA))
2424  {
2425    if (fdivides (B, A))
2426      return G/Lc(G);
2427  }
2428
2429  gcdlcAlcB= gcd (lcA, lcB);
2430  int skelSize= size (skel, skel.mvar());
2431
2432  int j= 0;
2433  int biggestSize= 0;
2434
2435  for (CFIterator i= skel; i.hasTerms(); i++, j++)
2436    biggestSize= tmax (biggestSize, size (i.coeff()));
2437
2438  CanonicalForm g, Aeval, Beval;
2439
2440  CFList evalPoints;
2441  bool evalFail= false;
2442  CFList list;
2443  bool GF= false;
2444  CanonicalForm LCA= LC (A);
2445  CanonicalForm tmp;
2446  CFArray gcds= CFArray (biggestSize);
2447  CFList * pEvalPoints= new CFList [biggestSize];
2448  Variable V_buf= alpha;
2449  CFList source, dest;
2450  CanonicalForm prim_elem, im_prim_elem;
2451  for (int i= 0; i < biggestSize; i++)
2452  {
2453    if (i == 0)
2454      evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf, evalFail,
2455                                    list);
2456    else
2457    {
2458      mult (evalPoints, pEvalPoints [0]);
2459      eval (A, B, Aeval, Beval, evalPoints);
2460    }
2461
2462    if (evalFail)
2463    {
2464      if (V_buf.level() != 1)
2465      {
2466        do
2467        {
2468          Variable V_buf2= chooseExtension (V_buf);
2469          source= CFList();
2470          dest= CFList();
2471
2472          bool prim_fail= false;
2473          Variable V_buf3;
2474          prim_elem= primitiveElement (V_buf, V_buf3, prim_fail);
2475
2476          ASSERT (!prim_fail, "failure in integer factorizer");
2477          if (prim_fail)
2478            ; //ERROR
2479          else
2480            im_prim_elem= mapPrimElem (prim_elem, V_buf, V_buf2);
2481
2482          DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf));
2483          DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf2));
2484
2485          for (CFListIterator j= list; j.hasItem(); j++)
2486            j.getItem()= mapUp (j.getItem(), V_buf, V_buf2, prim_elem,
2487                                im_prim_elem, source, dest);
2488          for (int k= 0; k < i; k++)
2489          {
2490            for (CFListIterator j= pEvalPoints[k]; j.hasItem(); j++)
2491              j.getItem()= mapUp (j.getItem(), V_buf, V_buf2, prim_elem,
2492                                  im_prim_elem, source, dest);
2493            gcds[k]= mapUp (gcds[k], V_buf, V_buf2, prim_elem, im_prim_elem,
2494                            source, dest);
2495          }
2496
2497          if (alpha.level() != 1)
2498          {
2499            A= mapUp (A, V_buf, V_buf2, prim_elem, im_prim_elem, source,dest);
2500            B= mapUp (B, V_buf, V_buf2, prim_elem, im_prim_elem, source,dest);
2501          }
2502          evalFail= false;
2503          evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf,
2504                                        evalFail, list);
2505        } while (evalFail);
2506      }
2507      else
2508      {
2509        CanonicalForm mipo;
2510        int deg= 2;
2511        do {
2512          mipo= randomIrredpoly (deg, x);
2513          V_buf= rootOf (mipo);
2514          evalFail= false;
2515          evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf,
2516                                        evalFail, list);
2517          deg++;
2518        } while (evalFail);
2519      }
2520    }
2521
2522    g= gcd (Aeval, Beval);
2523    g /= Lc (g);
2524
2525    if (degree (g) != degree (skel) || (skelSize != size (g)))
2526    {
2527      delete[] pEvalPoints;
2528      fail= true;
2529      return 0;
2530    }
2531    CFIterator l= skel;
2532    for (CFIterator k= g; k.hasTerms(); k++, l++)
2533    {
2534      if (k.exp() != l.exp())
2535      {
2536        delete[] pEvalPoints;
2537        fail= true;
2538        return 0;
2539      }
2540    }
2541    pEvalPoints[i]= evalPoints;
2542    gcds[i]= g;
2543
2544    tmp= 0;
2545    int j= 0;
2546    for (CFListIterator k= evalPoints; k.hasItem(); k++, j++)
2547      tmp += k.getItem()*power (x, j);
2548    list.append (tmp);
2549  }
2550
2551  if (Monoms.size() == 0)
2552    Monoms= getMonoms (skel);
2553  if (coeffMonoms == NULL)
2554    coeffMonoms= new CFArray [skelSize];
2555  j= 0;
2556  for (CFIterator i= skel; i.hasTerms(); i++, j++)
2557    coeffMonoms[j]= getMonoms (i.coeff());
2558
2559  CFArray* pL= new CFArray [skelSize];
2560  CFArray* pM= new CFArray [skelSize];
2561  for (int i= 0; i < biggestSize; i++)
2562  {
2563    CFIterator l= gcds [i];
2564    evalPoints= pEvalPoints [i];
2565    for (int k= 0; k < skelSize; k++, l++)
2566    {
2567      if (i == 0)
2568        pL[k]= CFArray (biggestSize);
2569      pL[k] [i]= l.coeff();
2570
2571      if (i == 0)
2572        pM[k]= evaluate (coeffMonoms [k], evalPoints);
2573    }
2574  }
2575
2576  CFArray solution;
2577  CanonicalForm result= 0;
2578  int ind= 0;
2579  CFArray bufArray;
2580  CFMatrix Mat;
2581  for (int k= 0; k < skelSize; k++)
2582  {
2583    if (biggestSize != coeffMonoms[k].size())
2584    {
2585      bufArray= CFArray (coeffMonoms[k].size());
2586      for (int i= 0; i < coeffMonoms[k].size(); i++)
2587        bufArray [i]= pL[k] [i];
2588      solution= solveGeneralVandermonde (pM [k], bufArray);
2589    }
2590    else
2591      solution= solveGeneralVandermonde (pM [k], pL[k]);
2592
2593    if (solution.size() == 0)
2594    {
2595      delete[] pEvalPoints;
2596      delete[] pM;
2597      delete[] pL;
2598      delete[] coeffMonoms;
2599      fail= true;
2600      return 0;
2601    }
2602    for (int l= 0; l < solution.size(); l++)
2603      result += solution[l]*Monoms [ind + l];
2604    ind += solution.size();
2605  }
2606
2607  delete[] pEvalPoints;
2608  delete[] pM;
2609  delete[] pL;
2610
2611  if (alpha.level() != 1 && V_buf != alpha)
2612  {
2613    CFList u, v;
2614    result= mapDown (result, prim_elem, im_prim_elem, alpha, u, v);
2615  }
2616
2617  result= N(result);
2618  if (fdivides (result, F) && fdivides (result, G))
2619    return result;
2620  else
2621  {
2622    delete[] coeffMonoms;
2623    fail= true;
2624    return 0;
2625  }
2626}
2627
2628CanonicalForm
2629nonMonicSparseInterpol (const CanonicalForm& F, const CanonicalForm& G,
2630                        const CanonicalForm& skeleton, const Variable& alpha,
2631                        bool& fail, CFArray*& coeffMonoms, CFArray& Monoms
2632                       )
2633{
2634  CanonicalForm A= F;
2635  CanonicalForm B= G;
2636  if (F.isZero() && degree(G) > 0) return G/Lc(G);
2637  else if (G.isZero() && degree (F) > 0) return F/Lc(F);
2638  else if (F.isZero() && G.isZero()) return F.genOne();
2639  if (F.inBaseDomain() || G.inBaseDomain()) return F.genOne();
2640  if (F.isUnivariate() && fdivides(F, G)) return F/Lc(F);
2641  if (G.isUnivariate() && fdivides(G, F)) return G/Lc(G);
2642  if (F == G) return F/Lc(F);
2643
2644  CFMap M,N;
2645  int best_level= myCompress (A, B, M, N, false);
2646
2647  if (best_level == 0)
2648    return B.genOne();
2649
2650  A= M(A);
2651  B= M(B);
2652
2653  Variable x= Variable (1);
2654  ASSERT (degree (A, x) == 0, "expected degree (F, 1) == 0");
2655  ASSERT (degree (B, x) == 0, "expected degree (G, 1) == 0");
2656
2657  //univariate case
2658  if (A.isUnivariate() && B.isUnivariate())
2659    return N (gcd (A, B));
2660
2661  CanonicalForm skel= M(skeleton);
2662
2663  CanonicalForm cA, cB;    // content of A and B
2664  CanonicalForm ppA, ppB;    // primitive part of A and B
2665  CanonicalForm gcdcAcB;
2666  cA = uni_content (A);
2667  cB = uni_content (B);
2668  gcdcAcB= gcd (cA, cB);
2669  ppA= A/cA;
2670  ppB= B/cB;
2671
2672  CanonicalForm lcA, lcB;  // leading coefficients of A and B
2673  CanonicalForm gcdlcAlcB;
2674  lcA= uni_lcoeff (ppA);
2675  lcB= uni_lcoeff (ppB);
2676
2677  if (fdivides (lcA, lcB))
2678  {
2679    if (fdivides (A, B))
2680      return F/Lc(F);
2681  }
2682  if (fdivides (lcB, lcA))
2683  {
2684    if (fdivides (B, A))
2685      return G/Lc(G);
2686  }
2687
2688  gcdlcAlcB= gcd (lcA, lcB);
2689  int skelSize= size (skel, skel.mvar());
2690
2691  int j= 0;
2692  int biggestSize= 0;
2693  int bufSize;
2694  int numberUni= 0;
2695  for (CFIterator i= skel; i.hasTerms(); i++, j++)
2696  {
2697    bufSize= size (i.coeff());
2698    biggestSize= tmax (biggestSize, bufSize);
2699    numberUni += bufSize;
2700  }
2701  numberUni--;
2702  numberUni= (int) ceil ( (double) numberUni / (skelSize - 1));
2703  biggestSize= tmax (biggestSize , numberUni);
2704
2705  numberUni= biggestSize + size (LC(skel)) - 1;
2706  int biggestSize2= tmax (numberUni, biggestSize);
2707
2708  CanonicalForm g, Aeval, Beval;
2709
2710  CFList evalPoints;
2711  CFArray coeffEval;
2712  bool evalFail= false;
2713  CFList list;
2714  bool GF= false;
2715  CanonicalForm LCA= LC (A);
2716  CanonicalForm tmp;
2717  CFArray gcds= CFArray (biggestSize);
2718  CFList * pEvalPoints= new CFList [biggestSize];
2719  Variable V_buf= alpha;
2720  CFList source, dest;
2721  CanonicalForm prim_elem, im_prim_elem;
2722  for (int i= 0; i < biggestSize; i++)
2723  {
2724    if (i == 0)
2725    {
2726      if (getCharacteristic() > 3)
2727        evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf,
2728                                    evalFail, list);
2729      else
2730        evalFail= true;
2731
2732      if (evalFail)
2733      {
2734        if (V_buf.level() != 1)
2735        {
2736          do
2737          {
2738            Variable V_buf2= chooseExtension (V_buf);
2739            source= CFList();
2740            dest= CFList();
2741
2742            bool prim_fail= false;
2743            Variable V_buf3;
2744            prim_elem= primitiveElement (V_buf, V_buf3, prim_fail);
2745
2746            ASSERT (!prim_fail, "failure in integer factorizer");
2747            if (prim_fail)
2748              ; //ERROR
2749            else
2750              im_prim_elem= mapPrimElem (prim_elem, V_buf, V_buf2);
2751
2752            DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf));
2753            DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf2));
2754
2755            for (CFListIterator i= list; i.hasItem(); i++)
2756              i.getItem()= mapUp (i.getItem(), V_buf, V_buf2, prim_elem,
2757                                im_prim_elem, source, dest);
2758            evalFail= false;
2759            evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf,
2760                                          evalFail, list);
2761          } while (evalFail);
2762        }
2763        else
2764        {
2765          CanonicalForm mipo;
2766          int deg= 2;
2767          do {
2768            mipo= randomIrredpoly (deg, x);
2769            V_buf= rootOf (mipo);
2770            evalFail= false;
2771            evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf,
2772                                          evalFail, list);
2773            deg++;
2774          } while (evalFail);
2775        }
2776      }
2777    }
2778    else
2779    {
2780      mult (evalPoints, pEvalPoints[0]);
2781      eval (A, B, Aeval, Beval, evalPoints);
2782    }
2783
2784    g= gcd (Aeval, Beval);
2785    g /= Lc (g);
2786
2787    if (degree (g) != degree (skel) || (skelSize != size (g)))
2788    {
2789      delete[] pEvalPoints;
2790      fail= true;
2791      return 0;
2792    }
2793    CFIterator l= skel;
2794    for (CFIterator k= g; k.hasTerms(); k++, l++)
2795    {
2796      if (k.exp() != l.exp())
2797      {
2798        delete[] pEvalPoints;
2799        fail= true;
2800        return 0;
2801      }
2802    }
2803    pEvalPoints[i]= evalPoints;
2804    gcds[i]= g;
2805
2806    tmp= 0;
2807    int j= 0;
2808    for (CFListIterator k= evalPoints; k.hasItem(); k++, j++)
2809      tmp += k.getItem()*power (x, j);
2810    list.append (tmp);
2811  }
2812
2813  if (Monoms.size() == 0)
2814    Monoms= getMonoms (skel);
2815
2816  if (coeffMonoms == NULL)
2817    coeffMonoms= new CFArray [skelSize];
2818
2819  j= 0;
2820  for (CFIterator i= skel; i.hasTerms(); i++, j++)
2821    coeffMonoms[j]= getMonoms (i.coeff());
2822
2823  int minimalColumnsIndex;
2824  if (skelSize > 1)
2825    minimalColumnsIndex= 1;
2826  else
2827    minimalColumnsIndex= 0;
2828  int minimalColumns=-1;
2829
2830  CFArray* pM= new CFArray [skelSize];
2831  CFMatrix Mat;
2832  // find the Matrix with minimal number of columns
2833  for (int i= 0; i < skelSize; i++)
2834  {
2835    pM[i]= CFArray (coeffMonoms[i].size());
2836    if (i == 1)
2837      minimalColumns= coeffMonoms[i].size();
2838    if (i > 1)
2839    {
2840      if (minimalColumns > coeffMonoms[i].size())
2841      {
2842        minimalColumns= coeffMonoms[i].size();
2843        minimalColumnsIndex= i;
2844      }
2845    }
2846  }
2847  CFMatrix* pMat= new CFMatrix [2];
2848  pMat[0]= CFMatrix (biggestSize, coeffMonoms[0].size());
2849  pMat[1]= CFMatrix (biggestSize, coeffMonoms[minimalColumnsIndex].size());
2850  CFArray* pL= new CFArray [skelSize];
2851  for (int i= 0; i < biggestSize; i++)
2852  {
2853    CFIterator l= gcds [i];
2854    evalPoints= pEvalPoints [i];
2855    for (int k= 0; k < skelSize; k++, l++)
2856    {
2857      if (i == 0)
2858        pL[k]= CFArray (biggestSize);
2859      pL[k] [i]= l.coeff();
2860
2861      if (i == 0 && k != 0 && k != minimalColumnsIndex)
2862        pM[k]= evaluate (coeffMonoms[k], evalPoints);
2863      else if (i == 0 && (k == 0 || k == minimalColumnsIndex))
2864        coeffEval= evaluate (coeffMonoms[k], evalPoints);
2865      if (i > 0 && (k == 0 || k == minimalColumnsIndex))
2866        coeffEval= evaluate (coeffMonoms[k], evalPoints);
2867
2868      if (k == 0)
2869      {
2870        if (pMat[k].rows() >= i + 1)
2871        {
2872          for (int ind= 1; ind < coeffEval.size() + 1; ind++)
2873            pMat[k] (i + 1, ind)= coeffEval[ind - 1];
2874        }
2875      }
2876      if (k == minimalColumnsIndex)
2877      {
2878        if (pMat[1].rows() >= i + 1)
2879        {
2880          for (int ind= 1; ind < coeffEval.size() + 1; ind++)
2881            pMat[1] (i + 1, ind)= coeffEval[ind - 1];
2882        }
2883      }
2884    }
2885  }
2886
2887  CFArray solution;
2888  CanonicalForm result= 0;
2889  int ind= 1;
2890  int matRows, matColumns;
2891  matRows= pMat[1].rows();
2892  matColumns= pMat[0].columns() - 1;
2893  matColumns += pMat[1].columns();
2894
2895  Mat= CFMatrix (matRows, matColumns);
2896  for (int i= 1; i <= matRows; i++)
2897    for (int j= 1; j <= pMat[1].columns(); j++)
2898      Mat (i, j)= pMat[1] (i, j);
2899
2900  for (int j= pMat[1].columns() + 1; j <= matColumns;
2901       j++, ind++)
2902  {
2903    for (int i= 1; i <= matRows; i++)
2904      Mat (i, j)= (-pMat [0] (i, ind + 1))*pL[minimalColumnsIndex] [i - 1];
2905  }
2906
2907  CFArray firstColumn= CFArray (pMat[0].rows());
2908  for (int i= 0; i < pMat[0].rows(); i++)
2909    firstColumn [i]= pMat[0] (i + 1, 1);
2910
2911  CFArray bufArray= pL[minimalColumnsIndex];
2912
2913  for (int i= 0; i < biggestSize; i++)
2914    pL[minimalColumnsIndex] [i] *= firstColumn[i];
2915
2916  CFMatrix bufMat= pMat[1];
2917  pMat[1]= Mat;
2918
2919  if (V_buf.level() != 1)
2920    solution= solveSystemFq (pMat[1],
2921                             pL[minimalColumnsIndex], V_buf);
2922  else
2923    solution= solveSystemFp (pMat[1],
2924                             pL[minimalColumnsIndex]);
2925
2926  if (solution.size() == 0)
2927  { //Javadi's method failed, try deKleine, Monagan, Wittkopf instead
2928    CFMatrix bufMat0= pMat[0];
2929    delete [] pMat;
2930    pMat= new CFMatrix [skelSize];
2931    pL[minimalColumnsIndex]= bufArray;
2932    CFList* bufpEvalPoints= NULL;
2933    CFArray bufGcds;
2934    if (biggestSize != biggestSize2)
2935    {
2936      bufpEvalPoints= pEvalPoints;
2937      pEvalPoints= new CFList [biggestSize2];
2938      bufGcds= gcds;
2939      gcds= CFArray (biggestSize2);
2940      for (int i= 0; i < biggestSize; i++)
2941      {
2942        pEvalPoints[i]= bufpEvalPoints [i];
2943        gcds[i]= bufGcds[i];
2944      }
2945      for (int i= 0; i < biggestSize2 - biggestSize; i++)
2946      {
2947        mult (evalPoints, pEvalPoints[0]);
2948        eval (A, B, Aeval, Beval, evalPoints);
2949        g= gcd (Aeval, Beval);
2950        g /= Lc (g);
2951        if (degree (g) != degree (skel) || (skelSize != size (g)))
2952        {
2953          delete[] pEvalPoints;
2954          delete[] pMat;
2955          delete[] pL;
2956          delete[] coeffMonoms;
2957          delete[] pM;
2958          if (bufpEvalPoints != NULL)
2959            delete [] bufpEvalPoints;
2960          fail= true;
2961          return 0;
2962        }
2963        CFIterator l= skel;
2964        for (CFIterator k= g; k.hasTerms(); k++, l++)
2965        {
2966          if (k.exp() != l.exp())
2967          {
2968            delete[] pEvalPoints;
2969            delete[] pMat;
2970            delete[] pL;
2971            delete[] coeffMonoms;
2972            delete[] pM;
2973            if (bufpEvalPoints != NULL)
2974              delete [] bufpEvalPoints;
2975            fail= true;
2976            return 0;
2977          }
2978        }
2979        pEvalPoints[i + biggestSize]= evalPoints;
2980        gcds[i + biggestSize]= g;
2981      }
2982    }
2983    for (int i= 0; i < biggestSize; i++)
2984    {
2985      CFIterator l= gcds [i];
2986      evalPoints= pEvalPoints [i];
2987      for (int k= 1; k < skelSize; k++, l++)
2988      {
2989        if (i == 0)
2990          pMat[k]= CFMatrix (biggestSize2,coeffMonoms[k].size()+biggestSize2-1);
2991        if (k == minimalColumnsIndex)
2992          continue;
2993        coeffEval= evaluate (coeffMonoms[k], evalPoints);
2994        if (pMat[k].rows() >= i + 1)
2995        {
2996          for (int ind= 1; ind < coeffEval.size() + 1; ind++)
2997            pMat[k] (i + 1, ind)= coeffEval[ind - 1];
2998        }
2999      }
3000    }
3001    Mat= bufMat0;
3002    pMat[0]= CFMatrix (biggestSize2, coeffMonoms[0].size() + biggestSize2 - 1);
3003    for (int j= 1; j <= Mat.rows(); j++)
3004      for (int k= 1; k <= Mat.columns(); k++)
3005        pMat [0] (j,k)= Mat (j,k);
3006    Mat= bufMat;
3007    for (int j= 1; j <= Mat.rows(); j++)
3008      for (int k= 1; k <= Mat.columns(); k++)
3009        pMat [minimalColumnsIndex] (j,k)= Mat (j,k);
3010    // write old matrix entries into new matrices
3011    for (int i= 0; i < skelSize; i++)
3012    {
3013      bufArray= pL[i];
3014      pL[i]= CFArray (biggestSize2);
3015      for (int j= 0; j < bufArray.size(); j++)
3016        pL[i] [j]= bufArray [j];
3017    }
3018    //write old vector entries into new and add new entries to old matrices
3019    for (int i= 0; i < biggestSize2 - biggestSize; i++)
3020    {
3021      CFIterator l= gcds [i + biggestSize];
3022      evalPoints= pEvalPoints [i + biggestSize];
3023      for (int k= 0; k < skelSize; k++, l++)
3024      {
3025        pL[k] [i + biggestSize]= l.coeff();
3026        coeffEval= evaluate (coeffMonoms[k], evalPoints);
3027        if (pMat[k].rows() >= i + biggestSize + 1)
3028        {
3029          for (int ind= 1; ind < coeffEval.size() + 1; ind++)
3030            pMat[k] (i + biggestSize + 1, ind)= coeffEval[ind - 1];
3031        }
3032      }
3033    }
3034    // begin new
3035    for (int i= 0; i < skelSize; i++)
3036    {
3037      if (pL[i].size() > 1)
3038      {
3039        for (int j= 2; j <= pMat[i].rows(); j++)
3040          pMat[i] (j, coeffMonoms[i].size() + j - 1)=
3041              -pL[i] [j - 1];
3042      }
3043    }
3044
3045    matColumns= biggestSize2 - 1;
3046    matRows= 0;
3047    for (int i= 0; i < skelSize; i++)
3048    {
3049      if (V_buf.level() == 1)
3050        (void) gaussianElimFp (pMat[i], pL[i]);
3051      else
3052        (void) gaussianElimFq (pMat[i], pL[i], V_buf);
3053
3054      if (pMat[i] (coeffMonoms[i].size(), coeffMonoms[i].size()) == 0)
3055      {
3056        delete[] pEvalPoints;
3057        delete[] pMat;
3058        delete[] pL;
3059        delete[] coeffMonoms;
3060        delete[] pM;
3061        if (bufpEvalPoints != NULL)
3062          delete [] bufpEvalPoints;
3063        fail= true;
3064        return 0;
3065      }
3066      matRows += pMat[i].rows() - coeffMonoms[i].size();
3067    }
3068
3069    CFMatrix bufMat;
3070    Mat= CFMatrix (matRows, matColumns);
3071    ind= 0;
3072    bufArray= CFArray (matRows);
3073    CFArray bufArray2;
3074    for (int i= 0; i < skelSize; i++)
3075    {
3076      bufMat= pMat[i] (coeffMonoms[i].size() + 1, pMat[i].rows(),
3077                       coeffMonoms[i].size() + 1, pMat[i].columns());
3078
3079      for (int j= 1; j <= bufMat.rows(); j++)
3080        for (int k= 1; k <= bufMat.columns(); k++)
3081          Mat (j + ind, k)= bufMat(j, k);
3082      bufArray2= coeffMonoms[i].size();
3083      for (int j= 1; j <= pMat[i].rows(); j++)
3084      {
3085        if (j > coeffMonoms[i].size())
3086          bufArray [j-coeffMonoms[i].size() + ind - 1]= pL[i] [j - 1];
3087        else
3088          bufArray2 [j - 1]= pL[i] [j - 1];
3089      }
3090      pL[i]= bufArray2;
3091      ind += bufMat.rows();
3092      pMat[i]= pMat[i] (1, coeffMonoms[i].size(), 1, pMat[i].columns());
3093    }
3094
3095    if (V_buf.level() != 1)
3096      solution= solveSystemFq (Mat, bufArray, V_buf);
3097    else
3098      solution= solveSystemFp (Mat, bufArray);
3099
3100    if (solution.size() == 0)
3101    {
3102      delete[] pEvalPoints;
3103      delete[] pMat;
3104      delete[] pL;
3105      delete[] coeffMonoms;
3106      delete[] pM;
3107      if (bufpEvalPoints != NULL)
3108        delete [] bufpEvalPoints;
3109      fail= true;
3110      return 0;
3111    }
3112
3113    ind= 0;
3114    result= 0;
3115    CFArray bufSolution;
3116    for (int i= 0; i < skelSize; i++)
3117    {
3118      bufSolution= readOffSolution (pMat[i], pL[i], solution);
3119      for (int i= 0; i < bufSolution.size(); i++, ind++)
3120        result += Monoms [ind]*bufSolution[i];
3121    }
3122    if (alpha.level() != 1 && V_buf != alpha)
3123    {
3124      CFList u, v;
3125      result= mapDown (result, prim_elem, im_prim_elem, alpha, u, v);
3126    }
3127    result= N(result);
3128    if (fdivides (result, F) && fdivides (result, G))
3129    {
3130      delete[] pEvalPoints;
3131      delete[] pMat;
3132      delete[] pL;
3133      delete[] pM;
3134      if (bufpEvalPoints != NULL)
3135        delete [] bufpEvalPoints;
3136      return result;
3137    }
3138    else
3139    {
3140      delete[] pEvalPoints;
3141      delete[] pMat;
3142      delete[] pL;
3143      delete[] coeffMonoms;
3144      delete[] pM;
3145      if (bufpEvalPoints != NULL)
3146        delete [] bufpEvalPoints;
3147      fail= true;
3148      return 0;
3149    }
3150  } // end of deKleine, Monagan & Wittkopf
3151
3152  result += Monoms[0];
3153  int ind2= 0, ind3= 2;
3154  ind= 0;
3155  for (int l= 0; l < minimalColumnsIndex; l++)
3156    ind += coeffMonoms[l].size();
3157  for (int l= solution.size() - pMat[0].columns() + 1; l < solution.size();
3158       l++, ind2++, ind3++)
3159  {
3160    result += solution[l]*Monoms [1 + ind2];
3161    for (int i= 0; i < pMat[0].rows(); i++)
3162      firstColumn[i] += solution[l]*pMat[0] (i + 1, ind3);
3163  }
3164  for (int l= 0; l < solution.size() + 1 - pMat[0].columns(); l++)
3165    result += solution[l]*Monoms [ind + l];
3166  ind= coeffMonoms[0].size();
3167  for (int k= 1; k < skelSize; k++)
3168  {
3169    if (k == minimalColumnsIndex)
3170    {
3171      ind += coeffMonoms[k].size();
3172      continue;
3173    }
3174    if (k != minimalColumnsIndex)
3175    {
3176      for (int i= 0; i < biggestSize; i++)
3177        pL[k] [i] *= firstColumn [i];
3178    }
3179
3180    if (biggestSize != coeffMonoms[k].size() && k != minimalColumnsIndex)
3181    {
3182      bufArray= CFArray (coeffMonoms[k].size());
3183      for (int i= 0; i < bufArray.size(); i++)
3184        bufArray [i]= pL[k] [i];
3185      solution= solveGeneralVandermonde (pM[k], bufArray);
3186    }
3187    else
3188      solution= solveGeneralVandermonde (pM[k], pL[k]);
3189
3190    if (solution.size() == 0)
3191    {
3192      delete[] pEvalPoints;
3193      delete[] pMat;
3194      delete[] pL;
3195      delete[] coeffMonoms;
3196      delete[] pM;
3197      fail= true;
3198      return 0;
3199    }
3200    if (k != minimalColumnsIndex)
3201    {
3202      for (int l= 0; l < solution.size(); l++)
3203        result += solution[l]*Monoms [ind + l];
3204      ind += solution.size();
3205    }
3206  }
3207
3208  delete[] pEvalPoints;
3209  delete[] pMat;
3210  delete[] pL;
3211  delete[] pM;
3212
3213  if (alpha.level() != 1 && V_buf != alpha)
3214  {
3215    CFList u, v;
3216    result= mapDown (result, prim_elem, im_prim_elem, alpha, u, v);
3217  }
3218  result= N(result);
3219
3220  if (fdivides (result, F) && fdivides (result, G))
3221    return result;
3222  else
3223  {
3224    delete[] coeffMonoms;
3225    fail= true;
3226    return 0;
3227  }
3228}
3229
3230CanonicalForm sparseGCDFq (const CanonicalForm& F, const CanonicalForm& G,
3231                           const Variable & alpha, CFList& l, bool& topLevel)
3232{
3233  CanonicalForm A= F;
3234  CanonicalForm B= G;
3235  if (F.isZero() && degree(G) > 0) return G/Lc(G);
3236  else if (G.isZero() && degree (F) > 0) return F/Lc(F);
3237  else if (F.isZero() && G.isZero()) return F.genOne();
3238  if (F.inBaseDomain() || G.inBaseDomain()) return F.genOne();
3239  if (F.isUnivariate() && fdivides(F, G)) return F/Lc(F);
3240  if (G.isUnivariate() && fdivides(G, F)) return G/Lc(G);
3241  if (F == G) return F/Lc(F);
3242
3243  CFMap M,N;
3244  int best_level= myCompress (A, B, M, N, topLevel);
3245
3246  if (best_level == 0) return B.genOne();
3247
3248  A= M(A);
3249  B= M(B);
3250
3251  Variable x= Variable (1);
3252
3253  //univariate case
3254  if (A.isUnivariate() && B.isUnivariate())
3255    return N (gcd (A, B));
3256
3257  CanonicalForm cA, cB;    // content of A and B
3258  CanonicalForm ppA, ppB;    // primitive part of A and B
3259  CanonicalForm gcdcAcB;
3260
3261  cA = uni_content (A);
3262  cB = uni_content (B);
3263  gcdcAcB= gcd (cA, cB);
3264  ppA= A/cA;
3265  ppB= B/cB;
3266
3267  CanonicalForm lcA, lcB;  // leading coefficients of A and B
3268  CanonicalForm gcdlcAlcB;
3269  lcA= uni_lcoeff (ppA);
3270  lcB= uni_lcoeff (ppB);
3271
3272  if (fdivides (lcA, lcB))
3273  {
3274    if (fdivides (A, B))
3275      return F/Lc(F);
3276  }
3277  if (fdivides (lcB, lcA))
3278  {
3279    if (fdivides (B, A))
3280      return G/Lc(G);
3281  }
3282
3283  gcdlcAlcB= gcd (lcA, lcB);
3284
3285  int d= totaldegree (ppA, Variable (2), Variable (ppA.level()));
3286  int d0;
3287
3288  if (d == 0)
3289    return N(gcdcAcB);
3290  d0= totaldegree (ppB, Variable (2), Variable (ppB.level()));
3291
3292  if (d0 < d)
3293    d= d0;
3294
3295  if (d == 0)
3296    return N(gcdcAcB);
3297
3298  CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH, skeleton;
3299  CanonicalForm newtonPoly= 1;
3300  m= gcdlcAlcB;
3301  G_m= 0;
3302  H= 0;
3303  bool fail= false;
3304  topLevel= false;
3305  bool inextension= false;
3306  Variable V_buf= alpha;
3307  CanonicalForm prim_elem, im_prim_elem;
3308  CFList source, dest;
3309  do // first do
3310  {
3311    random_element= randomElement (m, V_buf, l, fail);
3312    if (random_element == 0 && !fail)
3313    {
3314      if (!find (l, random_element))
3315        l.append (random_element);
3316      continue;
3317    }
3318    if (fail)
3319    {
3320      source= CFList();
3321      dest= CFList();
3322
3323      Variable V_buf3= V_buf;
3324      V_buf= chooseExtension (V_buf);
3325      bool prim_fail= false;
3326      Variable V_buf2;
3327      prim_elem= primitiveElement (alpha, V_buf2, prim_fail);
3328
3329      if (V_buf3 != alpha)
3330      {
3331        m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
3332        G_m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
3333        newtonPoly= mapDown (newtonPoly, prim_elem, im_prim_elem, alpha,
3334                             source, dest);
3335        ppA= mapDown (ppA, prim_elem, im_prim_elem, alpha, source, dest);
3336        ppB= mapDown (ppB, prim_elem, im_prim_elem, alpha, source, dest);
3337        gcdlcAlcB= mapDown (gcdlcAlcB, prim_elem, im_prim_elem, alpha, source,
3338                            dest);
3339        for (CFListIterator i= l; i.hasItem(); i++)
3340          i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, alpha,
3341                                source, dest);
3342      }
3343
3344      ASSERT (!prim_fail, "failure in integer factorizer");
3345      if (prim_fail)
3346        ; //ERROR
3347      else
3348        im_prim_elem= mapPrimElem (prim_elem, alpha, V_buf);
3349
3350      DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (alpha));
3351      DEBOUTLN (cerr, "getMipo (V_buf2)= " << getMipo (V_buf2));
3352      inextension= true;
3353      for (CFListIterator i= l; i.hasItem(); i++)
3354        i.getItem()= mapUp (i.getItem(), alpha, V_buf, prim_elem,
3355                             im_prim_elem, source, dest);
3356      m= mapUp (m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3357      G_m= mapUp (G_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3358      newtonPoly= mapUp (newtonPoly, alpha, V_buf, prim_elem, im_prim_elem,
3359                          source, dest);
3360      ppA= mapUp (ppA, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3361      ppB= mapUp (ppB, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3362      gcdlcAlcB= mapUp (gcdlcAlcB, alpha, V_buf, prim_elem, im_prim_elem,
3363                         source, dest);
3364
3365      fail= false;
3366      random_element= randomElement (m, V_buf, l, fail );
3367      DEBOUTLN (cerr, "fail= " << fail);
3368      CFList list;
3369      TIMING_START (gcd_recursion);
3370      G_random_element=
3371      sparseGCDFq (ppA (random_element, x), ppB (random_element, x), V_buf,
3372                        list, topLevel);
3373      TIMING_END_AND_PRINT (gcd_recursion,
3374                            "time for recursive call: ");
3375      DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3376    }
3377    else
3378    {
3379      CFList list;
3380      TIMING_START (gcd_recursion);
3381      G_random_element=
3382      sparseGCDFq (ppA(random_element, x), ppB(random_element, x), V_buf,
3383                        list, topLevel);
3384      TIMING_END_AND_PRINT (gcd_recursion,
3385                            "time for recursive call: ");
3386      DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3387    }
3388
3389    if (!G_random_element.inCoeffDomain())
3390      d0= totaldegree (G_random_element, Variable(2),
3391                       Variable (G_random_element.level()));
3392    else
3393      d0= 0;
3394
3395    if (d0 == 0)
3396      return N(gcdcAcB);
3397    if (d0 >  d)
3398    {
3399      if (!find (l, random_element))
3400        l.append (random_element);
3401      continue;
3402    }
3403
3404    G_random_element=
3405    (gcdlcAlcB(random_element, x)/uni_lcoeff (G_random_element))
3406    * G_random_element;
3407
3408    skeleton= G_random_element;
3409    if (!G_random_element.inCoeffDomain())
3410      d0= totaldegree (G_random_element, Variable(2),
3411                       Variable (G_random_element.level()));
3412    else
3413      d0= 0;
3414
3415    if (d0 <  d)
3416    {
3417      m= gcdlcAlcB;
3418      newtonPoly= 1;
3419      G_m= 0;
3420      d= d0;
3421    }
3422
3423    H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
3424    if (uni_lcoeff (H) == gcdlcAlcB)
3425    {
3426      cH= uni_content (H);
3427      ppH= H/cH;
3428      if (inextension)
3429      {
3430        CFList u, v;
3431        //maybe it's better to test if ppH is an element of F(\alpha) before
3432        //mapping down
3433        if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
3434        {
3435          DEBOUTLN (cerr, "ppH before mapDown= " << ppH);
3436          ppH= mapDown (ppH, prim_elem, im_prim_elem, alpha, u, v);
3437          ppH /= Lc(ppH);
3438          DEBOUTLN (cerr, "ppH after mapDown= " << ppH);
3439          return N(gcdcAcB*ppH);
3440        }
3441      }
3442      else if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
3443        return N(gcdcAcB*ppH);
3444    }
3445    G_m= H;
3446    newtonPoly= newtonPoly*(x - random_element);
3447    m= m*(x - random_element);
3448    if (!find (l, random_element))
3449      l.append (random_element);
3450
3451    if (getCharacteristic () > 3 && size (skeleton) < 100)
3452    {
3453      CFArray Monoms;
3454      CFArray *coeffMonoms= NULL;
3455      do //second do
3456      {
3457        random_element= randomElement (m, V_buf, l, fail);
3458        if (random_element == 0 && !fail)
3459        {
3460          if (!find (l, random_element))
3461            l.append (random_element);
3462          continue;
3463        }
3464        if (fail)
3465        {
3466          source= CFList();
3467          dest= CFList();
3468
3469          Variable V_buf3= V_buf;
3470          V_buf= chooseExtension (V_buf);
3471          bool prim_fail= false;
3472          Variable V_buf2;
3473          prim_elem= primitiveElement (alpha, V_buf2, prim_fail);
3474
3475          if (V_buf3 != alpha)
3476          {
3477            m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
3478            G_m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
3479            newtonPoly= mapDown (newtonPoly, prim_elem, im_prim_elem, alpha,
3480                                 source, dest);
3481            ppA= mapDown (ppA, prim_elem, im_prim_elem, alpha, source, dest);
3482            ppB= mapDown (ppB, prim_elem, im_prim_elem, alpha, source, dest);
3483            gcdlcAlcB= mapDown (gcdlcAlcB, prim_elem, im_prim_elem, alpha,
3484                                source, dest);
3485            for (CFListIterator i= l; i.hasItem(); i++)
3486              i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, alpha,
3487                                    source, dest);
3488          }
3489
3490          ASSERT (!prim_fail, "failure in integer factorizer");
3491          if (prim_fail)
3492            ; //ERROR
3493          else
3494            im_prim_elem= mapPrimElem (prim_elem, alpha, V_buf);
3495
3496          DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (alpha));
3497          DEBOUTLN (cerr, "getMipo (V_buf2)= " << getMipo (V_buf2));
3498          inextension= true;
3499          for (CFListIterator i= l; i.hasItem(); i++)
3500            i.getItem()= mapUp (i.getItem(), alpha, V_buf, prim_elem,
3501                                im_prim_elem, source, dest);
3502          m= mapUp (m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3503          G_m= mapUp (G_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3504          newtonPoly= mapUp (newtonPoly, alpha, V_buf, prim_elem, im_prim_elem,
3505                              source, dest);
3506          ppA= mapUp (ppA, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3507          ppB= mapUp (ppB, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3508
3509          gcdlcAlcB= mapUp (gcdlcAlcB, alpha, V_buf, prim_elem, im_prim_elem,
3510                            source, dest);
3511
3512          fail= false;
3513          random_element= randomElement (m, V_buf, l, fail);
3514          DEBOUTLN (cerr, "fail= " << fail);
3515          CFList list;
3516          TIMING_START (gcd_recursion);
3517
3518          //sparseInterpolation
3519          bool sparseFail= false;
3520          if (LC (skeleton).inCoeffDomain())
3521            G_random_element=
3522            monicSparseInterpol (ppA (random_element, x), ppB(random_element,x),
3523                                skeleton,V_buf, sparseFail, coeffMonoms,Monoms);
3524          else
3525            G_random_element=
3526            nonMonicSparseInterpol (ppA(random_element,x),ppB(random_element,x),
3527                                    skeleton, V_buf, sparseFail, coeffMonoms,
3528                                    Monoms);
3529          if (sparseFail)
3530            break;
3531
3532          TIMING_END_AND_PRINT (gcd_recursion,
3533                                "time for recursive call: ");
3534          DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3535        }
3536        else
3537        {
3538          CFList list;
3539          TIMING_START (gcd_recursion);
3540          bool sparseFail= false;
3541          if (LC (skeleton).inCoeffDomain())
3542            G_random_element=
3543            monicSparseInterpol (ppA (random_element, x),ppB(random_element, x),
3544                                skeleton,V_buf, sparseFail,coeffMonoms, Monoms);
3545          else
3546            G_random_element=
3547            nonMonicSparseInterpol (ppA(random_element,x),ppB(random_element,x),
3548                                    skeleton, V_buf, sparseFail, coeffMonoms,
3549                                    Monoms);
3550          if (sparseFail)
3551            break;
3552
3553          TIMING_END_AND_PRINT (gcd_recursion,
3554                                "time for recursive call: ");
3555          DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3556        }
3557
3558        if (!G_random_element.inCoeffDomain())
3559          d0= totaldegree (G_random_element, Variable(2),
3560                           Variable (G_random_element.level()));
3561        else
3562          d0= 0;
3563
3564        if (d0 == 0)
3565          return N(gcdcAcB);
3566        if (d0 >  d)
3567        {
3568          if (!find (l, random_element))
3569            l.append (random_element);
3570          continue;
3571        }
3572
3573        G_random_element=
3574        (gcdlcAlcB(random_element, x)/uni_lcoeff (G_random_element))
3575        * G_random_element;
3576
3577        if (!G_random_element.inCoeffDomain())
3578          d0= totaldegree (G_random_element, Variable(2),
3579                          Variable (G_random_element.level()));
3580        else
3581          d0= 0;
3582
3583        if (d0 <  d)
3584        {
3585          m= gcdlcAlcB;
3586          newtonPoly= 1;
3587          G_m= 0;
3588          d= d0;
3589        }
3590
3591        TIMING_START (newton_interpolation);
3592        H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
3593        TIMING_END_AND_PRINT (newton_interpolation,
3594                              "time for newton interpolation: ");
3595
3596        //termination test
3597        if (uni_lcoeff (H) == gcdlcAlcB)
3598        {
3599          cH= uni_content (H);
3600          ppH= H/cH;
3601          if (inextension)
3602          {
3603            CFList u, v;
3604            //maybe it's better to test if ppH is an element of F(\alpha) before
3605            //mapping down
3606            if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
3607            {
3608              DEBOUTLN (cerr, "ppH before mapDown= " << ppH);
3609              ppH= mapDown (ppH, prim_elem, im_prim_elem, alpha, u, v);
3610              ppH /= Lc(ppH);
3611              DEBOUTLN (cerr, "ppH after mapDown= " << ppH);
3612              return N(gcdcAcB*ppH);
3613            }
3614          }
3615          else if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
3616          {
3617            return N(gcdcAcB*ppH);
3618          }
3619        }
3620
3621        G_m= H;
3622        newtonPoly= newtonPoly*(x - random_element);
3623        m= m*(x - random_element);
3624        if (!find (l, random_element))
3625          l.append (random_element);
3626
3627      } while (1);
3628    }
3629  } while (1);
3630}
3631
3632CanonicalForm sparseGCDFp (const CanonicalForm& F, const CanonicalForm& G,
3633                           bool& topLevel, CFList& l)
3634{
3635  CanonicalForm A= F;
3636  CanonicalForm B= G;
3637  if (F.isZero() && degree(G) > 0) return G/Lc(G);
3638  else if (G.isZero() && degree (F) > 0) return F/Lc(F);
3639  else if (F.isZero() && G.isZero()) return F.genOne();
3640  if (F.inBaseDomain() || G.inBaseDomain()) return F.genOne();
3641  if (F.isUnivariate() && fdivides(F, G)) return F/Lc(F);
3642  if (G.isUnivariate() && fdivides(G, F)) return G/Lc(G);
3643  if (F == G) return F/Lc(F);
3644
3645  CFMap M,N;
3646  int best_level= myCompress (A, B, M, N, topLevel);
3647
3648  if (best_level == 0) return B.genOne();
3649
3650  A= M(A);
3651  B= M(B);
3652
3653  Variable x= Variable (1);
3654
3655  //univariate case
3656  if (A.isUnivariate() && B.isUnivariate())
3657    return N (gcd (A, B));
3658
3659  CanonicalForm cA, cB;    // content of A and B
3660  CanonicalForm ppA, ppB;    // primitive part of A and B
3661  CanonicalForm gcdcAcB;
3662
3663  cA = uni_content (A);
3664  cB = uni_content (B);
3665  gcdcAcB= gcd (cA, cB);
3666  ppA= A/cA;
3667  ppB= B/cB;
3668
3669  CanonicalForm lcA, lcB;  // leading coefficients of A and B
3670  CanonicalForm gcdlcAlcB;
3671  lcA= uni_lcoeff (ppA);
3672  lcB= uni_lcoeff (ppB);
3673
3674  if (fdivides (lcA, lcB))
3675  {
3676    if (fdivides (A, B))
3677      return F/Lc(F);
3678  }
3679  if (fdivides (lcB, lcA))
3680  {
3681    if (fdivides (B, A))
3682      return G/Lc(G);
3683  }
3684
3685  gcdlcAlcB= gcd (lcA, lcB);
3686
3687  int d= totaldegree (ppA, Variable (2), Variable (ppA.level()));
3688  int d0;
3689
3690  if (d == 0)
3691    return N(gcdcAcB);
3692
3693  d0= totaldegree (ppB, Variable (2), Variable (ppB.level()));
3694
3695  if (d0 < d)
3696    d= d0;
3697
3698  if (d == 0)
3699    return N(gcdcAcB);
3700
3701  CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH, skeleton;
3702  CanonicalForm newtonPoly= 1;
3703  m= gcdlcAlcB;
3704  G_m= 0;
3705  H= 0;
3706  bool fail= false;
3707  topLevel= false;
3708  bool inextension= false;
3709  Variable V_buf, alpha;
3710  CanonicalForm prim_elem, im_prim_elem;
3711  CFList source, dest;
3712  do //first do
3713  {
3714    if (inextension)
3715      random_element= randomElement (m, V_buf, l, fail);
3716    else
3717      random_element= FpRandomElement (m, l, fail);
3718    if (random_element == 0 && !fail)
3719    {
3720      if (!find (l, random_element))
3721        l.append (random_element);
3722      continue;
3723    }
3724
3725    if (!fail && !inextension)
3726    {
3727      CFList list;
3728      TIMING_START (gcd_recursion);
3729      G_random_element=
3730      sparseGCDFp (ppA (random_element,x), ppB (random_element,x), topLevel,
3731                   list);
3732      TIMING_END_AND_PRINT (gcd_recursion,
3733                            "time for recursive call: ");
3734      DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3735    }
3736    else if (!fail && inextension)
3737    {
3738      CFList list;
3739      TIMING_START (gcd_recursion);
3740      G_random_element=
3741      sparseGCDFq (ppA (random_element, x), ppB (random_element, x), alpha,
3742                        list, topLevel);
3743      TIMING_END_AND_PRINT (gcd_recursion,
3744                            "time for recursive call: ");
3745      DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3746    }
3747    else if (fail && !inextension)
3748    {
3749      source= CFList();
3750      dest= CFList();
3751      CFList list;
3752      CanonicalForm mipo;
3753      int deg= 2;
3754      do
3755      {
3756        mipo= randomIrredpoly (deg, x);
3757        alpha= rootOf (mipo);
3758        inextension= true;
3759        fail= false;
3760        random_element= randomElement (m, alpha, l, fail);
3761        deg++;
3762      } while (fail);
3763      V_buf= alpha;
3764      list= CFList();
3765      TIMING_START (gcd_recursion);
3766      G_random_element=
3767      sparseGCDFq (ppA (random_element, x), ppB (random_element, x), alpha,
3768                        list, topLevel);
3769      TIMING_END_AND_PRINT (gcd_recursion,
3770                            "time for recursive call: ");
3771      DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3772    }
3773    else if (fail && inextension)
3774    {
3775      source= CFList();
3776      dest= CFList();
3777
3778      Variable V_buf3= V_buf;
3779      V_buf= chooseExtension (V_buf);
3780      bool prim_fail= false;
3781      Variable V_buf2;
3782      CanonicalForm prim_elem, im_prim_elem;
3783      prim_elem= primitiveElement (alpha, V_buf2, prim_fail);
3784
3785      if (V_buf3 != alpha)
3786      {
3787        m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
3788        G_m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
3789        newtonPoly= mapDown (newtonPoly, prim_elem, im_prim_elem, alpha, source,
3790                             dest);
3791        ppA= mapDown (ppA, prim_elem, im_prim_elem, alpha, source, dest);
3792        ppB= mapDown (ppB, prim_elem, im_prim_elem, alpha, source, dest);
3793        gcdlcAlcB= mapDown (gcdlcAlcB, prim_elem, im_prim_elem, alpha, source,
3794                            dest);
3795        for (CFListIterator i= l; i.hasItem(); i++)
3796          i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, alpha,
3797                                source, dest);
3798      }
3799
3800      ASSERT (!prim_fail, "failure in integer factorizer");
3801      if (prim_fail)
3802        ; //ERROR
3803      else
3804        im_prim_elem= mapPrimElem (prim_elem, alpha, V_buf);
3805
3806      DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (alpha));
3807      DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf2));
3808
3809      for (CFListIterator i= l; i.hasItem(); i++)
3810        i.getItem()= mapUp (i.getItem(), alpha, V_buf, prim_elem,
3811                             im_prim_elem, source, dest);
3812      m= mapUp (m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3813      G_m= mapUp (G_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3814      newtonPoly= mapUp (newtonPoly, alpha, V_buf, prim_elem, im_prim_elem,
3815                          source, dest);
3816      ppA= mapUp (ppA, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3817      ppB= mapUp (ppB, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3818      gcdlcAlcB= mapUp (gcdlcAlcB, alpha, V_buf, prim_elem, im_prim_elem,
3819                         source, dest);
3820      fail= false;
3821      random_element= randomElement (m, V_buf, l, fail );
3822      DEBOUTLN (cerr, "fail= " << fail);
3823      CFList list;
3824      TIMING_START (gcd_recursion);
3825      G_random_element=
3826      sparseGCDFq (ppA (random_element, x), ppB (random_element, x), V_buf,
3827                        list, topLevel);
3828      TIMING_END_AND_PRINT (gcd_recursion,
3829                            "time for recursive call: ");
3830      DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3831    }
3832
3833    if (!G_random_element.inCoeffDomain())
3834      d0= totaldegree (G_random_element, Variable(2),
3835                       Variable (G_random_element.level()));
3836    else
3837      d0= 0;
3838
3839    if (d0 == 0)
3840      return N(gcdcAcB);
3841    if (d0 >  d)
3842    {
3843      if (!find (l, random_element))
3844        l.append (random_element);
3845      continue;
3846    }
3847
3848    G_random_element=
3849    (gcdlcAlcB(random_element, x)/uni_lcoeff (G_random_element))
3850    * G_random_element;
3851
3852    skeleton= G_random_element;
3853
3854    if (!G_random_element.inCoeffDomain())
3855      d0= totaldegree (G_random_element, Variable(2),
3856                       Variable (G_random_element.level()));
3857    else
3858      d0= 0;
3859
3860    if (d0 <  d)
3861    {
3862      m= gcdlcAlcB;
3863      newtonPoly= 1;
3864      G_m= 0;
3865      d= d0;
3866    }
3867
3868    H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
3869
3870    if (uni_lcoeff (H) == gcdlcAlcB)
3871    {
3872      cH= uni_content (H);
3873      ppH= H/cH;
3874      ppH /= Lc (ppH);
3875      DEBOUTLN (cerr, "ppH= " << ppH);
3876
3877      if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
3878        return N(gcdcAcB*ppH);
3879    }
3880    G_m= H;
3881    newtonPoly= newtonPoly*(x - random_element);
3882    m= m*(x - random_element);
3883    if (!find (l, random_element))
3884      l.append (random_element);
3885
3886    if ((getCharacteristic() > 3 && size (skeleton) < 200))
3887    {
3888      CFArray Monoms;
3889      CFArray* coeffMonoms= NULL;
3890
3891      do //second do
3892      {
3893        if (inextension)
3894          random_element= randomElement (m, alpha, l, fail);
3895        else
3896          random_element= FpRandomElement (m, l, fail);
3897        if (random_element == 0 && !fail)
3898        {
3899          if (!find (l, random_element))
3900            l.append (random_element);
3901          continue;
3902        }
3903
3904        bool sparseFail= false;
3905        if (!fail && !inextension)
3906        {
3907          CFList list;
3908          TIMING_START (gcd_recursion);
3909          if (LC (skeleton).inCoeffDomain())
3910            G_random_element=
3911            monicSparseInterpol(ppA(random_element, x), ppB (random_element, x),
3912                                skeleton, Variable(1), sparseFail, coeffMonoms,
3913                                Monoms);
3914          else
3915            G_random_element=
3916            nonMonicSparseInterpol(ppA(random_element,x), ppB(random_element,x),
3917                                    skeleton, Variable (1), sparseFail,
3918                                    coeffMonoms, Monoms);
3919          TIMING_END_AND_PRINT (gcd_recursion,
3920                                "time for recursive call: ");
3921          DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3922        }
3923        else if (!fail && inextension)
3924        {
3925          CFList list;
3926          TIMING_START (gcd_recursion);
3927          if (LC (skeleton).inCoeffDomain())
3928            G_random_element=
3929            monicSparseInterpol(ppA (random_element,x), ppB (random_element, x),
3930                                skeleton, alpha, sparseFail, coeffMonoms,
3931                                Monoms);
3932          else
3933            G_random_element=
3934            nonMonicSparseInterpol(ppA(random_element,x), ppB(random_element,x),
3935                                   skeleton, alpha, sparseFail, coeffMonoms,
3936                                   Monoms);
3937          TIMING_END_AND_PRINT (gcd_recursion,
3938                                "time for recursive call: ");
3939          DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3940        }
3941        else if (fail && !inextension)
3942        {
3943          source= CFList();
3944          dest= CFList();
3945          CFList list;
3946          CanonicalForm mipo;
3947          int deg= 2;
3948          do
3949          {
3950            mipo= randomIrredpoly (deg, x);
3951            alpha= rootOf (mipo);
3952            inextension= true;
3953            fail= false;
3954            random_element= randomElement (m, alpha, l, fail);
3955            deg++;
3956          } while (fail);
3957          V_buf= alpha;
3958          list= CFList();
3959          TIMING_START (gcd_recursion);
3960          if (LC (skeleton).inCoeffDomain())
3961            G_random_element=
3962            monicSparseInterpol (ppA (random_element,x), ppB (random_element,x),
3963                                skeleton, alpha, sparseFail, coeffMonoms,
3964                                Monoms);
3965          else
3966            G_random_element=
3967            nonMonicSparseInterpol(ppA(random_element,x), ppB(random_element,x),
3968                                   skeleton, alpha, sparseFail, coeffMonoms,
3969                                   Monoms);
3970          TIMING_END_AND_PRINT (gcd_recursion,
3971                                "time for recursive call: ");
3972          DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3973        }
3974        else if (fail && inextension)
3975        {
3976          source= CFList();
3977          dest= CFList();
3978
3979          Variable V_buf3= V_buf;
3980          V_buf= chooseExtension (V_buf);
3981          bool prim_fail= false;
3982          Variable V_buf2;
3983          CanonicalForm prim_elem, im_prim_elem;
3984          prim_elem= primitiveElement (alpha, V_buf2, prim_fail);
3985
3986          if (V_buf3 != alpha)
3987          {
3988            m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
3989            G_m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
3990            newtonPoly= mapDown (newtonPoly, prim_elem, im_prim_elem, alpha,
3991                                 source, dest);
3992            ppA= mapDown (ppA, prim_elem, im_prim_elem, alpha, source, dest);
3993            ppB= mapDown (ppB, prim_elem, im_prim_elem, alpha, source, dest);
3994            gcdlcAlcB= mapDown (gcdlcAlcB, prim_elem, im_prim_elem, alpha,
3995                                source, dest);
3996            for (CFListIterator i= l; i.hasItem(); i++)
3997              i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, alpha,
3998                                    source, dest);
3999          }
4000
4001          ASSERT (!prim_fail, "failure in integer factorizer");
4002          if (prim_fail)
4003            ; //ERROR
4004          else
4005            im_prim_elem= mapPrimElem (prim_elem, alpha, V_buf);
4006
4007          DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (alpha));
4008          DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf2));
4009
4010          for (CFListIterator i= l; i.hasItem(); i++)
4011            i.getItem()= mapUp (i.getItem(), alpha, V_buf, prim_elem,
4012                                im_prim_elem, source, dest);
4013          m= mapUp (m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
4014          G_m= mapUp (G_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
4015          newtonPoly= mapUp (newtonPoly, alpha, V_buf, prim_elem, im_prim_elem,
4016                              source, dest);
4017          ppA= mapUp (ppA, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
4018          ppB= mapUp (ppB, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
4019          gcdlcAlcB= mapUp (gcdlcAlcB, alpha, V_buf, prim_elem, im_prim_elem,
4020                            source, dest);
4021          fail= false;
4022          random_element= randomElement (m, V_buf, l, fail );
4023          DEBOUTLN (cerr, "fail= " << fail);
4024          CFList list;
4025          TIMING_START (gcd_recursion);
4026          if (LC (skeleton).inCoeffDomain())
4027            G_random_element=
4028            monicSparseInterpol (ppA (random_element, x), ppB (random_element, x),
4029                                skeleton, V_buf, sparseFail, coeffMonoms,
4030                                Monoms);
4031          else
4032            G_random_element=
4033            nonMonicSparseInterpol (ppA(random_element,x), ppB(random_element,x),
4034                                    skeleton, V_buf, sparseFail,
4035                                    coeffMonoms, Monoms);
4036          TIMING_END_AND_PRINT (gcd_recursion,
4037                                "time for recursive call: ");
4038          DEBOUTLN (cerr, "G_random_element= " << G_random_element);
4039        }
4040
4041        if (sparseFail)
4042          break;
4043
4044        if (!G_random_element.inCoeffDomain())
4045          d0= totaldegree (G_random_element, Variable(2),
4046                           Variable (G_random_element.level()));
4047        else
4048          d0= 0;
4049
4050        if (d0 == 0)
4051          return N(gcdcAcB);
4052        if (d0 >  d)
4053        {
4054          if (!find (l, random_element))
4055            l.append (random_element);
4056          continue;
4057        }
4058
4059        G_random_element=
4060        (gcdlcAlcB(random_element, x)/uni_lcoeff (G_random_element))
4061        * G_random_element;
4062
4063        if (!G_random_element.inCoeffDomain())
4064          d0= totaldegree (G_random_element, Variable(2),
4065                           Variable (G_random_element.level()));
4066        else
4067          d0= 0;
4068
4069        if (d0 <  d)
4070        {
4071          m= gcdlcAlcB;
4072          newtonPoly= 1;
4073          G_m= 0;
4074          d= d0;
4075        }
4076
4077        TIMING_START (newton_interpolation);
4078        H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
4079        TIMING_END_AND_PRINT (newton_interpolation,
4080                              "time for newton interpolation: ");
4081
4082        //termination test
4083        if (uni_lcoeff (H) == gcdlcAlcB)
4084        {
4085          cH= uni_content (H);
4086          ppH= H/cH;
4087          ppH /= Lc (ppH);
4088          DEBOUTLN (cerr, "ppH= " << ppH);
4089          if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
4090            return N(gcdcAcB*ppH);
4091        }
4092
4093        G_m= H;
4094        newtonPoly= newtonPoly*(x - random_element);
4095        m= m*(x - random_element);
4096        if (!find (l, random_element))
4097          l.append (random_element);
4098
4099      } while (1); //end of second do
4100    }
4101    else
4102      return N(gcdcAcB*GCD_small_p (ppA, ppB));
4103  } while (1); //end of first do
4104}
4105
4106static inline
4107int compress4EZGCD (const CanonicalForm& F, const CanonicalForm& G, CFMap & M,
4108                    CFMap & N, int& both_non_zero)
4109{
4110  int n= tmax (F.level(), G.level());
4111  int * degsf= new int [n + 1];
4112  int * degsg= new int [n + 1];
4113
4114  for (int i = 0; i <= n; i++)
4115    degsf[i]= degsg[i]= 0;
4116
4117  degsf= degrees (F, degsf);
4118  degsg= degrees (G, degsg);
4119
4120  both_non_zero= 0;
4121  int f_zero= 0;
4122  int g_zero= 0;
4123
4124  for (int i= 1; i <= n; i++)
4125  {
4126    if (degsf[i] != 0 && degsg[i] != 0)
4127    {
4128      both_non_zero++;
4129      continue;
4130    }
4131    if (degsf[i] == 0 && degsg[i] != 0 && i <= G.level())
4132    {
4133      f_zero++;
4134      continue;
4135    }
4136    if (degsg[i] == 0 && degsf[i] && i <= F.level())
4137    {
4138      g_zero++;
4139      continue;
4140    }
4141  }
4142
4143  if (both_non_zero == 0)
4144  {
4145    delete [] degsf;
4146    delete [] degsg;
4147    return 0;
4148  }
4149
4150  // map Variables which do not occur in both polynomials to higher levels
4151  int k= 1;
4152  int l= 1;
4153  for (int i= 1; i <= n; i++)
4154  {
4155    if (degsf[i] != 0 && degsg[i] == 0 && i <= F.level())
4156    {
4157      if (k + both_non_zero != i)
4158      {
4159        M.newpair (Variable (i), Variable (k + both_non_zero));
4160        N.newpair (Variable (k + both_non_zero), Variable (i));
4161      }
4162      k++;
4163    }
4164    if (degsf[i] == 0 && degsg[i] != 0 && i <= G.level())
4165    {
4166      if (l + g_zero + both_non_zero != i)
4167      {
4168        M.newpair (Variable (i), Variable (l + g_zero + both_non_zero));
4169        N.newpair (Variable (l + g_zero + both_non_zero), Variable (i));
4170      }
4171      l++;
4172    }
4173  }
4174
4175  // sort Variables x_{i} in decreasing order of
4176  // min(deg_{x_{i}}(f),deg_{x_{i}}(g))
4177  int m= tmin (F.level(), G.level());
4178  int max_min_deg;
4179  k= both_non_zero;
4180  l= 0;
4181  int i= 1;
4182  while (k > 0)
4183  {
4184    max_min_deg= tmin (degsf[i], degsg[i]);
4185    while (max_min_deg == 0)
4186    {
4187      i++;
4188      max_min_deg= tmin (degsf[i], degsg[i]);
4189    }
4190    for (int j= i + 1; j <= m; j++)
4191    {
4192      if ((tmin (degsf[j],degsg[j]) < max_min_deg) &&
4193          (tmin (degsf[j], degsg[j]) != 0))
4194      {
4195        max_min_deg= tmin (degsf[j], degsg[j]);
4196        l= j;
4197      }
4198    }
4199
4200    if (l != 0)
4201    {
4202      if (l != k)
4203      {
4204        M.newpair (Variable (l), Variable(k));
4205        N.newpair (Variable (k), Variable(l));
4206        degsf[l]= 0;
4207        degsg[l]= 0;
4208        l= 0;
4209      }
4210      else
4211      {
4212        degsf[l]= 0;
4213        degsg[l]= 0;
4214        l= 0;
4215      }
4216    }
4217    else if (l == 0)
4218    {
4219      if (i != k)
4220      {
4221        M.newpair (Variable (i), Variable (k));
4222        N.newpair (Variable (k), Variable (i));
4223        degsf[i]= 0;
4224        degsg[i]= 0;
4225      }
4226      else
4227      {
4228        degsf[i]= 0;
4229        degsg[i]= 0;
4230      }
4231      i++;
4232    }
4233    k--;
4234  }
4235
4236  delete [] degsf;
4237  delete [] degsg;
4238
4239  return both_non_zero;
4240}
4241
4242static inline
4243CanonicalForm myShift2Zero (const CanonicalForm& F, CFList& Feval,
4244                            const CFList& evaluation)
4245{
4246  CanonicalForm A= F;
4247  int k= 2;
4248  for (CFListIterator i= evaluation; i.hasItem(); i++, k++)
4249    A= A (Variable (k) + i.getItem(), k);
4250
4251  CanonicalForm buf= A;
4252  Feval= CFList();
4253  Feval.append (buf);
4254  for (k= evaluation.length() + 1; k > 2; k--)
4255  {
4256    buf= mod (buf, Variable (k));
4257    Feval.insert (buf);
4258  }
4259  return A;
4260}
4261
4262static inline
4263CanonicalForm myReverseShift (const CanonicalForm& F, const CFList& evaluation)
4264{
4265  int l= evaluation.length() + 1;
4266  CanonicalForm result= F;
4267  CFListIterator j= evaluation;
4268  for (int i= 2; i < l + 1; i++, j++)
4269  {
4270    if (F.level() < i)
4271      continue;
4272    result= result (Variable (i) - j.getItem(), i);
4273  }
4274  return result;
4275}
4276
4277static inline
4278Evaluation optimize4Lift (const CanonicalForm& F, CFMap & M,
4279                    CFMap & N, const Evaluation& A)
4280{
4281  int n= F.level();
4282  int * degsf= new int [n + 1];
4283
4284  for (int i = 0; i <= n; i++)
4285    degsf[i]= 0;
4286
4287  degsf= degrees (F, degsf);
4288
4289  Evaluation result= Evaluation (A.min(), A.max());
4290  ASSERT (A.min() == 2, "expected A.min() == 2");
4291  int max_deg;
4292  int k= n;
4293  int l= 1;
4294  int i= 2;
4295  int pos= 2;
4296  while (k > 1)
4297  {
4298    max_deg= degsf [i];
4299    while (max_deg == 0)
4300    {
4301      i++;
4302      max_deg= degsf [i];
4303    }
4304    l= i;
4305    for (int j= i + 1; j <=  n; j++)
4306    {
4307      if (degsf[j] > max_deg)
4308      {
4309        max_deg= degsf[j];
4310        l= j;
4311      }
4312    }
4313
4314    if (l <= n)
4315    {
4316      if (l != pos)
4317      {
4318        result.setValue (pos, A [l]);
4319        M.newpair (Variable (l), Variable (pos));
4320        N.newpair (Variable (pos), Variable (l));
4321        degsf[l]= 0;
4322        l= 2;
4323        if (k == 2 && n == 3)
4324        {
4325          result.setValue (l, A [pos]);
4326          M.newpair (Variable (pos), Variable (l));
4327          N.newpair (Variable (l), Variable (pos));
4328          degsf[pos]= 0;
4329        }
4330      }
4331      else
4332      {
4333        result.setValue (l, A [l]);
4334        degsf [l]= 0;
4335      }
4336    }
4337    pos++;
4338    k--;
4339    l= 2;
4340  }
4341
4342  delete [] degsf;
4343
4344  return result;
4345}
4346
4347static inline
4348int Hensel_P (const CanonicalForm & UU, CFArray & G, const Evaluation & AA,
4349              const CFArray& LeadCoeffs )
4350{
4351  CFList factors;
4352  factors.append (G[1]);
4353  factors.append (G[2]);
4354
4355  CFMap NN, MM;
4356  Evaluation A= optimize4Lift (UU, MM, NN, AA);
4357
4358  CanonicalForm U= MM (UU);
4359  CFArray LCs= CFArray (1,2);
4360  LCs [1]= MM (LeadCoeffs [1]);
4361  LCs [2]= MM (LeadCoeffs [2]);
4362
4363  CFList evaluation;
4364  long termEstimate= size (U);
4365  for (int i= A.min(); i <= A.max(); i++)
4366  {
4367    if (!A[i].isZero() && (getCharacteristic() > degree (U,i))) //TODO find a good estimate for getCharacteristic() <= degree (U,i)
4368    {
4369      termEstimate *= degree (U,i)*2;
4370      termEstimate /= 3;
4371    }
4372    evaluation.append (A [i]);
4373  }
4374  if (termEstimate/getNumVars(U) > 500)
4375    return -1;
4376  CFList UEval;
4377  CanonicalForm shiftedU= myShift2Zero (U, UEval, evaluation);
4378
4379  if (size (shiftedU)/getNumVars (U) > 500)
4380    return -1;
4381
4382  CFArray shiftedLCs= CFArray (2);
4383  CFList shiftedLCsEval1, shiftedLCsEval2;
4384  shiftedLCs[0]= myShift2Zero (LCs[1], shiftedLCsEval1, evaluation);
4385  shiftedLCs[1]= myShift2Zero (LCs[2], shiftedLCsEval2, evaluation);
4386  factors.insert (1);
4387  int liftBound= degree (UEval.getLast(), 2) + 1;
4388  CFArray Pi;
4389  CFMatrix M= CFMatrix (liftBound, factors.length() - 1);
4390  CFList diophant;
4391  CFArray lcs= CFArray (2);
4392  lcs [0]= shiftedLCsEval1.getFirst();
4393  lcs [1]= shiftedLCsEval2.getFirst();
4394  nonMonicHenselLift12 (UEval.getFirst(), factors, liftBound, Pi, diophant, M,
4395                        lcs, false);
4396
4397  for (CFListIterator i= factors; i.hasItem(); i++)
4398  {
4399    if (!fdivides (i.getItem(), UEval.getFirst()))
4400      return 0;
4401  }
4402
4403  int * liftBounds;
4404  bool noOneToOne= false;
4405  if (U.level() > 2)
4406  {
4407    liftBounds= new int [U.level() - 1]; /* index: 0.. U.level()-2 */
4408    liftBounds[0]= liftBound;
4409    for (int i= 1; i < U.level() - 1; i++)
4410      liftBounds[i]= degree (shiftedU, Variable (i + 2)) + 1;
4411    factors= nonMonicHenselLift2 (UEval, factors, liftBounds, U.level() - 1,
4412                                  false, shiftedLCsEval1, shiftedLCsEval2, Pi,
4413                                  diophant, noOneToOne);
4414    delete [] liftBounds;
4415    if (noOneToOne)
4416      return 0;
4417  }
4418  G[1]= factors.getFirst();
4419  G[2]= factors.getLast();
4420  G[1]= myReverseShift (G[1], evaluation);
4421  G[2]= myReverseShift (G[2], evaluation);
4422  G[1]= NN (G[1]);
4423  G[2]= NN (G[2]);
4424  return 1;
4425}
4426
4427static inline
4428bool findeval_P (const CanonicalForm & F, const CanonicalForm & G,
4429                 CanonicalForm & Fb, CanonicalForm & Gb, CanonicalForm & Db,
4430                 REvaluation & b, int delta, int degF, int degG, int maxeval,
4431                 int & count, int& k, int bound, int& l)
4432{
4433  if( count == 0 && delta != 0)
4434  {
4435    if( count++ > maxeval )
4436      return false;
4437  }
4438  if (count > 0)
4439  {
4440    b.nextpoint(k);
4441    if (k == 0)
4442      k++;
4443    l++;
4444    if (l > bound)
4445    {
4446      l= 1;
4447      k++;
4448      if (k > tmax (F.level(), G.level()) - 1)
4449        return false;
4450      b.nextpoint (k);
4451    }
4452    if (count++ > maxeval)
4453      return false;
4454  }
4455  while( true )
4456  {
4457    Fb = b( F );
4458    if( degree( Fb, 1 ) == degF )
4459    {
4460      Gb = b( G );
4461      if( degree( Gb, 1 ) == degG )
4462      {
4463        Db = gcd( Fb, Gb );
4464        if( delta > 0 )
4465        {
4466          if( degree( Db, 1 ) <= delta )
4467            return true;
4468        }
4469        else
4470          return true;
4471      }
4472    }
4473    if (k == 0)
4474      k++;
4475    b.nextpoint(k);
4476    l++;
4477    if (l > bound)
4478    {
4479      l= 1;
4480      k++;
4481      if (k > tmax (F.level(), G.level()) - 1)
4482        return false;
4483      b.nextpoint (k);
4484    }
4485    if( count++ > maxeval )
4486      return false;
4487  }
4488}
4489
4490// parameters for heuristic
4491static int maxNumEval= 200;
4492static int sizePerVars1= 500; //try dense gcd if size/#variables is bigger
4493
4494/// Extended Zassenhaus GCD for finite fields
4495CanonicalForm EZGCD_P( const CanonicalForm & FF, const CanonicalForm & GG )
4496{
4497  if (FF.isZero() && degree(GG) > 0) return GG/Lc(GG);
4498  else if (GG.isZero() && degree (FF) > 0) return FF/Lc(FF);
4499  else if (FF.isZero() && GG.isZero()) return FF.genOne();
4500  if (FF.inBaseDomain() || GG.inBaseDomain()) return FF.genOne();
4501  if (FF.isUnivariate() && fdivides(FF, GG)) return FF/Lc(FF);
4502  if (GG.isUnivariate() && fdivides(GG, FF)) return GG/Lc(GG);
4503  if (FF == GG) return FF/Lc(FF);
4504
4505  CanonicalForm F, G, f, g, d, Fb, Gb, Db, Fbt, Gbt, Dbt, B0, B, D0, lcF, lcG,
4506                lcD;
4507  CFArray DD( 1, 2 ), lcDD( 1, 2 );
4508  int degF, degG, delta, count;
4509  int maxeval;
4510  maxeval= tmin((getCharacteristic()/
4511                (int)(ilog2(getCharacteristic())*log2exp))*2, maxNumEval);
4512  count= 0; // number of eval. used
4513  REvaluation b, bt;
4514  int gcdfound = 0;
4515  Variable x = Variable(1);
4516
4517  F= FF;
4518  G= GG;
4519
4520  CFMap M,N;
4521  int smallestDegLev;
4522  TIMING_START (ez_p_compress)
4523  int best_level= compress4EZGCD (F, G, M, N, smallestDegLev);
4524
4525  if (best_level == 0) return G.genOne();
4526
4527  F= M (F);
4528  G= M (G);
4529  TIMING_END_AND_PRINT (ez_p_compress, "time for compression in EZ_P: ")
4530
4531  TIMING_START (ez_p_content)
4532  f = content( F, x ); g = content( G, x );
4533  d = gcd( f, g );
4534  F /= f; G /= g;
4535  TIMING_END_AND_PRINT (ez_p_content, "time to extract content in EZ_P: ")
4536
4537  if( F.isUnivariate() && G.isUnivariate() )
4538  {
4539    if( F.mvar() == G.mvar() )
4540      d *= gcd( F, G );
4541    else
4542      return N (d);
4543    return N (d);
4544  }
4545  if ( F.isUnivariate())
4546  {
4547    g= content (G,G.mvar());
4548    return N(d*gcd(F,g));
4549  }
4550  if ( G.isUnivariate())
4551  {
4552    f= content (F,F.mvar());
4553    return N(d*gcd(G,f));
4554  }
4555
4556  int maxNumVars= tmax (getNumVars (F), getNumVars (G));
4557  Variable a, oldA;
4558  int sizeF= size (F);
4559  int sizeG= size (G);
4560
4561  if (sizeF/maxNumVars > sizePerVars1 && sizeG/maxNumVars > sizePerVars1)
4562  {
4563    if (hasFirstAlgVar (F, a) || hasFirstAlgVar (G, a))
4564      return N (d*GCD_Fp_extension (F, G, a));
4565    else if (CFFactory::gettype() == GaloisFieldDomain)
4566      return N (d*GCD_GF (F, G));
4567    else
4568      return N (d*GCD_small_p (F, G));
4569  }
4570
4571  int dummy= 0;
4572  if( gcd_test_one( F, G, false, dummy ) )
4573  {
4574    return N (d);
4575  }
4576
4577  bool passToGF= false;
4578  bool extOfExt= false;
4579  int p= getCharacteristic();
4580  bool algExtension= (hasFirstAlgVar(F,a) || hasFirstAlgVar(G,a));
4581  int k= 1;
4582  CanonicalForm primElem, imPrimElem;
4583  CFList source, dest;
4584  if (p < 50 && CFFactory::gettype() != GaloisFieldDomain && !algExtension)
4585  {
4586    if (p == 2)
4587      setCharacteristic (2, 12, 'Z');
4588    else if (p == 3)
4589      setCharacteristic (3, 4, 'Z');
4590    else if (p == 5 || p == 7)
4591      setCharacteristic (p, 3, 'Z');
4592    else
4593      setCharacteristic (p, 2, 'Z');
4594    passToGF= true;
4595    F= F.mapinto();
4596    G= G.mapinto();
4597    maxeval= 2*ipower (p, getGFDegree());
4598  }
4599  else if (CFFactory::gettype() == GaloisFieldDomain &&
4600           ipower (p , getGFDegree()) < 50)
4601  {
4602    k= getGFDegree();
4603    if (ipower (p, 2*k) > 50)
4604      setCharacteristic (p, 2*k, gf_name);
4605    else
4606      setCharacteristic (p, 3*k, gf_name);
4607    F= GFMapUp (F, k);
4608    G= GFMapUp (G, k);
4609    maxeval= tmin (2*ipower (p, getGFDegree()), maxNumEval);
4610  }
4611  else if (p < 50 && algExtension && CFFactory::gettype() != GaloisFieldDomain)
4612  {
4613    int d= degree (getMipo (a));
4614    oldA= a;
4615    Variable v2;
4616    if (p == 2 && d < 6)
4617    {
4618      if (fac_NTL_char != p)
4619      {
4620        fac_NTL_char= p;
4621        zz_p::init (p);
4622      }
4623      bool primFail= false;
4624      Variable vBuf;
4625      primElem= primitiveElement (a, vBuf, primFail);
4626      ASSERT (!primFail, "failure in integer factorizer");
4627      if (d < 3)
4628      {
4629        zz_pX NTLIrredpoly;
4630        BuildIrred (NTLIrredpoly, d*3);
4631        CanonicalForm newMipo= convertNTLzzpX2CF (NTLIrredpoly, Variable (1));
4632        v2= rootOf (newMipo);
4633      }
4634      else
4635      {
4636        zz_pX NTLIrredpoly;
4637        BuildIrred (NTLIrredpoly, d*2);
4638        CanonicalForm newMipo= convertNTLzzpX2CF (NTLIrredpoly, Variable (1));
4639        v2= rootOf (newMipo);
4640      }
4641      imPrimElem= mapPrimElem (primElem, a, v2);
4642      extOfExt= true;
4643    }
4644    else if ((p == 3 && d < 4) || ((p == 5 || p == 7) && d < 3))
4645    {
4646      if (fac_NTL_char != p)
4647      {
4648        fac_NTL_char= p;
4649        zz_p::init (p);
4650      }
4651      bool primFail= false;
4652      Variable vBuf;
4653      primElem= primitiveElement (a, vBuf, primFail);
4654      ASSERT (!primFail, "failure in integer factorizer");
4655      zz_pX NTLIrredpoly;
4656      BuildIrred (NTLIrredpoly, d*2);
4657      CanonicalForm newMipo= convertNTLzzpX2CF (NTLIrredpoly, Variable (1));
4658      v2= rootOf (newMipo);
4659      imPrimElem= mapPrimElem (primElem, a, v2);
4660      extOfExt= true;
4661    }
4662    if (extOfExt)
4663    {
4664      maxeval= tmin (2*ipower (p, degree (getMipo (v2))), maxNumEval);
4665      F= mapUp (F, a, v2, primElem, imPrimElem, source, dest);
4666      G= mapUp (G, a, v2, primElem, imPrimElem, source, dest);
4667      a= v2;
4668    }
4669  }
4670
4671  lcF = LC( F, x ); lcG = LC( G, x );
4672  lcD = gcd( lcF, lcG );
4673
4674  delta = 0;
4675  degF = degree( F, x ); degG = degree( G, x );
4676
4677  if(hasFirstAlgVar(G,a))
4678    b = REvaluation( 2, tmax(F.level(), G.level()), AlgExtRandomF( a ) );
4679  else
4680  { // both not in extension given by algebraic variable
4681    if (CFFactory::gettype() != GaloisFieldDomain)
4682      b = REvaluation( 2, tmax(F.level(), G.level()), FFRandom() );
4683    else
4684      b = REvaluation( 2, tmax(F.level(), G.level()), GFRandom() );
4685  }
4686
4687  CanonicalForm cand, contcand;
4688  CanonicalForm result;
4689  int o, t;
4690  o= 0;
4691  t= 1;
4692  int goodPointCount= 0;
4693  while( !gcdfound )
4694  {
4695    TIMING_START (ez_p_eval);
4696    if( !findeval_P( F, G, Fb, Gb, Db, b, delta, degF, degG, maxeval, count, o,
4697         maxeval/maxNumVars, t ))
4698    { // too many eval. used --> try another method
4699      Off (SW_USE_EZGCD_P);
4700      result= gcd (F,G);
4701      On (SW_USE_EZGCD_P);
4702      if (passToGF)
4703      {
4704        CanonicalForm mipo= gf_mipo;
4705        setCharacteristic (p);
4706        Variable alpha= rootOf (mipo.mapinto());
4707        result= GF2FalphaRep (result, alpha);
4708      }
4709      if (k > 1)
4710      {
4711        result= GFMapDown (result, k);
4712        setCharacteristic (p, k, gf_name);
4713      }
4714      if (extOfExt)
4715        result= mapDown (result, primElem, imPrimElem, oldA, dest, source);
4716      return N (d*result);
4717    }
4718    TIMING_END_AND_PRINT (ez_p_eval, "time for eval point search in EZ_P1: ");
4719    delta = degree( Db );
4720    if( delta == 0 )
4721    {
4722      if (passToGF)
4723        setCharacteristic (p);
4724      if (k > 1)
4725        setCharacteristic (p, k, gf_name);
4726      return N (d);
4727    }
4728    while( true )
4729    {
4730      bt = b;
4731      TIMING_START (ez_p_eval);
4732      if( !findeval_P(F,G,Fbt,Gbt,Dbt, bt, delta, degF, degG, maxeval, count, o,
4733           maxeval/maxNumVars, t ))
4734      { // too many eval. used --> try another method
4735        Off (SW_USE_EZGCD_P);
4736        result= gcd (F,G);
4737        On (SW_USE_EZGCD_P);
4738        if (passToGF)
4739        {
4740          CanonicalForm mipo= gf_mipo;
4741          setCharacteristic (p);
4742          Variable alpha= rootOf (mipo.mapinto());
4743          result= GF2FalphaRep (result, alpha);
4744        }
4745        if (k > 1)
4746        {
4747          result= GFMapDown (result, k);
4748          setCharacteristic (p, k, gf_name);
4749        }
4750        if (extOfExt)
4751          result= mapDown (result, primElem, imPrimElem, oldA, dest, source);
4752        return N (d*result);
4753      }
4754      TIMING_END_AND_PRINT (ez_p_eval, "time for eval point search in EZ_P2: ");
4755      int dd = degree( Dbt );
4756      if( dd == 0 )
4757      {
4758        if (passToGF)
4759          setCharacteristic (p);
4760        if (k > 1)
4761          setCharacteristic (p, k, gf_name);
4762        return N (d);
4763      }
4764      if( dd == delta )
4765      {
4766        goodPointCount++;
4767        if (goodPointCount == 5)
4768          break;
4769      }
4770      if( dd < delta )
4771      {
4772        goodPointCount= 0;
4773        delta = dd;
4774        b = bt;
4775        Db = Dbt; Fb = Fbt; Gb = Gbt;
4776      }
4777      if (delta == degF)
4778      {
4779        if (degF <= degG && fdivides (F, G))
4780        {
4781          if (passToGF)
4782          {
4783            CanonicalForm mipo= gf_mipo;
4784            setCharacteristic (p);
4785            Variable alpha= rootOf (mipo.mapinto());
4786            F= GF2FalphaRep (F, alpha);
4787          }
4788          if (k > 1)
4789          {
4790            F= GFMapDown (F, k);
4791            setCharacteristic (p, k, gf_name);
4792          }
4793          if (extOfExt)
4794            F= mapDown (F, primElem, imPrimElem, oldA, dest, source);
4795          return N (d*F);
4796        }
4797        else
4798          delta--;
4799      }
4800      else if (delta == degG)
4801      {
4802        if (degG <= degF && fdivides (G, F))
4803        {
4804          if (passToGF)
4805          {
4806            CanonicalForm mipo= gf_mipo;
4807            setCharacteristic (p);
4808            Variable alpha= rootOf (mipo.mapinto());
4809            G= GF2FalphaRep (G, alpha);
4810          }
4811          if (k > 1)
4812          {
4813            G= GFMapDown (G, k);
4814            setCharacteristic (p, k, gf_name);
4815          }
4816          if (extOfExt)
4817            G= mapDown (G, primElem, imPrimElem, oldA, dest, source);
4818          return N (d*G);
4819        }
4820        else
4821          delta--;
4822      }
4823    }
4824    if( delta != degF && delta != degG )
4825    {
4826      bool B_is_F;
4827      CanonicalForm xxx1, xxx2;
4828      CanonicalForm buf;
4829      DD[1] = Fb / Db;
4830      buf= Gb/Db;
4831      xxx1 = gcd( DD[1], Db );
4832      xxx2 = gcd( buf, Db );
4833      if (((xxx1.inCoeffDomain() && xxx2.inCoeffDomain()) &&
4834          (size (F) <= size (G)))
4835          || (xxx1.inCoeffDomain() && !xxx2.inCoeffDomain()))
4836      {
4837        B = F;
4838        DD[2] = Db;
4839        lcDD[1] = lcF;
4840        lcDD[2] = lcD;
4841        B_is_F = true;
4842      }
4843      else if (((xxx1.inCoeffDomain() && xxx2.inCoeffDomain()) &&
4844               (size (G) < size (F)))
4845               || (!xxx1.inCoeffDomain() && xxx2.inCoeffDomain()))
4846      {
4847        DD[1] = buf;
4848        B = G;
4849        DD[2] = Db;
4850        lcDD[1] = lcG;
4851        lcDD[2] = lcD;
4852        B_is_F = false;
4853      }
4854      else // special case handling
4855      {
4856        Off (SW_USE_EZGCD_P);
4857        result= gcd (F,G);
4858        On (SW_USE_EZGCD_P);
4859        if (passToGF)
4860        {
4861          CanonicalForm mipo= gf_mipo;
4862          setCharacteristic (p);
4863          Variable alpha= rootOf (mipo.mapinto());
4864          result= GF2FalphaRep (result, alpha);
4865        }
4866        if (k > 1)
4867        {
4868          result= GFMapDown (result, k);
4869          setCharacteristic (p, k, gf_name);
4870        }
4871        if (extOfExt)
4872          result= mapDown (result, primElem, imPrimElem, oldA, dest, source);
4873        return N (d*result);
4874      }
4875      DD[2] = DD[2] * ( b( lcDD[2] ) / lc( DD[2] ) );
4876      DD[1] = DD[1] * ( b( lcDD[1] ) / lc( DD[1] ) );
4877
4878      if (size (B*lcDD[2])/maxNumVars > sizePerVars1)
4879      {
4880        if (algExtension)
4881        {
4882          result= GCD_Fp_extension (F, G, a);
4883          if (extOfExt)
4884            result= mapDown (result, primElem, imPrimElem, oldA, dest, source);
4885          return N (d*result);
4886        }
4887        if (CFFactory::gettype() == GaloisFieldDomain)
4888        {
4889          result= GCD_GF (F, G);
4890          if (passToGF)
4891          {
4892            CanonicalForm mipo= gf_mipo;
4893            setCharacteristic (p);
4894            Variable alpha= rootOf (mipo.mapinto());
4895            result= GF2FalphaRep (result, alpha);
4896          }
4897          if (k > 1)
4898          {
4899            result= GFMapDown (result, k);
4900            setCharacteristic (p, k, gf_name);
4901          }
4902          return N (d*result);
4903        }
4904        else
4905          return N (d*GCD_small_p (F,G));
4906      }
4907
4908      TIMING_START (ez_p_hensel_lift);
4909      gcdfound= Hensel_P (B*lcD, DD, b, lcDD);
4910      TIMING_END_AND_PRINT (ez_p_hensel_lift, "time for Hensel lift in EZ_P: ");
4911
4912      if (gcdfound == -1) //things became dense
4913      {
4914        if (algExtension)
4915        {
4916          result= GCD_Fp_extension (F, G, a);
4917          if (extOfExt)
4918            result= mapDown (result, primElem, imPrimElem, oldA, dest, source);
4919          return N (d*result);
4920        }
4921        if (CFFactory::gettype() == GaloisFieldDomain)
4922        {
4923          result= GCD_GF (F, G);
4924          if (passToGF)
4925          {
4926            CanonicalForm mipo= gf_mipo;
4927            setCharacteristic (p);
4928            Variable alpha= rootOf (mipo.mapinto());
4929            result= GF2FalphaRep (result, alpha);
4930          }
4931          if (k > 1)
4932          {
4933            result= GFMapDown (result, k);
4934            setCharacteristic (p, k, gf_name);
4935          }
4936          return N (d*result);
4937        }
4938        else
4939        {
4940          if (p >= cf_getBigPrime(0))
4941            return N (d*sparseGCDFp (F,G));
4942          else
4943            return N (d*GCD_small_p (F,G));
4944        }
4945      }
4946
4947      if (gcdfound == 1)
4948      {
4949        TIMING_START (termination_test);
4950        contcand= content (DD[2], Variable (1));
4951        cand = DD[2] / contcand;
4952        if (B_is_F)
4953          gcdfound = fdivides( cand, G ) && cand*(DD[1]/(lcD/contcand)) == F;
4954        else
4955          gcdfound = fdivides( cand, F ) && cand*(DD[1]/(lcD/contcand)) == G;
4956        TIMING_END_AND_PRINT (termination_test,
4957                              "time for termination test EZ_P: ");
4958
4959        if (passToGF && gcdfound)
4960        {
4961          CanonicalForm mipo= gf_mipo;
4962          setCharacteristic (p);
4963          Variable alpha= rootOf (mipo.mapinto());
4964          cand= GF2FalphaRep (cand, alpha);
4965        }
4966        if (k > 1 && gcdfound)
4967        {
4968          cand= GFMapDown (cand, k);
4969          setCharacteristic (p, k, gf_name);
4970        }
4971        if (extOfExt && gcdfound)
4972          cand= mapDown (cand, primElem, imPrimElem, oldA, dest, source);
4973      }
4974    }
4975    delta--;
4976    goodPointCount= 0;
4977  }
4978  return N (d*cand);
4979}
4980
4981
4982#endif
Note: See TracBrowser for help on using the repository browser.