source: git/factory/cf_gcd_smallp.cc @ 72bfc8

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