source: git/factory/cf_gcd_smallp.cc @ 3a3393b

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