source: git/factory/cf_gcd_smallp.cc @ 1a9dc5

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