source: git/factory/cf_gcd_smallp.cc @ f7a4e9

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