source: git/factory/cfModGcd.cc @ 7ecb780

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