source: git/factory/cf_gcd_smallp.cc @ 0dff6bc

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