source: git/factory/cf_gcd_smallp.cc @ 7c00d6

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