source: git/factory/cf_gcd_smallp.cc @ 4704674

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