source: git/factory/cf_gcd_smallp.cc @ ee668e

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