source: git/factory/cf_gcd_smallp.cc @ 82e0a7

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