source: git/factory/cf_gcd_smallp.cc @ ea5ff1d

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