source: git/factory/cf_gcd_smallp.cc @ 147b6d

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