source: git/factory/cf_gcd_smallp.cc @ ea095d

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