source: git/factory/cf_gcd_smallp.cc @ 015711

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