source: git/factory/cf_gcd_smallp.cc @ 2a95b2

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