source: git/factory/cf_gcd_smallp.cc @ 5c0bf0

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