source: git/factory/cf_gcd_smallp.cc @ e26667

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