source: git/factory/cf_gcd_smallp.cc @ 050d1b

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