source: git/factory/cf_gcd_smallp.cc @ a52291

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