source: git/factory/cf_gcd_smallp.cc @ fd80670

spielwiese
Last change on this file since fd80670 was bffe62d, checked in by Martin Lee <martinlee84@…>, 11 years ago
chg: skip zz_p::init() if it is already correctly initialized
  • Property mode set to 100644
File size: 133.4 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* dummy= new long [M.rows()];
2003  for (int i= 0; i < M.rows(); i++)
2004    dummy[i]= 0;
2005  long rk= nmod_mat_rref (dummy, FLINTN);
2006
2007  N= convertNmod_mat_t2FacCFMatrix (FLINTN);
2008  nmod_mat_clear (FLINTN);
2009  delete dummy;
2010#else
2011  int p= getCharacteristic ();
2012  if (fac_NTL_char != p)
2013  {
2014    fac_NTL_char= p;
2015    zz_p::init (p);
2016  }
2017  mat_zz_p *NTLN= convertFacCFMatrix2NTLmat_zz_p(*N);
2018  long rk= gauss (*NTLN);
2019
2020  N= convertNTLmat_zz_p2FacCFMatrix (*NTLN);
2021#endif
2022
2023  L= CFArray (M.rows());
2024  for (int i= 0; i < M.rows(); i++)
2025    L[i]= (*N) (i + 1, M.columns() + 1);
2026  M= (*N) (1, M.rows(), 1, M.columns());
2027  delete N;
2028  return rk;
2029}
2030
2031long
2032gaussianElimFq (CFMatrix& M, CFArray& L, const Variable& alpha)
2033{
2034  ASSERT (L.size() <= M.rows(), "dimension exceeded");
2035  CFMatrix *N;
2036  N= new CFMatrix (M.rows(), M.columns() + 1);
2037
2038  for (int i= 1; i <= M.rows(); i++)
2039    for (int j= 1; j <= M.columns(); j++)
2040      (*N) (i, j)= M (i, j);
2041
2042  int j= 1;
2043  for (int i= 0; i < L.size(); i++, j++)
2044    (*N) (j, M.columns() + 1)= L[i];
2045  int p= getCharacteristic ();
2046  if (fac_NTL_char != p)
2047  {
2048    fac_NTL_char= p;
2049    zz_p::init (p);
2050  }
2051  zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
2052  zz_pE::init (NTLMipo);
2053  mat_zz_pE *NTLN= convertFacCFMatrix2NTLmat_zz_pE(*N);
2054  long rk= gauss (*NTLN);
2055
2056  N= convertNTLmat_zz_pE2FacCFMatrix (*NTLN, alpha);
2057
2058  M= (*N) (1, M.rows(), 1, M.columns());
2059  L= CFArray (M.rows());
2060  for (int i= 0; i < M.rows(); i++)
2061    L[i]= (*N) (i + 1, M.columns() + 1);
2062
2063  delete N;
2064  return rk;
2065}
2066
2067CFArray
2068solveSystemFp (const CFMatrix& M, const CFArray& L)
2069{
2070  ASSERT (L.size() <= M.rows(), "dimension exceeded");
2071  CFMatrix *N;
2072  N= new CFMatrix (M.rows(), M.columns() + 1);
2073
2074  for (int i= 1; i <= M.rows(); i++)
2075    for (int j= 1; j <= M.columns(); j++)
2076      (*N) (i, j)= M (i, j);
2077
2078  int j= 1;
2079  for (int i= 0; i < L.size(); i++, j++)
2080    (*N) (j, M.columns() + 1)= L[i];
2081
2082#ifdef HAVE_FLINT
2083  nmod_mat_t FLINTN;
2084  convertFacCFMatrix2nmod_mat_t (FLINTN, *N);
2085  long* dummy= new long [M.rows()];
2086  for (int i= 0; i < M.rows(); i++)
2087    dummy[i]= 0;
2088  long rk= nmod_mat_rref (dummy, FLINTN);
2089#else
2090  int p= getCharacteristic ();
2091  if (fac_NTL_char != p)
2092  {
2093    fac_NTL_char= p;
2094    zz_p::init (p);
2095  }
2096  mat_zz_p *NTLN= convertFacCFMatrix2NTLmat_zz_p(*N);
2097  long rk= gauss (*NTLN);
2098#endif
2099  if (rk != M.columns())
2100  {
2101#ifdef HAVE_FLINT
2102    nmod_mat_clear (FLINTN);
2103    delete dummy;
2104#endif
2105    delete N;
2106    return CFArray();
2107  }
2108#ifdef HAVE_FLINT
2109  N= convertNmod_mat_t2FacCFMatrix (FLINTN);
2110  nmod_mat_clear (FLINTN);
2111  delete dummy;
2112#else
2113  N= convertNTLmat_zz_p2FacCFMatrix (*NTLN);
2114#endif
2115  CFArray A= readOffSolution (*N, rk);
2116
2117  delete N;
2118  return A;
2119}
2120
2121CFArray
2122solveSystemFq (const CFMatrix& M, const CFArray& L, const Variable& alpha)
2123{
2124  ASSERT (L.size() <= M.rows(), "dimension exceeded");
2125  CFMatrix *N;
2126  N= new CFMatrix (M.rows(), M.columns() + 1);
2127
2128  for (int i= 1; i <= M.rows(); i++)
2129    for (int j= 1; j <= M.columns(); j++)
2130      (*N) (i, j)= M (i, j);
2131  int j= 1;
2132  for (int i= 0; i < L.size(); i++, j++)
2133    (*N) (j, M.columns() + 1)= L[i];
2134  int p= getCharacteristic ();
2135  if (fac_NTL_char != p)
2136  {
2137    fac_NTL_char= p;
2138    zz_p::init (p);
2139  }
2140  zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
2141  zz_pE::init (NTLMipo);
2142  mat_zz_pE *NTLN= convertFacCFMatrix2NTLmat_zz_pE(*N);
2143  long rk= gauss (*NTLN);
2144  if (rk != M.columns())
2145  {
2146    delete N;
2147    return CFArray();
2148  }
2149  N= convertNTLmat_zz_pE2FacCFMatrix (*NTLN, alpha);
2150
2151  CFArray A= readOffSolution (*N, rk);
2152
2153  delete N;
2154  return A;
2155}
2156#endif
2157
2158CFArray
2159getMonoms (const CanonicalForm& F)
2160{
2161  if (F.inCoeffDomain())
2162  {
2163    CFArray result= CFArray (1);
2164    result [0]= 1;
2165    return result;
2166  }
2167  if (F.isUnivariate())
2168  {
2169    CFArray result= CFArray (size(F));
2170    int j= 0;
2171    for (CFIterator i= F; i.hasTerms(); i++, j++)
2172      result[j]= power (F.mvar(), i.exp());
2173    return result;
2174  }
2175  int numMon= size (F);
2176  CFArray result= CFArray (numMon);
2177  int j= 0;
2178  CFArray recResult;
2179  Variable x= F.mvar();
2180  CanonicalForm powX;
2181  for (CFIterator i= F; i.hasTerms(); i++)
2182  {
2183    powX= power (x, i.exp());
2184    recResult= getMonoms (i.coeff());
2185    for (int k= 0; k < recResult.size(); k++)
2186      result[j+k]= powX*recResult[k];
2187    j += recResult.size();
2188  }
2189  return result;
2190}
2191
2192#ifdef HAVE_NTL
2193CFArray
2194evaluateMonom (const CanonicalForm& F, const CFList& evalPoints)
2195{
2196  if (F.inCoeffDomain())
2197  {
2198    CFArray result= CFArray (1);
2199    result [0]= F;
2200    return result;
2201  }
2202  if (F.isUnivariate())
2203  {
2204    ASSERT (evalPoints.length() == 1,
2205            "expected an eval point with only one component");
2206    CFArray result= CFArray (size(F));
2207    int j= 0;
2208    CanonicalForm evalPoint= evalPoints.getLast();
2209    for (CFIterator i= F; i.hasTerms(); i++, j++)
2210      result[j]= power (evalPoint, i.exp());
2211    return result;
2212  }
2213  int numMon= size (F);
2214  CFArray result= CFArray (numMon);
2215  int j= 0;
2216  CanonicalForm evalPoint= evalPoints.getLast();
2217  CFList buf= evalPoints;
2218  buf.removeLast();
2219  CFArray recResult;
2220  CanonicalForm powEvalPoint;
2221  for (CFIterator i= F; i.hasTerms(); i++)
2222  {
2223    powEvalPoint= power (evalPoint, i.exp());
2224    recResult= evaluateMonom (i.coeff(), buf);
2225    for (int k= 0; k < recResult.size(); k++)
2226      result[j+k]= powEvalPoint*recResult[k];
2227    j += recResult.size();
2228  }
2229  return result;
2230}
2231
2232CFArray
2233evaluate (const CFArray& A, const CFList& evalPoints)
2234{
2235  CFArray result= A.size();
2236  CanonicalForm tmp;
2237  int k;
2238  for (int i= 0; i < A.size(); i++)
2239  {
2240    tmp= A[i];
2241    k= 1;
2242    for (CFListIterator j= evalPoints; j.hasItem(); j++, k++)
2243      tmp= tmp (j.getItem(), k);
2244    result[i]= tmp;
2245  }
2246  return result;
2247}
2248
2249CFList
2250evaluationPoints (const CanonicalForm& F, const CanonicalForm& G,
2251                  CanonicalForm& Feval, CanonicalForm& Geval,
2252                  const CanonicalForm& LCF, const bool& GF,
2253                  const Variable& alpha, bool& fail, CFList& list
2254                 )
2255{
2256  int k= tmax (F.level(), G.level()) - 1;
2257  Variable x= Variable (1);
2258  CFList result;
2259  FFRandom genFF;
2260  GFRandom genGF;
2261  int p= getCharacteristic ();
2262  int bound;
2263  if (alpha != Variable (1))
2264  {
2265    bound= ipower (p, degree (getMipo(alpha)));
2266    bound= ipower (bound, k);
2267  }
2268  else if (GF)
2269  {
2270    bound= ipower (p, getGFDegree());
2271    bound= ipower (bound, k);
2272  }
2273  else
2274    bound= ipower (p, k);
2275
2276  CanonicalForm random;
2277  int j;
2278  bool zeroOneOccured= false;
2279  bool allEqual= false;
2280  CanonicalForm buf;
2281  do
2282  {
2283    random= 0;
2284    // possible overflow if list.length() does not fit into a int
2285    if (list.length() >= bound)
2286    {
2287      fail= true;
2288      break;
2289    }
2290    for (int i= 0; i < k; i++)
2291    {
2292      if (GF)
2293      {
2294        result.append (genGF.generate());
2295        random += result.getLast()*power (x, i);
2296      }
2297      else if (alpha.level() != 1)
2298      {
2299        AlgExtRandomF genAlgExt (alpha);
2300        result.append (genAlgExt.generate());
2301        random += result.getLast()*power (x, i);
2302      }
2303      else
2304      {
2305        result.append (genFF.generate());
2306        random += result.getLast()*power (x, i);
2307      }
2308      if (result.getLast().isOne() || result.getLast().isZero())
2309        zeroOneOccured= true;
2310    }
2311    if (find (list, random))
2312    {
2313      zeroOneOccured= false;
2314      allEqual= false;
2315      result= CFList();
2316      continue;
2317    }
2318    if (zeroOneOccured)
2319    {
2320      list.append (random);
2321      zeroOneOccured= false;
2322      allEqual= false;
2323      result= CFList();
2324      continue;
2325    }
2326    // no zero at this point
2327    if (k > 1)
2328    {
2329      allEqual= true;
2330      CFIterator iter= random;
2331      buf= iter.coeff();
2332      iter++;
2333      for (; iter.hasTerms(); iter++)
2334        if (buf != iter.coeff())
2335          allEqual= false;
2336    }
2337    if (allEqual)
2338    {
2339      list.append (random);
2340      allEqual= false;
2341      zeroOneOccured= false;
2342      result= CFList();
2343      continue;
2344    }
2345
2346    Feval= F;
2347    Geval= G;
2348    CanonicalForm LCeval= LCF;
2349    j= 1;
2350    for (CFListIterator i= result; i.hasItem(); i++, j++)
2351    {
2352      Feval= Feval (i.getItem(), j);
2353      Geval= Geval (i.getItem(), j);
2354      LCeval= LCeval (i.getItem(), j);
2355    }
2356
2357    if (LCeval.isZero())
2358    {
2359      if (!find (list, random))
2360        list.append (random);
2361      zeroOneOccured= false;
2362      allEqual= false;
2363      result= CFList();
2364      continue;
2365    }
2366
2367    if (list.length() >= bound)
2368    {
2369      fail= true;
2370      break;
2371    }
2372  } while (find (list, random));
2373
2374  return result;
2375}
2376
2377/// multiply two lists componentwise
2378void mult (CFList& L1, const CFList& L2)
2379{
2380  ASSERT (L1.length() == L2.length(), "lists of the same size expected");
2381
2382  CFListIterator j= L2;
2383  for (CFListIterator i= L1; i.hasItem(); i++, j++)
2384    i.getItem() *= j.getItem();
2385}
2386
2387void eval (const CanonicalForm& A, const CanonicalForm& B, CanonicalForm& Aeval,
2388           CanonicalForm& Beval, const CFList& L)
2389{
2390  Aeval= A;
2391  Beval= B;
2392  int j= 1;
2393  for (CFListIterator i= L; i.hasItem(); i++, j++)
2394  {
2395    Aeval= Aeval (i.getItem(), j);
2396    Beval= Beval (i.getItem(), j);
2397  }
2398}
2399
2400CanonicalForm
2401monicSparseInterpol (const CanonicalForm& F, const CanonicalForm& G,
2402                     const CanonicalForm& skeleton, const Variable& alpha,
2403                     bool& fail, CFArray*& coeffMonoms, CFArray& Monoms
2404                    )
2405{
2406  CanonicalForm A= F;
2407  CanonicalForm B= G;
2408  if (F.isZero() && degree(G) > 0) return G/Lc(G);
2409  else if (G.isZero() && degree (F) > 0) return F/Lc(F);
2410  else if (F.isZero() && G.isZero()) return F.genOne();
2411  if (F.inBaseDomain() || G.inBaseDomain()) return F.genOne();
2412  if (F.isUnivariate() && fdivides(F, G)) return F/Lc(F);
2413  if (G.isUnivariate() && fdivides(G, F)) return G/Lc(G);
2414  if (F == G) return F/Lc(F);
2415
2416  CFMap M,N;
2417  int best_level= myCompress (A, B, M, N, false);
2418
2419  if (best_level == 0)
2420    return B.genOne();
2421
2422  A= M(A);
2423  B= M(B);
2424
2425  Variable x= Variable (1);
2426  ASSERT (degree (A, x) == 0, "expected degree (F, 1) == 0");
2427  ASSERT (degree (B, x) == 0, "expected degree (G, 1) == 0");
2428
2429  //univariate case
2430  if (A.isUnivariate() && B.isUnivariate())
2431    return N (gcd (A, B));
2432
2433  CanonicalForm skel= M(skeleton);
2434  CanonicalForm cA, cB;    // content of A and B
2435  CanonicalForm ppA, ppB;    // primitive part of A and B
2436  CanonicalForm gcdcAcB;
2437  cA = uni_content (A);
2438  cB = uni_content (B);
2439  gcdcAcB= gcd (cA, cB);
2440  ppA= A/cA;
2441  ppB= B/cB;
2442
2443  CanonicalForm lcA, lcB;  // leading coefficients of A and B
2444  CanonicalForm gcdlcAlcB;
2445  lcA= uni_lcoeff (ppA);
2446  lcB= uni_lcoeff (ppB);
2447
2448  if (fdivides (lcA, lcB))
2449  {
2450    if (fdivides (A, B))
2451      return F/Lc(F);
2452  }
2453  if (fdivides (lcB, lcA))
2454  {
2455    if (fdivides (B, A))
2456      return G/Lc(G);
2457  }
2458
2459  gcdlcAlcB= gcd (lcA, lcB);
2460  int skelSize= size (skel, skel.mvar());
2461
2462  int j= 0;
2463  int biggestSize= 0;
2464
2465  for (CFIterator i= skel; i.hasTerms(); i++, j++)
2466    biggestSize= tmax (biggestSize, size (i.coeff()));
2467
2468  CanonicalForm g, Aeval, Beval;
2469
2470  CFList evalPoints;
2471  bool evalFail= false;
2472  CFList list;
2473  bool GF= false;
2474  CanonicalForm LCA= LC (A);
2475  CanonicalForm tmp;
2476  CFArray gcds= CFArray (biggestSize);
2477  CFList * pEvalPoints= new CFList [biggestSize];
2478  Variable V_buf= alpha;
2479  CFList source, dest;
2480  CanonicalForm prim_elem, im_prim_elem;
2481  for (int i= 0; i < biggestSize; i++)
2482  {
2483    if (i == 0)
2484      evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf, evalFail,
2485                                    list);
2486    else
2487    {
2488      mult (evalPoints, pEvalPoints [0]);
2489      eval (A, B, Aeval, Beval, evalPoints);
2490    }
2491
2492    if (evalFail)
2493    {
2494      if (V_buf.level() != 1)
2495      {
2496        do
2497        {
2498          Variable V_buf2= chooseExtension (V_buf);
2499          source= CFList();
2500          dest= CFList();
2501
2502          bool prim_fail= false;
2503          Variable V_buf3;
2504          prim_elem= primitiveElement (V_buf, V_buf3, prim_fail);
2505
2506          ASSERT (!prim_fail, "failure in integer factorizer");
2507          if (prim_fail)
2508            ; //ERROR
2509          else
2510            im_prim_elem= mapPrimElem (prim_elem, V_buf, V_buf2);
2511
2512          DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf));
2513          DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf2));
2514
2515          for (CFListIterator j= list; j.hasItem(); j++)
2516            j.getItem()= mapUp (j.getItem(), V_buf, V_buf2, prim_elem,
2517                                im_prim_elem, source, dest);
2518          for (int k= 0; k < i; k++)
2519          {
2520            for (CFListIterator j= pEvalPoints[k]; j.hasItem(); j++)
2521              j.getItem()= mapUp (j.getItem(), V_buf, V_buf2, prim_elem,
2522                                  im_prim_elem, source, dest);
2523            gcds[k]= mapUp (gcds[k], V_buf, V_buf2, prim_elem, im_prim_elem,
2524                            source, dest);
2525          }
2526
2527          if (alpha.level() != 1)
2528          {
2529            A= mapUp (A, V_buf, V_buf2, prim_elem, im_prim_elem, source,dest);
2530            B= mapUp (B, V_buf, V_buf2, prim_elem, im_prim_elem, source,dest);
2531          }
2532          evalFail= false;
2533          evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf,
2534                                        evalFail, list);
2535        } while (evalFail);
2536      }
2537      else
2538      {
2539        CanonicalForm mipo;
2540        int deg= 2;
2541        do {
2542          mipo= randomIrredpoly (deg, x);
2543          V_buf= rootOf (mipo);
2544          evalFail= false;
2545          evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf,
2546                                        evalFail, list);
2547          deg++;
2548        } while (evalFail);
2549      }
2550    }
2551
2552    g= gcd (Aeval, Beval);
2553    g /= Lc (g);
2554
2555    if (degree (g) != degree (skel) || (skelSize != size (g)))
2556    {
2557      delete[] pEvalPoints;
2558      fail= true;
2559      return 0;
2560    }
2561    CFIterator l= skel;
2562    for (CFIterator k= g; k.hasTerms(); k++, l++)
2563    {
2564      if (k.exp() != l.exp())
2565      {
2566        delete[] pEvalPoints;
2567        fail= true;
2568        return 0;
2569      }
2570    }
2571    pEvalPoints[i]= evalPoints;
2572    gcds[i]= g;
2573
2574    tmp= 0;
2575    int j= 0;
2576    for (CFListIterator k= evalPoints; k.hasItem(); k++, j++)
2577      tmp += k.getItem()*power (x, j);
2578    list.append (tmp);
2579  }
2580
2581  if (Monoms.size() == 0)
2582    Monoms= getMonoms (skel);
2583  if (coeffMonoms == NULL)
2584    coeffMonoms= new CFArray [skelSize];
2585  j= 0;
2586  for (CFIterator i= skel; i.hasTerms(); i++, j++)
2587    coeffMonoms[j]= getMonoms (i.coeff());
2588
2589  CFArray* pL= new CFArray [skelSize];
2590  CFArray* pM= new CFArray [skelSize];
2591  for (int i= 0; i < biggestSize; i++)
2592  {
2593    CFIterator l= gcds [i];
2594    evalPoints= pEvalPoints [i];
2595    for (int k= 0; k < skelSize; k++, l++)
2596    {
2597      if (i == 0)
2598        pL[k]= CFArray (biggestSize);
2599      pL[k] [i]= l.coeff();
2600
2601      if (i == 0)
2602        pM[k]= evaluate (coeffMonoms [k], evalPoints);
2603    }
2604  }
2605
2606  CFArray solution;
2607  CanonicalForm result= 0;
2608  int ind= 0;
2609  CFArray bufArray;
2610  CFMatrix Mat;
2611  for (int k= 0; k < skelSize; k++)
2612  {
2613    if (biggestSize != coeffMonoms[k].size())
2614    {
2615      bufArray= CFArray (coeffMonoms[k].size());
2616      for (int i= 0; i < coeffMonoms[k].size(); i++)
2617        bufArray [i]= pL[k] [i];
2618      solution= solveGeneralVandermonde (pM [k], bufArray);
2619    }
2620    else
2621      solution= solveGeneralVandermonde (pM [k], pL[k]);
2622
2623    if (solution.size() == 0)
2624    {
2625      delete[] pEvalPoints;
2626      delete[] pM;
2627      delete[] pL;
2628      delete[] coeffMonoms;
2629      fail= true;
2630      return 0;
2631    }
2632    for (int l= 0; l < solution.size(); l++)
2633      result += solution[l]*Monoms [ind + l];
2634    ind += solution.size();
2635  }
2636
2637  delete[] pEvalPoints;
2638  delete[] pM;
2639  delete[] pL;
2640
2641  if (alpha.level() != 1 && V_buf != alpha)
2642  {
2643    CFList u, v;
2644    result= mapDown (result, prim_elem, im_prim_elem, alpha, u, v);
2645  }
2646
2647  result= N(result);
2648  if (fdivides (result, F) && fdivides (result, G))
2649    return result;
2650  else
2651  {
2652    delete[] coeffMonoms;
2653    fail= true;
2654    return 0;
2655  }
2656}
2657
2658CanonicalForm
2659nonMonicSparseInterpol (const CanonicalForm& F, const CanonicalForm& G,
2660                        const CanonicalForm& skeleton, const Variable& alpha,
2661                        bool& fail, CFArray*& coeffMonoms, CFArray& Monoms
2662                       )
2663{
2664  CanonicalForm A= F;
2665  CanonicalForm B= G;
2666  if (F.isZero() && degree(G) > 0) return G/Lc(G);
2667  else if (G.isZero() && degree (F) > 0) return F/Lc(F);
2668  else if (F.isZero() && G.isZero()) return F.genOne();
2669  if (F.inBaseDomain() || G.inBaseDomain()) return F.genOne();
2670  if (F.isUnivariate() && fdivides(F, G)) return F/Lc(F);
2671  if (G.isUnivariate() && fdivides(G, F)) return G/Lc(G);
2672  if (F == G) return F/Lc(F);
2673
2674  CFMap M,N;
2675  int best_level= myCompress (A, B, M, N, false);
2676
2677  if (best_level == 0)
2678    return B.genOne();
2679
2680  A= M(A);
2681  B= M(B);
2682
2683  Variable x= Variable (1);
2684  ASSERT (degree (A, x) == 0, "expected degree (F, 1) == 0");
2685  ASSERT (degree (B, x) == 0, "expected degree (G, 1) == 0");
2686
2687  //univariate case
2688  if (A.isUnivariate() && B.isUnivariate())
2689    return N (gcd (A, B));
2690
2691  CanonicalForm skel= M(skeleton);
2692
2693  CanonicalForm cA, cB;    // content of A and B
2694  CanonicalForm ppA, ppB;    // primitive part of A and B
2695  CanonicalForm gcdcAcB;
2696  cA = uni_content (A);
2697  cB = uni_content (B);
2698  gcdcAcB= gcd (cA, cB);
2699  ppA= A/cA;
2700  ppB= B/cB;
2701
2702  CanonicalForm lcA, lcB;  // leading coefficients of A and B
2703  CanonicalForm gcdlcAlcB;
2704  lcA= uni_lcoeff (ppA);
2705  lcB= uni_lcoeff (ppB);
2706
2707  if (fdivides (lcA, lcB))
2708  {
2709    if (fdivides (A, B))
2710      return F/Lc(F);
2711  }
2712  if (fdivides (lcB, lcA))
2713  {
2714    if (fdivides (B, A))
2715      return G/Lc(G);
2716  }
2717
2718  gcdlcAlcB= gcd (lcA, lcB);
2719  int skelSize= size (skel, skel.mvar());
2720
2721  int j= 0;
2722  int biggestSize= 0;
2723  int bufSize;
2724  int numberUni= 0;
2725  for (CFIterator i= skel; i.hasTerms(); i++, j++)
2726  {
2727    bufSize= size (i.coeff());
2728    biggestSize= tmax (biggestSize, bufSize);
2729    numberUni += bufSize;
2730  }
2731  numberUni--;
2732  numberUni= (int) ceil ( (double) numberUni / (skelSize - 1));
2733  biggestSize= tmax (biggestSize , numberUni);
2734
2735  numberUni= biggestSize + size (LC(skel)) - 1;
2736  int biggestSize2= tmax (numberUni, biggestSize);
2737
2738  CanonicalForm g, Aeval, Beval;
2739
2740  CFList evalPoints;
2741  CFArray coeffEval;
2742  bool evalFail= false;
2743  CFList list;
2744  bool GF= false;
2745  CanonicalForm LCA= LC (A);
2746  CanonicalForm tmp;
2747  CFArray gcds= CFArray (biggestSize);
2748  CFList * pEvalPoints= new CFList [biggestSize];
2749  Variable V_buf= alpha;
2750  CFList source, dest;
2751  CanonicalForm prim_elem, im_prim_elem;
2752  for (int i= 0; i < biggestSize; i++)
2753  {
2754    if (i == 0)
2755    {
2756      if (getCharacteristic() > 3)
2757        evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf,
2758                                    evalFail, list);
2759      else
2760        evalFail= true;
2761
2762      if (evalFail)
2763      {
2764        if (V_buf.level() != 1)
2765        {
2766          do
2767          {
2768            Variable V_buf2= chooseExtension (V_buf);
2769            source= CFList();
2770            dest= CFList();
2771
2772            bool prim_fail= false;
2773            Variable V_buf3;
2774            prim_elem= primitiveElement (V_buf, V_buf3, prim_fail);
2775
2776            ASSERT (!prim_fail, "failure in integer factorizer");
2777            if (prim_fail)
2778              ; //ERROR
2779            else
2780              im_prim_elem= mapPrimElem (prim_elem, V_buf, V_buf2);
2781
2782            DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf));
2783            DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf2));
2784
2785            for (CFListIterator i= list; i.hasItem(); i++)
2786              i.getItem()= mapUp (i.getItem(), V_buf, V_buf2, prim_elem,
2787                                im_prim_elem, source, dest);
2788            evalFail= false;
2789            evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf,
2790                                          evalFail, list);
2791          } while (evalFail);
2792        }
2793        else
2794        {
2795          CanonicalForm mipo;
2796          int deg= 2;
2797          do {
2798            mipo= randomIrredpoly (deg, x);
2799            V_buf= rootOf (mipo);
2800            evalFail= false;
2801            evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf,
2802                                          evalFail, list);
2803            deg++;
2804          } while (evalFail);
2805        }
2806      }
2807    }
2808    else
2809    {
2810      mult (evalPoints, pEvalPoints[0]);
2811      eval (A, B, Aeval, Beval, evalPoints);
2812    }
2813
2814    g= gcd (Aeval, Beval);
2815    g /= Lc (g);
2816
2817    if (degree (g) != degree (skel) || (skelSize != size (g)))
2818    {
2819      delete[] pEvalPoints;
2820      fail= true;
2821      return 0;
2822    }
2823    CFIterator l= skel;
2824    for (CFIterator k= g; k.hasTerms(); k++, l++)
2825    {
2826      if (k.exp() != l.exp())
2827      {
2828        delete[] pEvalPoints;
2829        fail= true;
2830        return 0;
2831      }
2832    }
2833    pEvalPoints[i]= evalPoints;
2834    gcds[i]= g;
2835
2836    tmp= 0;
2837    int j= 0;
2838    for (CFListIterator k= evalPoints; k.hasItem(); k++, j++)
2839      tmp += k.getItem()*power (x, j);
2840    list.append (tmp);
2841  }
2842
2843  if (Monoms.size() == 0)
2844    Monoms= getMonoms (skel);
2845
2846  if (coeffMonoms == NULL)
2847    coeffMonoms= new CFArray [skelSize];
2848
2849  j= 0;
2850  for (CFIterator i= skel; i.hasTerms(); i++, j++)
2851    coeffMonoms[j]= getMonoms (i.coeff());
2852
2853  int minimalColumnsIndex;
2854  if (skelSize > 1)
2855    minimalColumnsIndex= 1;
2856  else
2857    minimalColumnsIndex= 0;
2858  int minimalColumns=-1;
2859
2860  CFArray* pM= new CFArray [skelSize];
2861  CFMatrix Mat;
2862  // find the Matrix with minimal number of columns
2863  for (int i= 0; i < skelSize; i++)
2864  {
2865    pM[i]= CFArray (coeffMonoms[i].size());
2866    if (i == 1)
2867      minimalColumns= coeffMonoms[i].size();
2868    if (i > 1)
2869    {
2870      if (minimalColumns > coeffMonoms[i].size())
2871      {
2872        minimalColumns= coeffMonoms[i].size();
2873        minimalColumnsIndex= i;
2874      }
2875    }
2876  }
2877  CFMatrix* pMat= new CFMatrix [2];
2878  pMat[0]= CFMatrix (biggestSize, coeffMonoms[0].size());
2879  pMat[1]= CFMatrix (biggestSize, coeffMonoms[minimalColumnsIndex].size());
2880  CFArray* pL= new CFArray [skelSize];
2881  for (int i= 0; i < biggestSize; i++)
2882  {
2883    CFIterator l= gcds [i];
2884    evalPoints= pEvalPoints [i];
2885    for (int k= 0; k < skelSize; k++, l++)
2886    {
2887      if (i == 0)
2888        pL[k]= CFArray (biggestSize);
2889      pL[k] [i]= l.coeff();
2890
2891      if (i == 0 && k != 0 && k != minimalColumnsIndex)
2892        pM[k]= evaluate (coeffMonoms[k], evalPoints);
2893      else if (i == 0 && (k == 0 || k == minimalColumnsIndex))
2894        coeffEval= evaluate (coeffMonoms[k], evalPoints);
2895      if (i > 0 && (k == 0 || k == minimalColumnsIndex))
2896        coeffEval= evaluate (coeffMonoms[k], evalPoints);
2897
2898      if (k == 0)
2899      {
2900        if (pMat[k].rows() >= i + 1)
2901        {
2902          for (int ind= 1; ind < coeffEval.size() + 1; ind++)
2903            pMat[k] (i + 1, ind)= coeffEval[ind - 1];
2904        }
2905      }
2906      if (k == minimalColumnsIndex)
2907      {
2908        if (pMat[1].rows() >= i + 1)
2909        {
2910          for (int ind= 1; ind < coeffEval.size() + 1; ind++)
2911            pMat[1] (i + 1, ind)= coeffEval[ind - 1];
2912        }
2913      }
2914    }
2915  }
2916
2917  CFArray solution;
2918  CanonicalForm result= 0;
2919  int ind= 1;
2920  int matRows, matColumns;
2921  matRows= pMat[1].rows();
2922  matColumns= pMat[0].columns() - 1;
2923  matColumns += pMat[1].columns();
2924
2925  Mat= CFMatrix (matRows, matColumns);
2926  for (int i= 1; i <= matRows; i++)
2927    for (int j= 1; j <= pMat[1].columns(); j++)
2928      Mat (i, j)= pMat[1] (i, j);
2929
2930  for (int j= pMat[1].columns() + 1; j <= matColumns;
2931       j++, ind++)
2932  {
2933    for (int i= 1; i <= matRows; i++)
2934      Mat (i, j)= (-pMat [0] (i, ind + 1))*pL[minimalColumnsIndex] [i - 1];
2935  }
2936
2937  CFArray firstColumn= CFArray (pMat[0].rows());
2938  for (int i= 0; i < pMat[0].rows(); i++)
2939    firstColumn [i]= pMat[0] (i + 1, 1);
2940
2941  CFArray bufArray= pL[minimalColumnsIndex];
2942
2943  for (int i= 0; i < biggestSize; i++)
2944    pL[minimalColumnsIndex] [i] *= firstColumn[i];
2945
2946  CFMatrix bufMat= pMat[1];
2947  pMat[1]= Mat;
2948
2949  if (V_buf.level() != 1)
2950    solution= solveSystemFq (pMat[1],
2951                             pL[minimalColumnsIndex], V_buf);
2952  else
2953    solution= solveSystemFp (pMat[1],
2954                             pL[minimalColumnsIndex]);
2955
2956  if (solution.size() == 0)
2957  { //Javadi's method failed, try deKleine, Monagan, Wittkopf instead
2958    CFMatrix bufMat0= pMat[0];
2959    delete [] pMat;
2960    pMat= new CFMatrix [skelSize];
2961    pL[minimalColumnsIndex]= bufArray;
2962    CFList* bufpEvalPoints= NULL;
2963    CFArray bufGcds;
2964    if (biggestSize != biggestSize2)
2965    {
2966      bufpEvalPoints= pEvalPoints;
2967      pEvalPoints= new CFList [biggestSize2];
2968      bufGcds= gcds;
2969      gcds= CFArray (biggestSize2);
2970      for (int i= 0; i < biggestSize; i++)
2971      {
2972        pEvalPoints[i]= bufpEvalPoints [i];
2973        gcds[i]= bufGcds[i];
2974      }
2975      for (int i= 0; i < biggestSize2 - biggestSize; i++)
2976      {
2977        mult (evalPoints, pEvalPoints[0]);
2978        eval (A, B, Aeval, Beval, evalPoints);
2979        g= gcd (Aeval, Beval);
2980        g /= Lc (g);
2981        if (degree (g) != degree (skel) || (skelSize != size (g)))
2982        {
2983          delete[] pEvalPoints;
2984          delete[] pMat;
2985          delete[] pL;
2986          delete[] coeffMonoms;
2987          delete[] pM;
2988          if (bufpEvalPoints != NULL)
2989            delete [] bufpEvalPoints;
2990          fail= true;
2991          return 0;
2992        }
2993        CFIterator l= skel;
2994        for (CFIterator k= g; k.hasTerms(); k++, l++)
2995        {
2996          if (k.exp() != l.exp())
2997          {
2998            delete[] pEvalPoints;
2999            delete[] pMat;
3000            delete[] pL;
3001            delete[] coeffMonoms;
3002            delete[] pM;
3003            if (bufpEvalPoints != NULL)
3004              delete [] bufpEvalPoints;
3005            fail= true;
3006            return 0;
3007          }
3008        }
3009        pEvalPoints[i + biggestSize]= evalPoints;
3010        gcds[i + biggestSize]= g;
3011      }
3012    }
3013    for (int i= 0; i < biggestSize; i++)
3014    {
3015      CFIterator l= gcds [i];
3016      evalPoints= pEvalPoints [i];
3017      for (int k= 1; k < skelSize; k++, l++)
3018      {
3019        if (i == 0)
3020          pMat[k]= CFMatrix (biggestSize2,coeffMonoms[k].size()+biggestSize2-1);
3021        if (k == minimalColumnsIndex)
3022          continue;
3023        coeffEval= evaluate (coeffMonoms[k], evalPoints);
3024        if (pMat[k].rows() >= i + 1)
3025        {
3026          for (int ind= 1; ind < coeffEval.size() + 1; ind++)
3027            pMat[k] (i + 1, ind)= coeffEval[ind - 1];
3028        }
3029      }
3030    }
3031    Mat= bufMat0;
3032    pMat[0]= CFMatrix (biggestSize2, coeffMonoms[0].size() + biggestSize2 - 1);
3033    for (int j= 1; j <= Mat.rows(); j++)
3034      for (int k= 1; k <= Mat.columns(); k++)
3035        pMat [0] (j,k)= Mat (j,k);
3036    Mat= bufMat;
3037    for (int j= 1; j <= Mat.rows(); j++)
3038      for (int k= 1; k <= Mat.columns(); k++)
3039        pMat [minimalColumnsIndex] (j,k)= Mat (j,k);
3040    // write old matrix entries into new matrices
3041    for (int i= 0; i < skelSize; i++)
3042    {
3043      bufArray= pL[i];
3044      pL[i]= CFArray (biggestSize2);
3045      for (int j= 0; j < bufArray.size(); j++)
3046        pL[i] [j]= bufArray [j];
3047    }
3048    //write old vector entries into new and add new entries to old matrices
3049    for (int i= 0; i < biggestSize2 - biggestSize; i++)
3050    {
3051      CFIterator l= gcds [i + biggestSize];
3052      evalPoints= pEvalPoints [i + biggestSize];
3053      for (int k= 0; k < skelSize; k++, l++)
3054      {
3055        pL[k] [i + biggestSize]= l.coeff();
3056        coeffEval= evaluate (coeffMonoms[k], evalPoints);
3057        if (pMat[k].rows() >= i + biggestSize + 1)
3058        {
3059          for (int ind= 1; ind < coeffEval.size() + 1; ind++)
3060            pMat[k] (i + biggestSize + 1, ind)= coeffEval[ind - 1];
3061        }
3062      }
3063    }
3064    // begin new
3065    for (int i= 0; i < skelSize; i++)
3066    {
3067      if (pL[i].size() > 1)
3068      {
3069        for (int j= 2; j <= pMat[i].rows(); j++)
3070          pMat[i] (j, coeffMonoms[i].size() + j - 1)=
3071              -pL[i] [j - 1];
3072      }
3073    }
3074
3075    matColumns= biggestSize2 - 1;
3076    matRows= 0;
3077    for (int i= 0; i < skelSize; i++)
3078    {
3079      if (V_buf.level() == 1)
3080        (void) gaussianElimFp (pMat[i], pL[i]);
3081      else
3082        (void) gaussianElimFq (pMat[i], pL[i], V_buf);
3083
3084      if (pMat[i] (coeffMonoms[i].size(), coeffMonoms[i].size()) == 0)
3085      {
3086        delete[] pEvalPoints;
3087        delete[] pMat;
3088        delete[] pL;
3089        delete[] coeffMonoms;
3090        delete[] pM;
3091        if (bufpEvalPoints != NULL)
3092          delete [] bufpEvalPoints;
3093        fail= true;
3094        return 0;
3095      }
3096      matRows += pMat[i].rows() - coeffMonoms[i].size();
3097    }
3098
3099    CFMatrix bufMat;
3100    Mat= CFMatrix (matRows, matColumns);
3101    ind= 0;
3102    bufArray= CFArray (matRows);
3103    CFArray bufArray2;
3104    for (int i= 0; i < skelSize; i++)
3105    {
3106      bufMat= pMat[i] (coeffMonoms[i].size() + 1, pMat[i].rows(),
3107                       coeffMonoms[i].size() + 1, pMat[i].columns());
3108
3109      for (int j= 1; j <= bufMat.rows(); j++)
3110        for (int k= 1; k <= bufMat.columns(); k++)
3111          Mat (j + ind, k)= bufMat(j, k);
3112      bufArray2= coeffMonoms[i].size();
3113      for (int j= 1; j <= pMat[i].rows(); j++)
3114      {
3115        if (j > coeffMonoms[i].size())
3116          bufArray [j-coeffMonoms[i].size() + ind - 1]= pL[i] [j - 1];
3117        else
3118          bufArray2 [j - 1]= pL[i] [j - 1];
3119      }
3120      pL[i]= bufArray2;
3121      ind += bufMat.rows();
3122      pMat[i]= pMat[i] (1, coeffMonoms[i].size(), 1, pMat[i].columns());
3123    }
3124
3125    if (V_buf.level() != 1)
3126      solution= solveSystemFq (Mat, bufArray, V_buf);
3127    else
3128      solution= solveSystemFp (Mat, bufArray);
3129
3130    if (solution.size() == 0)
3131    {
3132      delete[] pEvalPoints;
3133      delete[] pMat;
3134      delete[] pL;
3135      delete[] coeffMonoms;
3136      delete[] pM;
3137      if (bufpEvalPoints != NULL)
3138        delete [] bufpEvalPoints;
3139      fail= true;
3140      return 0;
3141    }
3142
3143    ind= 0;
3144    result= 0;
3145    CFArray bufSolution;
3146    for (int i= 0; i < skelSize; i++)
3147    {
3148      bufSolution= readOffSolution (pMat[i], pL[i], solution);
3149      for (int i= 0; i < bufSolution.size(); i++, ind++)
3150        result += Monoms [ind]*bufSolution[i];
3151    }
3152    if (alpha.level() != 1 && V_buf != alpha)
3153    {
3154      CFList u, v;
3155      result= mapDown (result, prim_elem, im_prim_elem, alpha, u, v);
3156    }
3157    result= N(result);
3158    if (fdivides (result, F) && fdivides (result, G))
3159    {
3160      delete[] pEvalPoints;
3161      delete[] pMat;
3162      delete[] pL;
3163      delete[] pM;
3164      if (bufpEvalPoints != NULL)
3165        delete [] bufpEvalPoints;
3166      return result;
3167    }
3168    else
3169    {
3170      delete[] pEvalPoints;
3171      delete[] pMat;
3172      delete[] pL;
3173      delete[] coeffMonoms;
3174      delete[] pM;
3175      if (bufpEvalPoints != NULL)
3176        delete [] bufpEvalPoints;
3177      fail= true;
3178      return 0;
3179    }
3180  } // end of deKleine, Monagan & Wittkopf
3181
3182  result += Monoms[0];
3183  int ind2= 0, ind3= 2;
3184  ind= 0;
3185  for (int l= 0; l < minimalColumnsIndex; l++)
3186    ind += coeffMonoms[l].size();
3187  for (int l= solution.size() - pMat[0].columns() + 1; l < solution.size();
3188       l++, ind2++, ind3++)
3189  {
3190    result += solution[l]*Monoms [1 + ind2];
3191    for (int i= 0; i < pMat[0].rows(); i++)
3192      firstColumn[i] += solution[l]*pMat[0] (i + 1, ind3);
3193  }
3194  for (int l= 0; l < solution.size() + 1 - pMat[0].columns(); l++)
3195    result += solution[l]*Monoms [ind + l];
3196  ind= coeffMonoms[0].size();
3197  for (int k= 1; k < skelSize; k++)
3198  {
3199    if (k == minimalColumnsIndex)
3200    {
3201      ind += coeffMonoms[k].size();
3202      continue;
3203    }
3204    if (k != minimalColumnsIndex)
3205    {
3206      for (int i= 0; i < biggestSize; i++)
3207        pL[k] [i] *= firstColumn [i];
3208    }
3209
3210    if (biggestSize != coeffMonoms[k].size() && k != minimalColumnsIndex)
3211    {
3212      bufArray= CFArray (coeffMonoms[k].size());
3213      for (int i= 0; i < bufArray.size(); i++)
3214        bufArray [i]= pL[k] [i];
3215      solution= solveGeneralVandermonde (pM[k], bufArray);
3216    }
3217    else
3218      solution= solveGeneralVandermonde (pM[k], pL[k]);
3219
3220    if (solution.size() == 0)
3221    {
3222      delete[] pEvalPoints;
3223      delete[] pMat;
3224      delete[] pL;
3225      delete[] coeffMonoms;
3226      delete[] pM;
3227      fail= true;
3228      return 0;
3229    }
3230    if (k != minimalColumnsIndex)
3231    {
3232      for (int l= 0; l < solution.size(); l++)
3233        result += solution[l]*Monoms [ind + l];
3234      ind += solution.size();
3235    }
3236  }
3237
3238  delete[] pEvalPoints;
3239  delete[] pMat;
3240  delete[] pL;
3241  delete[] pM;
3242
3243  if (alpha.level() != 1 && V_buf != alpha)
3244  {
3245    CFList u, v;
3246    result= mapDown (result, prim_elem, im_prim_elem, alpha, u, v);
3247  }
3248  result= N(result);
3249
3250  if (fdivides (result, F) && fdivides (result, G))
3251    return result;
3252  else
3253  {
3254    delete[] coeffMonoms;
3255    fail= true;
3256    return 0;
3257  }
3258}
3259
3260CanonicalForm sparseGCDFq (const CanonicalForm& F, const CanonicalForm& G,
3261                           const Variable & alpha, CFList& l, bool& topLevel)
3262{
3263  CanonicalForm A= F;
3264  CanonicalForm B= G;
3265  if (F.isZero() && degree(G) > 0) return G/Lc(G);
3266  else if (G.isZero() && degree (F) > 0) return F/Lc(F);
3267  else if (F.isZero() && G.isZero()) return F.genOne();
3268  if (F.inBaseDomain() || G.inBaseDomain()) return F.genOne();
3269  if (F.isUnivariate() && fdivides(F, G)) return F/Lc(F);
3270  if (G.isUnivariate() && fdivides(G, F)) return G/Lc(G);
3271  if (F == G) return F/Lc(F);
3272
3273  CFMap M,N;
3274  int best_level= myCompress (A, B, M, N, topLevel);
3275
3276  if (best_level == 0) return B.genOne();
3277
3278  A= M(A);
3279  B= M(B);
3280
3281  Variable x= Variable (1);
3282
3283  //univariate case
3284  if (A.isUnivariate() && B.isUnivariate())
3285    return N (gcd (A, B));
3286
3287  CanonicalForm cA, cB;    // content of A and B
3288  CanonicalForm ppA, ppB;    // primitive part of A and B
3289  CanonicalForm gcdcAcB;
3290
3291  cA = uni_content (A);
3292  cB = uni_content (B);
3293  gcdcAcB= gcd (cA, cB);
3294  ppA= A/cA;
3295  ppB= B/cB;
3296
3297  CanonicalForm lcA, lcB;  // leading coefficients of A and B
3298  CanonicalForm gcdlcAlcB;
3299  lcA= uni_lcoeff (ppA);
3300  lcB= uni_lcoeff (ppB);
3301
3302  if (fdivides (lcA, lcB))
3303  {
3304    if (fdivides (A, B))
3305      return F/Lc(F);
3306  }
3307  if (fdivides (lcB, lcA))
3308  {
3309    if (fdivides (B, A))
3310      return G/Lc(G);
3311  }
3312
3313  gcdlcAlcB= gcd (lcA, lcB);
3314
3315  int d= totaldegree (ppA, Variable (2), Variable (ppA.level()));
3316  int d0;
3317
3318  if (d == 0)
3319    return N(gcdcAcB);
3320  d0= totaldegree (ppB, Variable (2), Variable (ppB.level()));
3321
3322  if (d0 < d)
3323    d= d0;
3324
3325  if (d == 0)
3326    return N(gcdcAcB);
3327
3328  CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH, skeleton;
3329  CanonicalForm newtonPoly= 1;
3330  m= gcdlcAlcB;
3331  G_m= 0;
3332  H= 0;
3333  bool fail= false;
3334  topLevel= false;
3335  bool inextension= false;
3336  Variable V_buf= alpha;
3337  CanonicalForm prim_elem, im_prim_elem;
3338  CFList source, dest;
3339  do // first do
3340  {
3341    random_element= randomElement (m, V_buf, l, fail);
3342    if (random_element == 0 && !fail)
3343    {
3344      if (!find (l, random_element))
3345        l.append (random_element);
3346      continue;
3347    }
3348    if (fail)
3349    {
3350      source= CFList();
3351      dest= CFList();
3352
3353      Variable V_buf3= V_buf;
3354      V_buf= chooseExtension (V_buf);
3355      bool prim_fail= false;
3356      Variable V_buf2;
3357      prim_elem= primitiveElement (alpha, V_buf2, prim_fail);
3358
3359      if (V_buf3 != alpha)
3360      {
3361        m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
3362        G_m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
3363        newtonPoly= mapDown (newtonPoly, prim_elem, im_prim_elem, alpha,
3364                             source, dest);
3365        ppA= mapDown (ppA, prim_elem, im_prim_elem, alpha, source, dest);
3366        ppB= mapDown (ppB, prim_elem, im_prim_elem, alpha, source, dest);
3367        gcdlcAlcB= mapDown (gcdlcAlcB, prim_elem, im_prim_elem, alpha, source,
3368                            dest);
3369        for (CFListIterator i= l; i.hasItem(); i++)
3370          i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, alpha,
3371                                source, dest);
3372      }
3373
3374      ASSERT (!prim_fail, "failure in integer factorizer");
3375      if (prim_fail)
3376        ; //ERROR
3377      else
3378        im_prim_elem= mapPrimElem (prim_elem, alpha, V_buf);
3379
3380      DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (alpha));
3381      DEBOUTLN (cerr, "getMipo (V_buf2)= " << getMipo (V_buf2));
3382      inextension= true;
3383      for (CFListIterator i= l; i.hasItem(); i++)
3384        i.getItem()= mapUp (i.getItem(), alpha, V_buf, prim_elem,
3385                             im_prim_elem, source, dest);
3386      m= mapUp (m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3387      G_m= mapUp (G_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3388      newtonPoly= mapUp (newtonPoly, alpha, V_buf, prim_elem, im_prim_elem,
3389                          source, dest);
3390      ppA= mapUp (ppA, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3391      ppB= mapUp (ppB, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3392      gcdlcAlcB= mapUp (gcdlcAlcB, alpha, V_buf, prim_elem, im_prim_elem,
3393                         source, dest);
3394
3395      fail= false;
3396      random_element= randomElement (m, V_buf, l, fail );
3397      DEBOUTLN (cerr, "fail= " << fail);
3398      CFList list;
3399      TIMING_START (gcd_recursion);
3400      G_random_element=
3401      sparseGCDFq (ppA (random_element, x), ppB (random_element, x), V_buf,
3402                        list, topLevel);
3403      TIMING_END_AND_PRINT (gcd_recursion,
3404                            "time for recursive call: ");
3405      DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3406    }
3407    else
3408    {
3409      CFList list;
3410      TIMING_START (gcd_recursion);
3411      G_random_element=
3412      sparseGCDFq (ppA(random_element, x), ppB(random_element, x), V_buf,
3413                        list, topLevel);
3414      TIMING_END_AND_PRINT (gcd_recursion,
3415                            "time for recursive call: ");
3416      DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3417    }
3418
3419    if (!G_random_element.inCoeffDomain())
3420      d0= totaldegree (G_random_element, Variable(2),
3421                       Variable (G_random_element.level()));
3422    else
3423      d0= 0;
3424
3425    if (d0 == 0)
3426      return N(gcdcAcB);
3427    if (d0 >  d)
3428    {
3429      if (!find (l, random_element))
3430        l.append (random_element);
3431      continue;
3432    }
3433
3434    G_random_element=
3435    (gcdlcAlcB(random_element, x)/uni_lcoeff (G_random_element))
3436    * G_random_element;
3437
3438    skeleton= G_random_element;
3439    if (!G_random_element.inCoeffDomain())
3440      d0= totaldegree (G_random_element, Variable(2),
3441                       Variable (G_random_element.level()));
3442    else
3443      d0= 0;
3444
3445    if (d0 <  d)
3446    {
3447      m= gcdlcAlcB;
3448      newtonPoly= 1;
3449      G_m= 0;
3450      d= d0;
3451    }
3452
3453    H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
3454    if (uni_lcoeff (H) == gcdlcAlcB)
3455    {
3456      cH= uni_content (H);
3457      ppH= H/cH;
3458      if (inextension)
3459      {
3460        CFList u, v;
3461        //maybe it's better to test if ppH is an element of F(\alpha) before
3462        //mapping down
3463        if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
3464        {
3465          DEBOUTLN (cerr, "ppH before mapDown= " << ppH);
3466          ppH= mapDown (ppH, prim_elem, im_prim_elem, alpha, u, v);
3467          ppH /= Lc(ppH);
3468          DEBOUTLN (cerr, "ppH after mapDown= " << ppH);
3469          return N(gcdcAcB*ppH);
3470        }
3471      }
3472      else if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
3473        return N(gcdcAcB*ppH);
3474    }
3475    G_m= H;
3476    newtonPoly= newtonPoly*(x - random_element);
3477    m= m*(x - random_element);
3478    if (!find (l, random_element))
3479      l.append (random_element);
3480
3481    if (getCharacteristic () > 3 && size (skeleton) < 100)
3482    {
3483      CFArray Monoms;
3484      CFArray *coeffMonoms= NULL;
3485      do //second do
3486      {
3487        random_element= randomElement (m, V_buf, l, fail);
3488        if (random_element == 0 && !fail)
3489        {
3490          if (!find (l, random_element))
3491            l.append (random_element);
3492          continue;
3493        }
3494        if (fail)
3495        {
3496          source= CFList();
3497          dest= CFList();
3498
3499          Variable V_buf3= V_buf;
3500          V_buf= chooseExtension (V_buf);
3501          bool prim_fail= false;
3502          Variable V_buf2;
3503          prim_elem= primitiveElement (alpha, V_buf2, prim_fail);
3504
3505          if (V_buf3 != alpha)
3506          {
3507            m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
3508            G_m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
3509            newtonPoly= mapDown (newtonPoly, prim_elem, im_prim_elem, alpha,
3510                                 source, dest);
3511            ppA= mapDown (ppA, prim_elem, im_prim_elem, alpha, source, dest);
3512            ppB= mapDown (ppB, prim_elem, im_prim_elem, alpha, source, dest);
3513            gcdlcAlcB= mapDown (gcdlcAlcB, prim_elem, im_prim_elem, alpha,
3514                                source, dest);
3515            for (CFListIterator i= l; i.hasItem(); i++)
3516              i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, alpha,
3517                                    source, dest);
3518          }
3519
3520          ASSERT (!prim_fail, "failure in integer factorizer");
3521          if (prim_fail)
3522            ; //ERROR
3523          else
3524            im_prim_elem= mapPrimElem (prim_elem, alpha, V_buf);
3525
3526          DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (alpha));
3527          DEBOUTLN (cerr, "getMipo (V_buf2)= " << getMipo (V_buf2));
3528          inextension= true;
3529          for (CFListIterator i= l; i.hasItem(); i++)
3530            i.getItem()= mapUp (i.getItem(), alpha, V_buf, prim_elem,
3531                                im_prim_elem, source, dest);
3532          m= mapUp (m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3533          G_m= mapUp (G_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3534          newtonPoly= mapUp (newtonPoly, alpha, V_buf, prim_elem, im_prim_elem,
3535                              source, dest);
3536          ppA= mapUp (ppA, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3537          ppB= mapUp (ppB, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3538
3539          gcdlcAlcB= mapUp (gcdlcAlcB, alpha, V_buf, prim_elem, im_prim_elem,
3540                            source, dest);
3541
3542          fail= false;
3543          random_element= randomElement (m, V_buf, l, fail);
3544          DEBOUTLN (cerr, "fail= " << fail);
3545          CFList list;
3546          TIMING_START (gcd_recursion);
3547
3548          //sparseInterpolation
3549          bool sparseFail= false;
3550          if (LC (skeleton).inCoeffDomain())
3551            G_random_element=
3552            monicSparseInterpol (ppA (random_element, x), ppB(random_element,x),
3553                                skeleton,V_buf, sparseFail, coeffMonoms,Monoms);
3554          else
3555            G_random_element=
3556            nonMonicSparseInterpol (ppA(random_element,x),ppB(random_element,x),
3557                                    skeleton, V_buf, sparseFail, coeffMonoms,
3558                                    Monoms);
3559          if (sparseFail)
3560            break;
3561
3562          TIMING_END_AND_PRINT (gcd_recursion,
3563                                "time for recursive call: ");
3564          DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3565        }
3566        else
3567        {
3568          CFList list;
3569          TIMING_START (gcd_recursion);
3570          bool sparseFail= false;
3571          if (LC (skeleton).inCoeffDomain())
3572            G_random_element=
3573            monicSparseInterpol (ppA (random_element, x),ppB(random_element, x),
3574                                skeleton,V_buf, sparseFail,coeffMonoms, Monoms);
3575          else
3576            G_random_element=
3577            nonMonicSparseInterpol (ppA(random_element,x),ppB(random_element,x),
3578                                    skeleton, V_buf, sparseFail, coeffMonoms,
3579                                    Monoms);
3580          if (sparseFail)
3581            break;
3582
3583          TIMING_END_AND_PRINT (gcd_recursion,
3584                                "time for recursive call: ");
3585          DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3586        }
3587
3588        if (!G_random_element.inCoeffDomain())
3589          d0= totaldegree (G_random_element, Variable(2),
3590                           Variable (G_random_element.level()));
3591        else
3592          d0= 0;
3593
3594        if (d0 == 0)
3595          return N(gcdcAcB);
3596        if (d0 >  d)
3597        {
3598          if (!find (l, random_element))
3599            l.append (random_element);
3600          continue;
3601        }
3602
3603        G_random_element=
3604        (gcdlcAlcB(random_element, x)/uni_lcoeff (G_random_element))
3605        * G_random_element;
3606
3607        if (!G_random_element.inCoeffDomain())
3608          d0= totaldegree (G_random_element, Variable(2),
3609                          Variable (G_random_element.level()));
3610        else
3611          d0= 0;
3612
3613        if (d0 <  d)
3614        {
3615          m= gcdlcAlcB;
3616          newtonPoly= 1;
3617          G_m= 0;
3618          d= d0;
3619        }
3620
3621        TIMING_START (newton_interpolation);
3622        H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
3623        TIMING_END_AND_PRINT (newton_interpolation,
3624                              "time for newton interpolation: ");
3625
3626        //termination test
3627        if (uni_lcoeff (H) == gcdlcAlcB)
3628        {
3629          cH= uni_content (H);
3630          ppH= H/cH;
3631          if (inextension)
3632          {
3633            CFList u, v;
3634            //maybe it's better to test if ppH is an element of F(\alpha) before
3635            //mapping down
3636            if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
3637            {
3638              DEBOUTLN (cerr, "ppH before mapDown= " << ppH);
3639              ppH= mapDown (ppH, prim_elem, im_prim_elem, alpha, u, v);
3640              ppH /= Lc(ppH);
3641              DEBOUTLN (cerr, "ppH after mapDown= " << ppH);
3642              return N(gcdcAcB*ppH);
3643            }
3644          }
3645          else if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
3646          {
3647            return N(gcdcAcB*ppH);
3648          }
3649        }
3650
3651        G_m= H;
3652        newtonPoly= newtonPoly*(x - random_element);
3653        m= m*(x - random_element);
3654        if (!find (l, random_element))
3655          l.append (random_element);
3656
3657      } while (1);
3658    }
3659  } while (1);
3660}
3661
3662CanonicalForm sparseGCDFp (const CanonicalForm& F, const CanonicalForm& G,
3663                           bool& topLevel, CFList& l)
3664{
3665  CanonicalForm A= F;
3666  CanonicalForm B= G;
3667  if (F.isZero() && degree(G) > 0) return G/Lc(G);
3668  else if (G.isZero() && degree (F) > 0) return F/Lc(F);
3669  else if (F.isZero() && G.isZero()) return F.genOne();
3670  if (F.inBaseDomain() || G.inBaseDomain()) return F.genOne();
3671  if (F.isUnivariate() && fdivides(F, G)) return F/Lc(F);
3672  if (G.isUnivariate() && fdivides(G, F)) return G/Lc(G);
3673  if (F == G) return F/Lc(F);
3674
3675  CFMap M,N;
3676  int best_level= myCompress (A, B, M, N, topLevel);
3677
3678  if (best_level == 0) return B.genOne();
3679
3680  A= M(A);
3681  B= M(B);
3682
3683  Variable x= Variable (1);
3684
3685  //univariate case
3686  if (A.isUnivariate() && B.isUnivariate())
3687    return N (gcd (A, B));
3688
3689  CanonicalForm cA, cB;    // content of A and B
3690  CanonicalForm ppA, ppB;    // primitive part of A and B
3691  CanonicalForm gcdcAcB;
3692
3693  cA = uni_content (A);
3694  cB = uni_content (B);
3695  gcdcAcB= gcd (cA, cB);
3696  ppA= A/cA;
3697  ppB= B/cB;
3698
3699  CanonicalForm lcA, lcB;  // leading coefficients of A and B
3700  CanonicalForm gcdlcAlcB;
3701  lcA= uni_lcoeff (ppA);
3702  lcB= uni_lcoeff (ppB);
3703
3704  if (fdivides (lcA, lcB))
3705  {
3706    if (fdivides (A, B))
3707      return F/Lc(F);
3708  }
3709  if (fdivides (lcB, lcA))
3710  {
3711    if (fdivides (B, A))
3712      return G/Lc(G);
3713  }
3714
3715  gcdlcAlcB= gcd (lcA, lcB);
3716
3717  int d= totaldegree (ppA, Variable (2), Variable (ppA.level()));
3718  int d0;
3719
3720  if (d == 0)
3721    return N(gcdcAcB);
3722
3723  d0= totaldegree (ppB, Variable (2), Variable (ppB.level()));
3724
3725  if (d0 < d)
3726    d= d0;
3727
3728  if (d == 0)
3729    return N(gcdcAcB);
3730
3731  CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH, skeleton;
3732  CanonicalForm newtonPoly= 1;
3733  m= gcdlcAlcB;
3734  G_m= 0;
3735  H= 0;
3736  bool fail= false;
3737  topLevel= false;
3738  bool inextension= false;
3739  Variable V_buf, alpha;
3740  CanonicalForm prim_elem, im_prim_elem;
3741  CFList source, dest;
3742  do //first do
3743  {
3744    if (inextension)
3745      random_element= randomElement (m, V_buf, l, fail);
3746    else
3747      random_element= FpRandomElement (m, l, fail);
3748    if (random_element == 0 && !fail)
3749    {
3750      if (!find (l, random_element))
3751        l.append (random_element);
3752      continue;
3753    }
3754
3755    if (!fail && !inextension)
3756    {
3757      CFList list;
3758      TIMING_START (gcd_recursion);
3759      G_random_element=
3760      sparseGCDFp (ppA (random_element,x), ppB (random_element,x), topLevel,
3761                   list);
3762      TIMING_END_AND_PRINT (gcd_recursion,
3763                            "time for recursive call: ");
3764      DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3765    }
3766    else if (!fail && inextension)
3767    {
3768      CFList list;
3769      TIMING_START (gcd_recursion);
3770      G_random_element=
3771      sparseGCDFq (ppA (random_element, x), ppB (random_element, x), alpha,
3772                        list, topLevel);
3773      TIMING_END_AND_PRINT (gcd_recursion,
3774                            "time for recursive call: ");
3775      DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3776    }
3777    else if (fail && !inextension)
3778    {
3779      source= CFList();
3780      dest= CFList();
3781      CFList list;
3782      CanonicalForm mipo;
3783      int deg= 2;
3784      do
3785      {
3786        mipo= randomIrredpoly (deg, x);
3787        alpha= rootOf (mipo);
3788        inextension= true;
3789        fail= false;
3790        random_element= randomElement (m, alpha, l, fail);
3791        deg++;
3792      } while (fail);
3793      V_buf= alpha;
3794      list= CFList();
3795      TIMING_START (gcd_recursion);
3796      G_random_element=
3797      sparseGCDFq (ppA (random_element, x), ppB (random_element, x), alpha,
3798                        list, topLevel);
3799      TIMING_END_AND_PRINT (gcd_recursion,
3800                            "time for recursive call: ");
3801      DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3802    }
3803    else if (fail && inextension)
3804    {
3805      source= CFList();
3806      dest= CFList();
3807
3808      Variable V_buf3= V_buf;
3809      V_buf= chooseExtension (V_buf);
3810      bool prim_fail= false;
3811      Variable V_buf2;
3812      CanonicalForm prim_elem, im_prim_elem;
3813      prim_elem= primitiveElement (alpha, V_buf2, prim_fail);
3814
3815      if (V_buf3 != alpha)
3816      {
3817        m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
3818        G_m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
3819        newtonPoly= mapDown (newtonPoly, prim_elem, im_prim_elem, alpha, source,
3820                             dest);
3821        ppA= mapDown (ppA, prim_elem, im_prim_elem, alpha, source, dest);
3822        ppB= mapDown (ppB, prim_elem, im_prim_elem, alpha, source, dest);
3823        gcdlcAlcB= mapDown (gcdlcAlcB, prim_elem, im_prim_elem, alpha, source,
3824                            dest);
3825        for (CFListIterator i= l; i.hasItem(); i++)
3826          i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, alpha,
3827                                source, dest);
3828      }
3829
3830      ASSERT (!prim_fail, "failure in integer factorizer");
3831      if (prim_fail)
3832        ; //ERROR
3833      else
3834        im_prim_elem= mapPrimElem (prim_elem, alpha, V_buf);
3835
3836      DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (alpha));
3837      DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf2));
3838
3839      for (CFListIterator i= l; i.hasItem(); i++)
3840        i.getItem()= mapUp (i.getItem(), alpha, V_buf, prim_elem,
3841                             im_prim_elem, source, dest);
3842      m= mapUp (m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3843      G_m= mapUp (G_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3844      newtonPoly= mapUp (newtonPoly, alpha, V_buf, prim_elem, im_prim_elem,
3845                          source, dest);
3846      ppA= mapUp (ppA, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3847      ppB= mapUp (ppB, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3848      gcdlcAlcB= mapUp (gcdlcAlcB, alpha, V_buf, prim_elem, im_prim_elem,
3849                         source, dest);
3850      fail= false;
3851      random_element= randomElement (m, V_buf, l, fail );
3852      DEBOUTLN (cerr, "fail= " << fail);
3853      CFList list;
3854      TIMING_START (gcd_recursion);
3855      G_random_element=
3856      sparseGCDFq (ppA (random_element, x), ppB (random_element, x), V_buf,
3857                        list, topLevel);
3858      TIMING_END_AND_PRINT (gcd_recursion,
3859                            "time for recursive call: ");
3860      DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3861    }
3862
3863    if (!G_random_element.inCoeffDomain())
3864      d0= totaldegree (G_random_element, Variable(2),
3865                       Variable (G_random_element.level()));
3866    else
3867      d0= 0;
3868
3869    if (d0 == 0)
3870      return N(gcdcAcB);
3871    if (d0 >  d)
3872    {
3873      if (!find (l, random_element))
3874        l.append (random_element);
3875      continue;
3876    }
3877
3878    G_random_element=
3879    (gcdlcAlcB(random_element, x)/uni_lcoeff (G_random_element))
3880    * G_random_element;
3881
3882    skeleton= G_random_element;
3883
3884    if (!G_random_element.inCoeffDomain())
3885      d0= totaldegree (G_random_element, Variable(2),
3886                       Variable (G_random_element.level()));
3887    else
3888      d0= 0;
3889
3890    if (d0 <  d)
3891    {
3892      m= gcdlcAlcB;
3893      newtonPoly= 1;
3894      G_m= 0;
3895      d= d0;
3896    }
3897
3898    H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
3899
3900    if (uni_lcoeff (H) == gcdlcAlcB)
3901    {
3902      cH= uni_content (H);
3903      ppH= H/cH;
3904      ppH /= Lc (ppH);
3905      DEBOUTLN (cerr, "ppH= " << ppH);
3906
3907      if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
3908        return N(gcdcAcB*ppH);
3909    }
3910    G_m= H;
3911    newtonPoly= newtonPoly*(x - random_element);
3912    m= m*(x - random_element);
3913    if (!find (l, random_element))
3914      l.append (random_element);
3915
3916    if ((getCharacteristic() > 3 && size (skeleton) < 100))
3917    {
3918      CFArray Monoms;
3919      CFArray* coeffMonoms= NULL;
3920
3921      do //second do
3922      {
3923        if (inextension)
3924          random_element= randomElement (m, alpha, l, fail);
3925        else
3926          random_element= FpRandomElement (m, l, fail);
3927        if (random_element == 0 && !fail)
3928        {
3929          if (!find (l, random_element))
3930            l.append (random_element);
3931          continue;
3932        }
3933
3934        bool sparseFail= false;
3935        if (!fail && !inextension)
3936        {
3937          CFList list;
3938          TIMING_START (gcd_recursion);
3939          if (LC (skeleton).inCoeffDomain())
3940            G_random_element=
3941            monicSparseInterpol(ppA(random_element, x), ppB (random_element, x),
3942                                skeleton, Variable(1), sparseFail, coeffMonoms,
3943                                Monoms);
3944          else
3945            G_random_element=
3946            nonMonicSparseInterpol(ppA(random_element,x), ppB(random_element,x),
3947                                    skeleton, Variable (1), sparseFail,
3948                                    coeffMonoms, Monoms);
3949          TIMING_END_AND_PRINT (gcd_recursion,
3950                                "time for recursive call: ");
3951          DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3952        }
3953        else if (!fail && inextension)
3954        {
3955          CFList list;
3956          TIMING_START (gcd_recursion);
3957          if (LC (skeleton).inCoeffDomain())
3958            G_random_element=
3959            monicSparseInterpol(ppA (random_element,x), ppB (random_element, x),
3960                                skeleton, alpha, sparseFail, coeffMonoms,
3961                                Monoms);
3962          else
3963            G_random_element=
3964            nonMonicSparseInterpol(ppA(random_element,x), ppB(random_element,x),
3965                                   skeleton, alpha, sparseFail, coeffMonoms,
3966                                   Monoms);
3967          TIMING_END_AND_PRINT (gcd_recursion,
3968                                "time for recursive call: ");
3969          DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3970        }
3971        else if (fail && !inextension)
3972        {
3973          source= CFList();
3974          dest= CFList();
3975          CFList list;
3976          CanonicalForm mipo;
3977          int deg= 2;
3978          do
3979          {
3980            mipo= randomIrredpoly (deg, x);
3981            alpha= rootOf (mipo);
3982            inextension= true;
3983            fail= false;
3984            random_element= randomElement (m, alpha, l, fail);
3985            deg++;
3986          } while (fail);
3987          V_buf= alpha;
3988          list= CFList();
3989          TIMING_START (gcd_recursion);
3990          if (LC (skeleton).inCoeffDomain())
3991            G_random_element=
3992            monicSparseInterpol (ppA (random_element,x), ppB (random_element,x),
3993                                skeleton, alpha, sparseFail, coeffMonoms,
3994                                Monoms);
3995          else
3996            G_random_element=
3997            nonMonicSparseInterpol(ppA(random_element,x), ppB(random_element,x),
3998                                   skeleton, alpha, sparseFail, coeffMonoms,
3999                                   Monoms);
4000          TIMING_END_AND_PRINT (gcd_recursion,
4001                                "time for recursive call: ");
4002          DEBOUTLN (cerr, "G_random_element= " << G_random_element);
4003        }
4004        else if (fail && inextension)
4005        {
4006          source= CFList();
4007          dest= CFList();
4008
4009          Variable V_buf3= V_buf;
4010          V_buf= chooseExtension (V_buf);
4011          bool prim_fail= false;
4012          Variable V_buf2;
4013          CanonicalForm prim_elem, im_prim_elem;
4014          prim_elem= primitiveElement (alpha, V_buf2, prim_fail);
4015
4016          if (V_buf3 != alpha)
4017          {
4018            m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
4019            G_m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
4020            newtonPoly= mapDown (newtonPoly, prim_elem, im_prim_elem, alpha,
4021                                 source, dest);
4022            ppA= mapDown (ppA, prim_elem, im_prim_elem, alpha, source, dest);
4023            ppB= mapDown (ppB, prim_elem, im_prim_elem, alpha, source, dest);
4024            gcdlcAlcB= mapDown (gcdlcAlcB, prim_elem, im_prim_elem, alpha,
4025                                source, dest);
4026            for (CFListIterator i= l; i.hasItem(); i++)
4027              i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, alpha,
4028                                    source, dest);
4029          }
4030
4031          ASSERT (!prim_fail, "failure in integer factorizer");
4032          if (prim_fail)
4033            ; //ERROR
4034          else
4035            im_prim_elem= mapPrimElem (prim_elem, alpha, V_buf);
4036
4037          DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (alpha));
4038          DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf2));
4039
4040          for (CFListIterator i= l; i.hasItem(); i++)
4041            i.getItem()= mapUp (i.getItem(), alpha, V_buf, prim_elem,
4042                                im_prim_elem, source, dest);
4043          m= mapUp (m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
4044          G_m= mapUp (G_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
4045          newtonPoly= mapUp (newtonPoly, alpha, V_buf, prim_elem, im_prim_elem,
4046                              source, dest);
4047          ppA= mapUp (ppA, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
4048          ppB= mapUp (ppB, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
4049          gcdlcAlcB= mapUp (gcdlcAlcB, alpha, V_buf, prim_elem, im_prim_elem,
4050                            source, dest);
4051          fail= false;
4052          random_element= randomElement (m, V_buf, l, fail );
4053          DEBOUTLN (cerr, "fail= " << fail);
4054          CFList list;
4055          TIMING_START (gcd_recursion);
4056          if (LC (skeleton).inCoeffDomain())
4057            G_random_element=
4058            monicSparseInterpol (ppA (random_element, x), ppB (random_element, x),
4059                                skeleton, V_buf, sparseFail, coeffMonoms,
4060                                Monoms);
4061          else
4062            G_random_element=
4063            nonMonicSparseInterpol (ppA(random_element,x), ppB(random_element,x),
4064                                    skeleton, V_buf, sparseFail,
4065                                    coeffMonoms, Monoms);
4066          TIMING_END_AND_PRINT (gcd_recursion,
4067                                "time for recursive call: ");
4068          DEBOUTLN (cerr, "G_random_element= " << G_random_element);
4069        }
4070
4071        if (sparseFail)
4072          break;
4073
4074        if (!G_random_element.inCoeffDomain())
4075          d0= totaldegree (G_random_element, Variable(2),
4076                           Variable (G_random_element.level()));
4077        else
4078          d0= 0;
4079
4080        if (d0 == 0)
4081          return N(gcdcAcB);
4082        if (d0 >  d)
4083        {
4084          if (!find (l, random_element))
4085            l.append (random_element);
4086          continue;
4087        }
4088
4089        G_random_element=
4090        (gcdlcAlcB(random_element, x)/uni_lcoeff (G_random_element))
4091        * G_random_element;
4092
4093        if (!G_random_element.inCoeffDomain())
4094          d0= totaldegree (G_random_element, Variable(2),
4095                           Variable (G_random_element.level()));
4096        else
4097          d0= 0;
4098
4099        if (d0 <  d)
4100        {
4101          m= gcdlcAlcB;
4102          newtonPoly= 1;
4103          G_m= 0;
4104          d= d0;
4105        }
4106
4107        TIMING_START (newton_interpolation);
4108        H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
4109        TIMING_END_AND_PRINT (newton_interpolation,
4110                              "time for newton interpolation: ");
4111
4112        //termination test
4113        if (uni_lcoeff (H) == gcdlcAlcB)
4114        {
4115          cH= uni_content (H);
4116          ppH= H/cH;
4117          ppH /= Lc (ppH);
4118          DEBOUTLN (cerr, "ppH= " << ppH);
4119          if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
4120            return N(gcdcAcB*ppH);
4121        }
4122
4123        G_m= H;
4124        newtonPoly= newtonPoly*(x - random_element);
4125        m= m*(x - random_element);
4126        if (!find (l, random_element))
4127          l.append (random_element);
4128
4129      } while (1); //end of second do
4130    }
4131  } while (1); //end of first do
4132}
4133
4134static inline
4135int compress4EZGCD (const CanonicalForm& F, const CanonicalForm& G, CFMap & M,
4136                    CFMap & N, int& both_non_zero)
4137{
4138  int n= tmax (F.level(), G.level());
4139  int * degsf= new int [n + 1];
4140  int * degsg= new int [n + 1];
4141
4142  for (int i = 0; i <= n; i++)
4143    degsf[i]= degsg[i]= 0;
4144
4145  degsf= degrees (F, degsf);
4146  degsg= degrees (G, degsg);
4147
4148  both_non_zero= 0;
4149  int f_zero= 0;
4150  int g_zero= 0;
4151
4152  for (int i= 1; i <= n; i++)
4153  {
4154    if (degsf[i] != 0 && degsg[i] != 0)
4155    {
4156      both_non_zero++;
4157      continue;
4158    }
4159    if (degsf[i] == 0 && degsg[i] != 0 && i <= G.level())
4160    {
4161      f_zero++;
4162      continue;
4163    }
4164    if (degsg[i] == 0 && degsf[i] && i <= F.level())
4165    {
4166      g_zero++;
4167      continue;
4168    }
4169  }
4170
4171  if (both_non_zero == 0)
4172  {
4173    delete [] degsf;
4174    delete [] degsg;
4175    return 0;
4176  }
4177
4178  // map Variables which do not occur in both polynomials to higher levels
4179  int k= 1;
4180  int l= 1;
4181  for (int i= 1; i <= n; i++)
4182  {
4183    if (degsf[i] != 0 && degsg[i] == 0 && i <= F.level())
4184    {
4185      if (k + both_non_zero != i)
4186      {
4187        M.newpair (Variable (i), Variable (k + both_non_zero));
4188        N.newpair (Variable (k + both_non_zero), Variable (i));
4189      }
4190      k++;
4191    }
4192    if (degsf[i] == 0 && degsg[i] != 0 && i <= G.level())
4193    {
4194      if (l + g_zero + both_non_zero != i)
4195      {
4196        M.newpair (Variable (i), Variable (l + g_zero + both_non_zero));
4197        N.newpair (Variable (l + g_zero + both_non_zero), Variable (i));
4198      }
4199      l++;
4200    }
4201  }
4202
4203  // sort Variables x_{i} in decreasing order of
4204  // min(deg_{x_{i}}(f),deg_{x_{i}}(g))
4205  int m= tmin (F.level(), G.level());
4206  int max_min_deg;
4207  k= both_non_zero;
4208  l= 0;
4209  int i= 1;
4210  while (k > 0)
4211  {
4212    max_min_deg= tmin (degsf[i], degsg[i]);
4213    while (max_min_deg == 0)
4214    {
4215      i++;
4216      max_min_deg= tmin (degsf[i], degsg[i]);
4217    }
4218    for (int j= i + 1; j <= m; j++)
4219    {
4220      if ((tmin (degsf[j],degsg[j]) < max_min_deg) &&
4221          (tmin (degsf[j], degsg[j]) != 0))
4222      {
4223        max_min_deg= tmin (degsf[j], degsg[j]);
4224        l= j;
4225      }
4226    }
4227
4228    if (l != 0)
4229    {
4230      if (l != k)
4231      {
4232        M.newpair (Variable (l), Variable(k));
4233        N.newpair (Variable (k), Variable(l));
4234        degsf[l]= 0;
4235        degsg[l]= 0;
4236        l= 0;
4237      }
4238      else
4239      {
4240        degsf[l]= 0;
4241        degsg[l]= 0;
4242        l= 0;
4243      }
4244    }
4245    else if (l == 0)
4246    {
4247      if (i != k)
4248      {
4249        M.newpair (Variable (i), Variable (k));
4250        N.newpair (Variable (k), Variable (i));
4251        degsf[i]= 0;
4252        degsg[i]= 0;
4253      }
4254      else
4255      {
4256        degsf[i]= 0;
4257        degsg[i]= 0;
4258      }
4259      i++;
4260    }
4261    k--;
4262  }
4263
4264  delete [] degsf;
4265  delete [] degsg;
4266
4267  return both_non_zero;
4268}
4269
4270static inline
4271CanonicalForm myShift2Zero (const CanonicalForm& F, CFList& Feval,
4272                            const CFList& evaluation)
4273{
4274  CanonicalForm A= F;
4275  int k= 2;
4276  for (CFListIterator i= evaluation; i.hasItem(); i++, k++)
4277    A= A (Variable (k) + i.getItem(), k);
4278
4279  CanonicalForm buf= A;
4280  Feval= CFList();
4281  Feval.append (buf);
4282  for (k= evaluation.length() + 1; k > 2; k--)
4283  {
4284    buf= mod (buf, Variable (k));
4285    Feval.insert (buf);
4286  }
4287  return A;
4288}
4289
4290static inline
4291CanonicalForm myReverseShift (const CanonicalForm& F, const CFList& evaluation)
4292{
4293  int l= evaluation.length() + 1;
4294  CanonicalForm result= F;
4295  CFListIterator j= evaluation;
4296  for (int i= 2; i < l + 1; i++, j++)
4297  {
4298    if (F.level() < i)
4299      continue;
4300    result= result (Variable (i) - j.getItem(), i);
4301  }
4302  return result;
4303}
4304
4305static inline
4306Evaluation optimize4Lift (const CanonicalForm& F, CFMap & M,
4307                    CFMap & N, const Evaluation& A)
4308{
4309  int n= F.level();
4310  int * degsf= new int [n + 1];
4311
4312  for (int i = 0; i <= n; i++)
4313    degsf[i]= 0;
4314
4315  degsf= degrees (F, degsf);
4316
4317  Evaluation result= Evaluation (A.min(), A.max());
4318  ASSERT (A.min() == 2, "expected A.min() == 2");
4319  int max_deg;
4320  int k= n;
4321  int l= 1;
4322  int i= 2;
4323  int pos= 2;
4324  while (k > 1)
4325  {
4326    max_deg= degsf [i];
4327    while (max_deg == 0)
4328    {
4329      i++;
4330      max_deg= degsf [i];
4331    }
4332    l= i;
4333    for (int j= i + 1; j <=  n; j++)
4334    {
4335      if (degsf[j] > max_deg)
4336      {
4337        max_deg= degsf[j];
4338        l= j;
4339      }
4340    }
4341
4342    if (l <= n)
4343    {
4344      if (l != pos)
4345      {
4346        result.setValue (pos, A [l]);
4347        M.newpair (Variable (l), Variable (pos));
4348        N.newpair (Variable (pos), Variable (l));
4349        degsf[l]= 0;
4350        l= 2;
4351        if (k == 2 && n == 3)
4352        {
4353          result.setValue (l, A [pos]);
4354          M.newpair (Variable (pos), Variable (l));
4355          N.newpair (Variable (l), Variable (pos));
4356          degsf[pos]= 0;
4357        }
4358      }
4359      else
4360      {
4361        result.setValue (l, A [l]);
4362        degsf [l]= 0;
4363      }
4364    }
4365    pos++;
4366    k--;
4367    l= 2;
4368  }
4369
4370  delete [] degsf;
4371
4372  return result;
4373}
4374
4375static inline
4376int Hensel_P (const CanonicalForm & UU, CFArray & G, const Evaluation & AA,
4377              const CFArray& LeadCoeffs )
4378{
4379  CFList factors;
4380  factors.append (G[1]);
4381  factors.append (G[2]);
4382
4383  CFMap NN, MM;
4384  Evaluation A= optimize4Lift (UU, MM, NN, AA);
4385
4386  CanonicalForm U= MM (UU);
4387  CFArray LCs= CFArray (1,2);
4388  LCs [1]= MM (LeadCoeffs [1]);
4389  LCs [2]= MM (LeadCoeffs [2]);
4390
4391  CFList evaluation;
4392  for (int i= A.min(); i <= A.max(); i++)
4393    evaluation.append (A [i]);
4394  CFList UEval;
4395  CanonicalForm shiftedU= myShift2Zero (U, UEval, evaluation);
4396
4397  if (size (shiftedU)/getNumVars (U) > 500)
4398    return -1;
4399
4400  CFArray shiftedLCs= CFArray (2);
4401  CFList shiftedLCsEval1, shiftedLCsEval2;
4402  shiftedLCs[0]= myShift2Zero (LCs[1], shiftedLCsEval1, evaluation);
4403  shiftedLCs[1]= myShift2Zero (LCs[2], shiftedLCsEval2, evaluation);
4404  factors.insert (1);
4405  int liftBound= degree (UEval.getLast(), 2) + 1;
4406  CFArray Pi;
4407  CFMatrix M= CFMatrix (liftBound, factors.length() - 1);
4408  CFList diophant;
4409  CFArray lcs= CFArray (2);
4410  lcs [0]= shiftedLCsEval1.getFirst();
4411  lcs [1]= shiftedLCsEval2.getFirst();
4412  nonMonicHenselLift12 (UEval.getFirst(), factors, liftBound, Pi, diophant, M,
4413                        lcs, false);
4414
4415  for (CFListIterator i= factors; i.hasItem(); i++)
4416  {
4417    if (!fdivides (i.getItem(), UEval.getFirst()))
4418      return 0;
4419  }
4420
4421  int * liftBounds;
4422  bool noOneToOne= false;
4423  if (U.level() > 2)
4424  {
4425    liftBounds= new int [U.level() - 1]; /* index: 0.. U.level()-2 */
4426    liftBounds[0]= liftBound;
4427    for (int i= 1; i < U.level() - 1; i++)
4428      liftBounds[i]= degree (shiftedU, Variable (i + 2)) + 1;
4429    factors= nonMonicHenselLift2 (UEval, factors, liftBounds, U.level() - 1,
4430                                  false, shiftedLCsEval1, shiftedLCsEval2, Pi,
4431                                  diophant, noOneToOne);
4432    delete [] liftBounds;
4433    if (noOneToOne)
4434      return 0;
4435  }
4436  G[1]= factors.getFirst();
4437  G[2]= factors.getLast();
4438  G[1]= myReverseShift (G[1], evaluation);
4439  G[2]= myReverseShift (G[2], evaluation);
4440  G[1]= NN (G[1]);
4441  G[2]= NN (G[2]);
4442  return 1;
4443}
4444
4445static inline
4446bool findeval_P (const CanonicalForm & F, const CanonicalForm & G,
4447                 CanonicalForm & Fb, CanonicalForm & Gb, CanonicalForm & Db,
4448                 REvaluation & b, int delta, int degF, int degG, int maxeval,
4449                 int & count, int& k, int bound, int& l)
4450{
4451  if( count == 0 && delta != 0)
4452  {
4453    if( count++ > maxeval )
4454      return false;
4455  }
4456  if (count > 0)
4457  {
4458    b.nextpoint(k);
4459    if (k == 0)
4460      k++;
4461    l++;
4462    if (l > bound)
4463    {
4464      l= 1;
4465      k++;
4466      if (k > tmax (F.level(), G.level()) - 1)
4467        return false;
4468      b.nextpoint (k);
4469    }
4470    if (count++ > maxeval)
4471      return false;
4472  }
4473  while( true )
4474  {
4475    Fb = b( F );
4476    if( degree( Fb, 1 ) == degF )
4477    {
4478      Gb = b( G );
4479      if( degree( Gb, 1 ) == degG )
4480      {
4481        Db = gcd( Fb, Gb );
4482        if( delta > 0 )
4483        {
4484          if( degree( Db, 1 ) <= delta )
4485            return true;
4486        }
4487        else
4488          return true;
4489      }
4490    }
4491    if (k == 0)
4492      k++;
4493    b.nextpoint(k);
4494    l++;
4495    if (l > bound)
4496    {
4497      l= 1;
4498      k++;
4499      if (k > tmax (F.level(), G.level()) - 1)
4500        return false;
4501      b.nextpoint (k);
4502    }
4503    if( count++ > maxeval )
4504      return false;
4505  }
4506}
4507
4508// parameters for heuristic
4509static int maxNumEval= 200;
4510static int sizePerVars1= 500; //try dense gcd if size/#variables is bigger
4511
4512/// Extended Zassenhaus GCD for finite fields
4513CanonicalForm EZGCD_P( const CanonicalForm & FF, const CanonicalForm & GG )
4514{
4515  if (FF.isZero() && degree(GG) > 0) return GG/Lc(GG);
4516  else if (GG.isZero() && degree (FF) > 0) return FF/Lc(FF);
4517  else if (FF.isZero() && GG.isZero()) return FF.genOne();
4518  if (FF.inBaseDomain() || GG.inBaseDomain()) return FF.genOne();
4519  if (FF.isUnivariate() && fdivides(FF, GG)) return FF/Lc(FF);
4520  if (GG.isUnivariate() && fdivides(GG, FF)) return GG/Lc(GG);
4521  if (FF == GG) return FF/Lc(FF);
4522
4523  CanonicalForm F, G, f, g, d, Fb, Gb, Db, Fbt, Gbt, Dbt, B0, B, D0, lcF, lcG,
4524                lcD;
4525  CFArray DD( 1, 2 ), lcDD( 1, 2 );
4526  int degF, degG, delta, count;
4527  int maxeval;
4528  maxeval= tmin((getCharacteristic()/
4529                (int)(ilog2(getCharacteristic())*log2exp))*2, maxNumEval);
4530  count= 0; // number of eval. used
4531  REvaluation b, bt;
4532  int gcdfound = 0;
4533  Variable x = Variable(1);
4534
4535  F= FF;
4536  G= GG;
4537
4538  CFMap M,N;
4539  int smallestDegLev;
4540  TIMING_START (ez_p_compress)
4541  int best_level= compress4EZGCD (F, G, M, N, smallestDegLev);
4542
4543  if (best_level == 0) return G.genOne();
4544
4545  F= M (F);
4546  G= M (G);
4547  TIMING_END_AND_PRINT (ez_p_compress, "time for compression in EZ_P: ")
4548
4549  TIMING_START (ez_p_content)
4550  f = content( F, x ); g = content( G, x );
4551  d = gcd( f, g );
4552  F /= f; G /= g;
4553  TIMING_END_AND_PRINT (ez_p_content, "time to extract content in EZ_P: ")
4554
4555  if( F.isUnivariate() && G.isUnivariate() )
4556  {
4557    if( F.mvar() == G.mvar() )
4558      d *= gcd( F, G );
4559    return N (d);
4560  }
4561
4562  int maxNumVars= tmax (getNumVars (F), getNumVars (G));
4563  Variable a, oldA;
4564  int sizeF= size (F);
4565  int sizeG= size (G);
4566
4567  if (sizeF/maxNumVars > sizePerVars1 && sizeG/maxNumVars > sizePerVars1)
4568  {
4569    if (hasFirstAlgVar (F, a) || hasFirstAlgVar (G, a))
4570      return N (d*GCD_Fp_extension (F, G, a));
4571    else if (CFFactory::gettype() == GaloisFieldDomain)
4572      return N (d*GCD_GF (F, G));
4573    else
4574      return N (d*GCD_small_p (F, G));
4575  }
4576
4577  int dummy= 0;
4578  if( gcd_test_one( F, G, false, dummy ) )
4579  {
4580    return N (d);
4581  }
4582
4583  bool passToGF= false;
4584  bool extOfExt= false;
4585  int p= getCharacteristic();
4586  bool algExtension= (hasFirstAlgVar(F,a) || hasFirstAlgVar(G,a));
4587  int k= 1;
4588  CanonicalForm primElem, imPrimElem;
4589  CFList source, dest;
4590  if (p < 50 && CFFactory::gettype() != GaloisFieldDomain && !algExtension)
4591  {
4592    if (p == 2)
4593      setCharacteristic (2, 12, 'Z');
4594    else if (p == 3)
4595      setCharacteristic (3, 4, 'Z');
4596    else if (p == 5 || p == 7)
4597      setCharacteristic (p, 3, 'Z');
4598    else
4599      setCharacteristic (p, 2, 'Z');
4600    passToGF= true;
4601    F= F.mapinto();
4602    G= G.mapinto();
4603    maxeval= 2*ipower (p, getGFDegree());
4604  }
4605  else if (CFFactory::gettype() == GaloisFieldDomain &&
4606           ipower (p , getGFDegree()) < 50)
4607  {
4608    k= getGFDegree();
4609    if (ipower (p, 2*k) > 50)
4610      setCharacteristic (p, 2*k, gf_name);
4611    else
4612      setCharacteristic (p, 3*k, gf_name);
4613    F= GFMapUp (F, k);
4614    G= GFMapUp (G, k);
4615    maxeval= tmin (2*ipower (p, getGFDegree()), maxNumEval);
4616  }
4617  else if (p < 50 && algExtension && !CFFactory::gettype() == GaloisFieldDomain)
4618  {
4619    int d= degree (getMipo (a));
4620    oldA= a;
4621    Variable v2;
4622    if (p == 2 && d < 6)
4623    {
4624      if (fac_NTL_char != p)
4625      {
4626        fac_NTL_char= p;
4627        zz_p::init (p);
4628      }
4629      bool primFail= false;
4630      Variable vBuf;
4631      primElem= primitiveElement (a, vBuf, primFail);
4632      ASSERT (!primFail, "failure in integer factorizer");
4633      if (d < 3)
4634      {
4635        zz_pX NTLIrredpoly;
4636        BuildIrred (NTLIrredpoly, d*3);
4637        CanonicalForm newMipo= convertNTLzzpX2CF (NTLIrredpoly, Variable (1));
4638        v2= rootOf (newMipo);
4639      }
4640      else
4641      {
4642        zz_pX NTLIrredpoly;
4643        BuildIrred (NTLIrredpoly, d*2);
4644        CanonicalForm newMipo= convertNTLzzpX2CF (NTLIrredpoly, Variable (1));
4645        v2= rootOf (newMipo);
4646      }
4647      imPrimElem= mapPrimElem (primElem, a, v2);
4648      extOfExt= true;
4649    }
4650    else if ((p == 3 && d < 4) || ((p == 5 || p == 7) && d < 3))
4651    {
4652      if (fac_NTL_char != p)
4653      {
4654        fac_NTL_char= p;
4655        zz_p::init (p);
4656      }
4657      bool primFail= false;
4658      Variable vBuf;
4659      primElem= primitiveElement (a, vBuf, primFail);
4660      ASSERT (!primFail, "failure in integer factorizer");
4661      zz_pX NTLIrredpoly;
4662      BuildIrred (NTLIrredpoly, d*2);
4663      CanonicalForm newMipo= convertNTLzzpX2CF (NTLIrredpoly, Variable (1));
4664      v2= rootOf (newMipo);
4665      imPrimElem= mapPrimElem (primElem, a, v2);
4666      extOfExt= true;
4667    }
4668    if (extOfExt)
4669    {
4670      maxeval= tmin (2*ipower (p, degree (getMipo (v2))), maxNumEval);
4671      F= mapUp (F, a, v2, primElem, imPrimElem, source, dest);
4672      G= mapUp (G, a, v2, primElem, imPrimElem, source, dest);
4673      a= v2;
4674    }
4675  }
4676
4677  lcF = LC( F, x ); lcG = LC( G, x );
4678  lcD = gcd( lcF, lcG );
4679
4680  delta = 0;
4681  degF = degree( F, x ); degG = degree( G, x );
4682
4683  if(hasFirstAlgVar(G,a))
4684    b = REvaluation( 2, tmax(F.level(), G.level()), AlgExtRandomF( a ) );
4685  else
4686  { // both not in extension given by algebraic variable
4687    if (CFFactory::gettype() != GaloisFieldDomain)
4688      b = REvaluation( 2, tmax(F.level(), G.level()), FFRandom() );
4689    else
4690      b = REvaluation( 2, tmax(F.level(), G.level()), GFRandom() );
4691  }
4692
4693  CanonicalForm cand, contcand;
4694  CanonicalForm result;
4695  int o, t;
4696  o= 0;
4697  t= 1;
4698  int goodPointCount= 0;
4699  while( !gcdfound )
4700  {
4701    TIMING_START (ez_p_eval);
4702    if( !findeval_P( F, G, Fb, Gb, Db, b, delta, degF, degG, maxeval, count, o,
4703         maxeval/maxNumVars, t ))
4704    { // too many eval. used --> try another method
4705      Off (SW_USE_EZGCD_P);
4706      result= gcd (F,G);
4707      On (SW_USE_EZGCD_P);
4708      if (passToGF)
4709      {
4710        CanonicalForm mipo= gf_mipo;
4711        setCharacteristic (p);
4712        Variable alpha= rootOf (mipo.mapinto());
4713        result= GF2FalphaRep (result, alpha);
4714      }
4715      if (k > 1)
4716      {
4717        result= GFMapDown (result, k);
4718        setCharacteristic (p, k, gf_name);
4719      }
4720      if (extOfExt)
4721        result= mapDown (result, primElem, imPrimElem, oldA, dest, source);
4722      return N (d*result);
4723    }
4724    TIMING_END_AND_PRINT (ez_p_eval, "time for eval point search in EZ_P1: ");
4725    delta = degree( Db );
4726    if( delta == 0 )
4727    {
4728      if (passToGF)
4729        setCharacteristic (p);
4730      if (k > 1)
4731        setCharacteristic (p, k, gf_name);
4732      return N (d);
4733    }
4734    while( true )
4735    {
4736      bt = b;
4737      TIMING_START (ez_p_eval);
4738      if( !findeval_P(F,G,Fbt,Gbt,Dbt, bt, delta, degF, degG, maxeval, count, o,
4739           maxeval/maxNumVars, t ))
4740      { // too many eval. used --> try another method
4741        Off (SW_USE_EZGCD_P);
4742        result= gcd (F,G);
4743        On (SW_USE_EZGCD_P);
4744        if (passToGF)
4745        {
4746          CanonicalForm mipo= gf_mipo;
4747          setCharacteristic (p);
4748          Variable alpha= rootOf (mipo.mapinto());
4749          result= GF2FalphaRep (result, alpha);
4750        }
4751        if (k > 1)
4752        {
4753          result= GFMapDown (result, k);
4754          setCharacteristic (p, k, gf_name);
4755        }
4756        if (extOfExt)
4757          result= mapDown (result, primElem, imPrimElem, oldA, dest, source);
4758        return N (d*result);
4759      }
4760      TIMING_END_AND_PRINT (ez_p_eval, "time for eval point search in EZ_P2: ");
4761      int dd = degree( Dbt );
4762      if( dd == 0 )
4763      {
4764        if (passToGF)
4765          setCharacteristic (p);
4766        if (k > 1)
4767          setCharacteristic (p, k, gf_name);
4768        return N (d);
4769      }
4770      if( dd == delta )
4771      {
4772        goodPointCount++;
4773        if (goodPointCount == 5)
4774          break;
4775      }
4776      if( dd < delta )
4777      {
4778        goodPointCount= 0;
4779        delta = dd;
4780        b = bt;
4781        Db = Dbt; Fb = Fbt; Gb = Gbt;
4782      }
4783      if (delta == degF)
4784      {
4785        if (degF <= degG && fdivides (F, G))
4786        {
4787          if (passToGF)
4788          {
4789            CanonicalForm mipo= gf_mipo;
4790            setCharacteristic (p);
4791            Variable alpha= rootOf (mipo.mapinto());
4792            F= GF2FalphaRep (F, alpha);
4793          }
4794          if (k > 1)
4795          {
4796            F= GFMapDown (F, k);
4797            setCharacteristic (p, k, gf_name);
4798          }
4799          if (extOfExt)
4800            F= mapDown (F, primElem, imPrimElem, oldA, dest, source);
4801          return N (d*F);
4802        }
4803        else
4804          delta--;
4805      }
4806      else if (delta == degG)
4807      {
4808        if (degG <= degF && fdivides (G, F))
4809        {
4810          if (passToGF)
4811          {
4812            CanonicalForm mipo= gf_mipo;
4813            setCharacteristic (p);
4814            Variable alpha= rootOf (mipo.mapinto());
4815            G= GF2FalphaRep (G, alpha);
4816          }
4817          if (k > 1)
4818          {
4819            G= GFMapDown (G, k);
4820            setCharacteristic (p, k, gf_name);
4821          }
4822          if (extOfExt)
4823            G= mapDown (G, primElem, imPrimElem, oldA, dest, source);
4824          return N (d*G);
4825        }
4826        else
4827          delta--;
4828      }
4829    }
4830    if( delta != degF && delta != degG )
4831    {
4832      bool B_is_F;
4833      CanonicalForm xxx1, xxx2;
4834      CanonicalForm buf;
4835      DD[1] = Fb / Db;
4836      buf= Gb/Db;
4837      xxx1 = gcd( DD[1], Db );
4838      xxx2 = gcd( buf, Db );
4839      if (((xxx1.inCoeffDomain() && xxx2.inCoeffDomain()) &&
4840          (size (F) <= size (G)))
4841          || (xxx1.inCoeffDomain() && !xxx2.inCoeffDomain()))
4842      {
4843        B = F;
4844        DD[2] = Db;
4845        lcDD[1] = lcF;
4846        lcDD[2] = lcD;
4847        B_is_F = true;
4848      }
4849      else if (((xxx1.inCoeffDomain() && xxx2.inCoeffDomain()) &&
4850               (size (G) < size (F)))
4851               || (!xxx1.inCoeffDomain() && xxx2.inCoeffDomain()))
4852      {
4853        DD[1] = buf;
4854        B = G;
4855        DD[2] = Db;
4856        lcDD[1] = lcG;
4857        lcDD[2] = lcD;
4858        B_is_F = false;
4859      }
4860      else // special case handling
4861      {
4862        Off (SW_USE_EZGCD_P);
4863        result= gcd (F,G);
4864        On (SW_USE_EZGCD_P);
4865        if (passToGF)
4866        {
4867          CanonicalForm mipo= gf_mipo;
4868          setCharacteristic (p);
4869          Variable alpha= rootOf (mipo.mapinto());
4870          result= GF2FalphaRep (result, alpha);
4871        }
4872        if (k > 1)
4873        {
4874          result= GFMapDown (result, k);
4875          setCharacteristic (p, k, gf_name);
4876        }
4877        if (extOfExt)
4878          result= mapDown (result, primElem, imPrimElem, oldA, dest, source);
4879        return N (d*result);
4880      }
4881      DD[2] = DD[2] * ( b( lcDD[2] ) / lc( DD[2] ) );
4882      DD[1] = DD[1] * ( b( lcDD[1] ) / lc( DD[1] ) );
4883
4884      if (size (B*lcDD[2])/maxNumVars > sizePerVars1)
4885      {
4886        if (algExtension)
4887        {
4888          result= GCD_Fp_extension (F, G, a);
4889          if (extOfExt)
4890            result= mapDown (result, primElem, imPrimElem, oldA, dest, source);
4891          return N (d*result);
4892        }
4893        if (CFFactory::gettype() == GaloisFieldDomain)
4894        {
4895          result= GCD_GF (F, G);
4896          if (passToGF)
4897          {
4898            CanonicalForm mipo= gf_mipo;
4899            setCharacteristic (p);
4900            Variable alpha= rootOf (mipo.mapinto());
4901            result= GF2FalphaRep (result, alpha);
4902          }
4903          if (k > 1)
4904          {
4905            result= GFMapDown (result, k);
4906            setCharacteristic (p, k, gf_name);
4907          }
4908          return N (d*result);
4909        }
4910        else
4911          return N (d*GCD_small_p (F,G));
4912      }
4913
4914      TIMING_START (ez_p_hensel_lift);
4915      gcdfound= Hensel_P (B*lcD, DD, b, lcDD);
4916      TIMING_END_AND_PRINT (ez_p_hensel_lift, "time for Hensel lift in EZ_P: ");
4917
4918      if (gcdfound == -1) //things became dense
4919      {
4920        if (algExtension)
4921        {
4922          result= GCD_Fp_extension (F, G, a);
4923          if (extOfExt)
4924            result= mapDown (result, primElem, imPrimElem, oldA, dest, source);
4925          return N (d*result);
4926        }
4927        if (CFFactory::gettype() == GaloisFieldDomain)
4928        {
4929          result= GCD_GF (F, G);
4930          if (passToGF)
4931          {
4932            CanonicalForm mipo= gf_mipo;
4933            setCharacteristic (p);
4934            Variable alpha= rootOf (mipo.mapinto());
4935            result= GF2FalphaRep (result, alpha);
4936          }
4937          if (k > 1)
4938          {
4939            result= GFMapDown (result, k);
4940            setCharacteristic (p, k, gf_name);
4941          }
4942          return N (d*result);
4943        }
4944        else
4945          return N (d*GCD_small_p (F,G));
4946      }
4947
4948      if (gcdfound == 1)
4949      {
4950        TIMING_START (termination_test);
4951        contcand= content (DD[2], Variable (1));
4952        cand = DD[2] / contcand;
4953        if (B_is_F)
4954          gcdfound = fdivides( cand, G ) && cand*(DD[1]/(lcD/contcand)) == F;
4955        else
4956          gcdfound = fdivides( cand, F ) && cand*(DD[1]/(lcD/contcand)) == G;
4957        TIMING_END_AND_PRINT (termination_test,
4958                              "time for termination test EZ_P: ");
4959
4960        if (passToGF && gcdfound)
4961        {
4962          CanonicalForm mipo= gf_mipo;
4963          setCharacteristic (p);
4964          Variable alpha= rootOf (mipo.mapinto());
4965          cand= GF2FalphaRep (cand, alpha);
4966        }
4967        if (k > 1 && gcdfound)
4968        {
4969          cand= GFMapDown (cand, k);
4970          setCharacteristic (p, k, gf_name);
4971        }
4972        if (extOfExt && gcdfound)
4973          cand= mapDown (cand, primElem, imPrimElem, oldA, dest, source);
4974      }
4975    }
4976    delta--;
4977    goodPointCount= 0;
4978  }
4979  return N (d*cand);
4980}
4981
4982
4983#endif
Note: See TracBrowser for help on using the repository browser.