source: git/factory/cf_gcd_smallp.cc @ 7762549

spielwiese
Last change on this file since 7762549 was 60d63da, checked in by Martin Lee <martinlee84@…>, 11 years ago
fix: memory leak
  • Property mode set to 100644
File size: 127.0 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, Labahn
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/// computes a random monic irreducible univariate polynomial in x of random
1136/// degree < i
1137CanonicalForm
1138randomIrredpoly (int i, const Variable & x)
1139{
1140  int p= getCharacteristic();
1141  if (fac_NTL_char != p)
1142  {
1143    fac_NTL_char= p;
1144    zz_p::init (p);
1145  }
1146  zz_pX NTLirredpoly;
1147  CanonicalForm CFirredpoly;
1148  BuildIrred (NTLirredpoly, i + 1);
1149  CFirredpoly= convertNTLzzpX2CF (NTLirredpoly, x);
1150  return CFirredpoly;
1151}
1152
1153static inline
1154CanonicalForm
1155FpRandomElement (const CanonicalForm& F, CFList& list, bool& fail)
1156{
1157  fail= false;
1158  Variable x= F.mvar();
1159  FFRandom genFF;
1160  CanonicalForm random;
1161  int p= getCharacteristic();
1162  int bound= p;
1163  do
1164  {
1165    if (list.length() == bound)
1166    {
1167      fail= true;
1168      break;
1169    }
1170    if (list.length() < 1)
1171      random= 0;
1172    else
1173    {
1174      random= genFF.generate();
1175      while (find (list, random))
1176        random= genFF.generate();
1177    }
1178    if (F (random, x) == 0)
1179    {
1180      list.append (random);
1181      continue;
1182    }
1183  } while (find (list, random));
1184  return random;
1185}
1186
1187CanonicalForm
1188GCD_small_p (const CanonicalForm& F, const CanonicalForm&  G,
1189             CanonicalForm& coF, CanonicalForm& coG,
1190             bool& topLevel, CFList& l);
1191
1192CanonicalForm
1193GCD_small_p (const CanonicalForm& F, const CanonicalForm& G,
1194             bool& topLevel, CFList& l)
1195{
1196  CanonicalForm dummy1, dummy2;
1197  CanonicalForm result= GCD_small_p (F, G, dummy1, dummy2, topLevel, l);
1198  return result;
1199}
1200
1201CanonicalForm
1202GCD_small_p (const CanonicalForm& F, const CanonicalForm&  G,
1203             CanonicalForm& coF, CanonicalForm& coG,
1204             bool& topLevel, CFList& l)
1205{
1206  CanonicalForm A= F;
1207  CanonicalForm B= G;
1208  if (F.isZero() && degree(G) > 0)
1209  {
1210    coF= 0;
1211    coG= Lc (G);
1212    return G/Lc(G);
1213  }
1214  else if (G.isZero() && degree (F) > 0)
1215  {
1216    coF= Lc (F);
1217    coG= 0;
1218    return F/Lc(F);
1219  }
1220  else if (F.isZero() && G.isZero())
1221  {
1222    coF= coG= 0;
1223    return F.genOne();
1224  }
1225  if (F.inBaseDomain() || G.inBaseDomain())
1226  {
1227    coF= F;
1228    coG= G;
1229    return F.genOne();
1230  }
1231  if (F.isUnivariate() && fdivides(F, G, coG))
1232  {
1233    coF= Lc (F);
1234    return F/Lc(F);
1235  }
1236  if (G.isUnivariate() && fdivides(G, F, coF))
1237  {
1238    coG= Lc (G);
1239    return G/Lc(G);
1240  }
1241  if (F == G)
1242  {
1243    coF= coG= Lc (F);
1244    return F/Lc(F);
1245  }
1246  CFMap M,N;
1247  int best_level= myCompress (A, B, M, N, topLevel);
1248
1249  if (best_level == 0)
1250  {
1251    coF= F;
1252    coG= G;
1253    return B.genOne();
1254  }
1255
1256  A= M(A);
1257  B= M(B);
1258
1259  Variable x= Variable (1);
1260
1261  //univariate case
1262  if (A.isUnivariate() && B.isUnivariate())
1263  {
1264    CanonicalForm result= gcd (A, B);
1265    coF= N (A/result);
1266    coG= N (B/result);
1267    return N (result);
1268  }
1269
1270  CanonicalForm cA, cB;    // content of A and B
1271  CanonicalForm ppA, ppB;    // primitive part of A and B
1272  CanonicalForm gcdcAcB;
1273
1274  cA = uni_content (A);
1275  cB = uni_content (B);
1276  gcdcAcB= gcd (cA, cB);
1277  ppA= A/cA;
1278  ppB= B/cB;
1279
1280  CanonicalForm lcA, lcB;  // leading coefficients of A and B
1281  CanonicalForm gcdlcAlcB;
1282  lcA= uni_lcoeff (ppA);
1283  lcB= uni_lcoeff (ppB);
1284
1285  gcdlcAlcB= gcd (lcA, lcB);
1286
1287  int d= totaldegree (ppA, Variable (2), Variable (ppA.level()));
1288  int d0;
1289
1290  if (d == 0)
1291  {
1292    coF= N (ppA*(cA/gcdcAcB));
1293    coG= N (ppB*(cB/gcdcAcB));
1294    return N(gcdcAcB);
1295  }
1296
1297  d0= totaldegree (ppB, Variable (2), Variable (ppB.level()));
1298
1299  if (d0 < d)
1300    d= d0;
1301
1302  if (d == 0)
1303  {
1304    coF= N (ppA*(cA/gcdcAcB));
1305    coG= N (ppB*(cB/gcdcAcB));
1306    return N(gcdcAcB);
1307  }
1308
1309  CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH;
1310  CanonicalForm newtonPoly, coF_random_element, coG_random_element,
1311                coF_m, coG_m, ppCoF, ppCoG;
1312
1313  newtonPoly= 1;
1314  m= gcdlcAlcB;
1315  H= 0;
1316  coF= 0;
1317  coG= 0;
1318  G_m= 0;
1319  Variable alpha, V_buf;
1320  bool fail= false;
1321  bool inextension= false;
1322  topLevel= false;
1323  CFList source, dest;
1324  int bound1= degree (ppA, 1);
1325  int bound2= degree (ppB, 1);
1326  do
1327  {
1328    if (inextension)
1329      random_element= randomElement (m*lcA*lcB, V_buf, l, fail);
1330    else
1331      random_element= FpRandomElement (m*lcA*lcB, l, fail);
1332
1333    if (!fail && !inextension)
1334    {
1335      CFList list;
1336      TIMING_START (gcd_recursion);
1337      G_random_element=
1338      GCD_small_p (ppA (random_element,x), ppB (random_element,x),
1339                   coF_random_element, coG_random_element, topLevel,
1340                   list);
1341      TIMING_END_AND_PRINT (gcd_recursion,
1342                            "time for recursive call: ");
1343      DEBOUTLN (cerr, "G_random_element= " << G_random_element);
1344    }
1345    else if (!fail && inextension)
1346    {
1347      CFList list;
1348      TIMING_START (gcd_recursion);
1349      G_random_element=
1350      GCD_Fp_extension (ppA (random_element, x), ppB (random_element, x),
1351                        coF_random_element, coG_random_element, V_buf,
1352                        list, topLevel);
1353      TIMING_END_AND_PRINT (gcd_recursion,
1354                            "time for recursive call: ");
1355      DEBOUTLN (cerr, "G_random_element= " << G_random_element);
1356    }
1357    else if (fail && !inextension)
1358    {
1359      source= CFList();
1360      dest= CFList();
1361      CFList list;
1362      CanonicalForm mipo;
1363      int deg= 2;
1364      do {
1365        mipo= randomIrredpoly (deg, x);
1366        alpha= rootOf (mipo);
1367        inextension= true;
1368        fail= false;
1369        random_element= randomElement (m*lcA*lcB, alpha, l, fail);
1370        deg++;
1371      } while (fail);
1372      list= CFList();
1373      V_buf= alpha;
1374      TIMING_START (gcd_recursion);
1375      G_random_element=
1376      GCD_Fp_extension (ppA (random_element, x), ppB (random_element, x),
1377                        coF_random_element, coG_random_element, alpha,
1378                        list, topLevel);
1379      TIMING_END_AND_PRINT (gcd_recursion,
1380                            "time for recursive call: ");
1381      DEBOUTLN (cerr, "G_random_element= " << G_random_element);
1382    }
1383    else if (fail && inextension)
1384    {
1385      source= CFList();
1386      dest= CFList();
1387
1388      Variable V_buf3= V_buf;
1389      V_buf= chooseExtension (V_buf);
1390      bool prim_fail= false;
1391      Variable V_buf2;
1392      CanonicalForm prim_elem, im_prim_elem;
1393      prim_elem= primitiveElement (alpha, V_buf2, prim_fail);
1394
1395      if (V_buf3 != alpha)
1396      {
1397        m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
1398        G_m= mapDown (G_m, prim_elem, im_prim_elem, alpha, source, dest);
1399        coF_m= mapDown (coF_m, prim_elem, im_prim_elem, alpha, source, dest);
1400        coG_m= mapDown (coG_m, prim_elem, im_prim_elem, alpha, source, dest);
1401        newtonPoly= mapDown (newtonPoly, prim_elem, im_prim_elem, alpha,
1402                             source, dest);
1403        ppA= mapDown (ppA, prim_elem, im_prim_elem, alpha, source, dest);
1404        ppB= mapDown (ppB, prim_elem, im_prim_elem, alpha, source, dest);
1405        gcdlcAlcB= mapDown (gcdlcAlcB, prim_elem, im_prim_elem, alpha, source,
1406                            dest);
1407        lcA= mapDown (lcA, prim_elem, im_prim_elem, alpha, source, dest);
1408        lcB= mapDown (lcB, prim_elem, im_prim_elem, alpha, source, dest);
1409        for (CFListIterator i= l; i.hasItem(); i++)
1410          i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, alpha,
1411                                source, dest);
1412      }
1413
1414      ASSERT (!prim_fail, "failure in integer factorizer");
1415      if (prim_fail)
1416        ; //ERROR
1417      else
1418        im_prim_elem= mapPrimElem (prim_elem, alpha, V_buf);
1419
1420      DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (alpha));
1421      DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf2));
1422
1423      for (CFListIterator i= l; i.hasItem(); i++)
1424        i.getItem()= mapUp (i.getItem(), alpha, V_buf, prim_elem,
1425                             im_prim_elem, source, dest);
1426      m= mapUp (m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
1427      G_m= mapUp (G_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
1428      coF_m= mapUp (coF_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
1429      coG_m= mapUp (coG_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
1430      newtonPoly= mapUp (newtonPoly, alpha, V_buf, prim_elem, im_prim_elem,
1431                          source, dest);
1432      ppA= mapUp (ppA, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
1433      ppB= mapUp (ppB, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
1434      gcdlcAlcB= mapUp (gcdlcAlcB, alpha, V_buf, prim_elem, im_prim_elem,
1435                         source, dest);
1436      lcA= mapUp (lcA, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
1437      lcB= mapUp (lcB, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
1438      fail= false;
1439      random_element= randomElement (m*lcA*lcB, V_buf, l, fail );
1440      DEBOUTLN (cerr, "fail= " << fail);
1441      CFList list;
1442      TIMING_START (gcd_recursion);
1443      G_random_element=
1444      GCD_Fp_extension (ppA (random_element, x), ppB (random_element, x),
1445                        coF_random_element, coG_random_element, V_buf,
1446                        list, topLevel);
1447      TIMING_END_AND_PRINT (gcd_recursion,
1448                            "time for recursive call: ");
1449      DEBOUTLN (cerr, "G_random_element= " << G_random_element);
1450    }
1451
1452    if (!G_random_element.inCoeffDomain())
1453      d0= totaldegree (G_random_element, Variable(2),
1454                       Variable (G_random_element.level()));
1455    else
1456      d0= 0;
1457
1458    if (d0 == 0)
1459    {
1460      coF= N (ppA*(cA/gcdcAcB));
1461      coG= N (ppB*(cB/gcdcAcB));
1462      return N(gcdcAcB);
1463    }
1464
1465    if (d0 >  d)
1466    {
1467      if (!find (l, random_element))
1468        l.append (random_element);
1469      continue;
1470    }
1471
1472    G_random_element= (gcdlcAlcB(random_element,x)/uni_lcoeff(G_random_element))
1473                       *G_random_element;
1474
1475    coF_random_element= (lcA(random_element,x)/uni_lcoeff(coF_random_element))
1476                        *coF_random_element;
1477    coG_random_element= (lcB(random_element,x)/uni_lcoeff(coG_random_element))
1478                        *coG_random_element;
1479
1480    if (!G_random_element.inCoeffDomain())
1481      d0= totaldegree (G_random_element, Variable(2),
1482                       Variable (G_random_element.level()));
1483    else
1484      d0= 0;
1485
1486    if (d0 <  d)
1487    {
1488      m= gcdlcAlcB;
1489      newtonPoly= 1;
1490      G_m= 0;
1491      d= d0;
1492      coF_m= 0;
1493      coG_m= 0;
1494    }
1495
1496    TIMING_START (newton_interpolation);
1497    H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
1498    coF= newtonInterp (random_element, coF_random_element, newtonPoly, coF_m,x);
1499    coG= newtonInterp (random_element, coG_random_element, newtonPoly, coG_m,x);
1500    TIMING_END_AND_PRINT (newton_interpolation,
1501                          "time for newton_interpolation: ");
1502
1503    //termination test
1504    if ((uni_lcoeff (H) == gcdlcAlcB) || (G_m == H))
1505    {
1506      TIMING_START (termination_test);
1507      if (gcdlcAlcB.isOne())
1508        cH= 1;
1509      else
1510        cH= uni_content (H);
1511      ppH= H/cH;
1512      ppH /= Lc (ppH);
1513      CanonicalForm lcppH= gcdlcAlcB/cH;
1514      CanonicalForm ccoF= lcppH/Lc (lcppH);
1515      CanonicalForm ccoG= lcppH/Lc (lcppH);
1516      ppCoF= coF/ccoF;
1517      ppCoG= coG/ccoG;
1518      DEBOUTLN (cerr, "ppH= " << ppH);
1519      if (((degree (ppCoF,1)+degree (ppH,1) == bound1) &&
1520           (degree (ppCoG,1)+degree (ppH,1) == bound2) &&
1521           terminationTest (ppA, ppB, ppCoF, ppCoG, ppH)) ||
1522           (fdivides (ppH, ppA, ppCoF) && fdivides (ppH, ppB, ppCoG)))
1523      {
1524        coF= N ((cA/gcdcAcB)*ppCoF);
1525        coG= N ((cB/gcdcAcB)*ppCoG);
1526        TIMING_END_AND_PRINT (termination_test,
1527                              "time for successful termination Fp: ");
1528        return N(gcdcAcB*ppH);
1529      }
1530      TIMING_END_AND_PRINT (termination_test,
1531                            "time for unsuccessful termination Fp: ");
1532    }
1533
1534    G_m= H;
1535    coF_m= coF;
1536    coG_m= coG;
1537    newtonPoly= newtonPoly*(x - random_element);
1538    m= m*(x - random_element);
1539    if (!find (l, random_element))
1540      l.append (random_element);
1541  } while (1);
1542}
1543
1544CFArray
1545solveVandermonde (const CFArray& M, const CFArray& A)
1546{
1547  int r= M.size();
1548  ASSERT (A.size() == r, "vector does not have right size");
1549
1550  if (r == 1)
1551  {
1552    CFArray result= CFArray (1);
1553    result [0]= A [0] / M [0];
1554    return result;
1555  }
1556  // check solvability
1557  bool notDistinct= false;
1558  for (int i= 0; i < r - 1; i++)
1559  {
1560    for (int j= i + 1; j < r; j++)
1561    {
1562      if (M [i] == M [j])
1563      {
1564        notDistinct= true;
1565        break;
1566      }
1567    }
1568  }
1569  if (notDistinct)
1570    return CFArray();
1571
1572  CanonicalForm master= 1;
1573  Variable x= Variable (1);
1574  for (int i= 0; i < r; i++)
1575    master *= x - M [i];
1576  CFList Pj;
1577  CanonicalForm tmp;
1578  for (int i= 0; i < r; i++)
1579  {
1580    tmp= master/(x - M [i]);
1581    tmp /= tmp (M [i], 1);
1582    Pj.append (tmp);
1583  }
1584  CFArray result= CFArray (r);
1585
1586  CFListIterator j= Pj;
1587  for (int i= 1; i <= r; i++, j++)
1588  {
1589    tmp= 0;
1590    for (int l= 0; l < A.size(); l++)
1591      tmp += A[l]*j.getItem()[l];
1592    result[i - 1]= tmp;
1593  }
1594  return result;
1595}
1596
1597CFArray
1598solveGeneralVandermonde (const CFArray& M, const CFArray& A)
1599{
1600  int r= M.size();
1601  ASSERT (A.size() == r, "vector does not have right size");
1602  if (r == 1)
1603  {
1604    CFArray result= CFArray (1);
1605    result [0]= A[0] / M [0];
1606    return result;
1607  }
1608  // check solvability
1609  bool notDistinct= false;
1610  for (int i= 0; i < r - 1; i++)
1611  {
1612    for (int j= i + 1; j < r; j++)
1613    {
1614      if (M [i] == M [j])
1615      {
1616        notDistinct= true;
1617        break;
1618      }
1619    }
1620  }
1621  if (notDistinct)
1622    return CFArray();
1623
1624  CanonicalForm master= 1;
1625  Variable x= Variable (1);
1626  for (int i= 0; i < r; i++)
1627    master *= x - M [i];
1628  master *= x;
1629  CFList Pj;
1630  CanonicalForm tmp;
1631  for (int i= 0; i < r; i++)
1632  {
1633    tmp= master/(x - M [i]);
1634    tmp /= tmp (M [i], 1);
1635    Pj.append (tmp);
1636  }
1637
1638  CFArray result= CFArray (r);
1639
1640  CFListIterator j= Pj;
1641  for (int i= 1; i <= r; i++, j++)
1642  {
1643    tmp= 0;
1644
1645    for (int l= 1; l <= A.size(); l++)
1646      tmp += A[l - 1]*j.getItem()[l];
1647    result[i - 1]= tmp;
1648  }
1649  return result;
1650}
1651
1652/// M in row echolon form, rk rank of M
1653CFArray
1654readOffSolution (const CFMatrix& M, const long rk)
1655{
1656  CFArray result= CFArray (rk);
1657  CanonicalForm tmp1, tmp2, tmp3;
1658  for (int i= rk; i >= 1; i--)
1659  {
1660    tmp3= 0;
1661    tmp1= M (i, M.columns());
1662    for (int j= M.columns() - 1; j >= 1; j--)
1663    {
1664      tmp2= M (i, j);
1665      if (j == i)
1666        break;
1667      else
1668        tmp3 += tmp2*result[j - 1];
1669    }
1670    result[i - 1]= (tmp1 - tmp3)/tmp2;
1671  }
1672  return result;
1673}
1674
1675CFArray
1676readOffSolution (const CFMatrix& M, const CFArray& L, const CFArray& partialSol)
1677{
1678  CFArray result= CFArray (M.rows());
1679  CanonicalForm tmp1, tmp2, tmp3;
1680  int k;
1681  for (int i= M.rows(); i >= 1; i--)
1682  {
1683    tmp3= 0;
1684    tmp1= L[i - 1];
1685    k= 0;
1686    for (int j= M.columns(); j >= 1; j--, k++)
1687    {
1688      tmp2= M (i, j);
1689      if (j == i)
1690        break;
1691      else
1692      {
1693        if (k > partialSol.size() - 1)
1694          tmp3 += tmp2*result[j - 1];
1695        else
1696          tmp3 += tmp2*partialSol[partialSol.size() - k - 1];
1697      }
1698    }
1699    result [i - 1]= (tmp1 - tmp3)/tmp2;
1700  }
1701  return result;
1702}
1703
1704long
1705gaussianElimFp (CFMatrix& M, CFArray& L)
1706{
1707  ASSERT (L.size() <= M.rows(), "dimension exceeded");
1708  CFMatrix *N;
1709  N= new CFMatrix (M.rows(), M.columns() + 1);
1710
1711  for (int i= 1; i <= M.rows(); i++)
1712    for (int j= 1; j <= M.columns(); j++)
1713      (*N) (i, j)= M (i, j);
1714
1715  int j= 1;
1716  for (int i= 0; i < L.size(); i++, j++)
1717    (*N) (j, M.columns() + 1)= L[i];
1718#ifdef HAVE_FLINT
1719  nmod_mat_t FLINTN;
1720  convertFacCFMatrix2nmod_mat_t (FLINTN, *N);
1721  long rk= nmod_mat_rref (FLINTN);
1722
1723  N= convertNmod_mat_t2FacCFMatrix (FLINTN);
1724  nmod_mat_clear (FLINTN);
1725#else
1726  int p= getCharacteristic ();
1727  if (fac_NTL_char != p)
1728  {
1729    fac_NTL_char= p;
1730    zz_p::init (p);
1731  }
1732  mat_zz_p *NTLN= convertFacCFMatrix2NTLmat_zz_p(*N);
1733  long rk= gauss (*NTLN);
1734
1735  N= convertNTLmat_zz_p2FacCFMatrix (*NTLN);
1736#endif
1737
1738  L= CFArray (M.rows());
1739  for (int i= 0; i < M.rows(); i++)
1740    L[i]= (*N) (i + 1, M.columns() + 1);
1741  M= (*N) (1, M.rows(), 1, M.columns());
1742  delete N;
1743  return rk;
1744}
1745
1746long
1747gaussianElimFq (CFMatrix& M, CFArray& L, const Variable& alpha)
1748{
1749  ASSERT (L.size() <= M.rows(), "dimension exceeded");
1750  CFMatrix *N;
1751  N= new CFMatrix (M.rows(), M.columns() + 1);
1752
1753  for (int i= 1; i <= M.rows(); i++)
1754    for (int j= 1; j <= M.columns(); j++)
1755      (*N) (i, j)= M (i, j);
1756
1757  int j= 1;
1758  for (int i= 0; i < L.size(); i++, j++)
1759    (*N) (j, M.columns() + 1)= L[i];
1760  int p= getCharacteristic ();
1761  if (fac_NTL_char != p)
1762  {
1763    fac_NTL_char= p;
1764    zz_p::init (p);
1765  }
1766  zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
1767  zz_pE::init (NTLMipo);
1768  mat_zz_pE *NTLN= convertFacCFMatrix2NTLmat_zz_pE(*N);
1769  long rk= gauss (*NTLN);
1770
1771  N= convertNTLmat_zz_pE2FacCFMatrix (*NTLN, alpha);
1772
1773  M= (*N) (1, M.rows(), 1, M.columns());
1774  L= CFArray (M.rows());
1775  for (int i= 0; i < M.rows(); i++)
1776    L[i]= (*N) (i + 1, M.columns() + 1);
1777
1778  delete N;
1779  return rk;
1780}
1781
1782CFArray
1783solveSystemFp (const CFMatrix& M, const CFArray& L)
1784{
1785  ASSERT (L.size() <= M.rows(), "dimension exceeded");
1786  CFMatrix *N;
1787  N= new CFMatrix (M.rows(), M.columns() + 1);
1788
1789  for (int i= 1; i <= M.rows(); i++)
1790    for (int j= 1; j <= M.columns(); j++)
1791      (*N) (i, j)= M (i, j);
1792
1793  int j= 1;
1794  for (int i= 0; i < L.size(); i++, j++)
1795    (*N) (j, M.columns() + 1)= L[i];
1796
1797#ifdef HAVE_FLINT
1798  nmod_mat_t FLINTN;
1799  convertFacCFMatrix2nmod_mat_t (FLINTN, *N);
1800  long rk= nmod_mat_rref (FLINTN);
1801#else
1802  int p= getCharacteristic ();
1803  if (fac_NTL_char != p)
1804  {
1805    fac_NTL_char= p;
1806    zz_p::init (p);
1807  }
1808  mat_zz_p *NTLN= convertFacCFMatrix2NTLmat_zz_p(*N);
1809  long rk= gauss (*NTLN);
1810#endif
1811  if (rk != M.columns())
1812  {
1813#ifdef HAVE_FLINT
1814    nmod_mat_clear (FLINTN);
1815#endif
1816    delete N;
1817    return CFArray();
1818  }
1819#ifdef HAVE_FLINT
1820  N= convertNmod_mat_t2FacCFMatrix (FLINTN);
1821  nmod_mat_clear (FLINTN);
1822#else
1823  N= convertNTLmat_zz_p2FacCFMatrix (*NTLN);
1824#endif
1825  CFArray A= readOffSolution (*N, rk);
1826
1827  delete N;
1828  return A;
1829}
1830
1831CFArray
1832solveSystemFq (const CFMatrix& M, const CFArray& L, const Variable& alpha)
1833{
1834  ASSERT (L.size() <= M.rows(), "dimension exceeded");
1835  CFMatrix *N;
1836  N= new CFMatrix (M.rows(), M.columns() + 1);
1837
1838  for (int i= 1; i <= M.rows(); i++)
1839    for (int j= 1; j <= M.columns(); j++)
1840      (*N) (i, j)= M (i, j);
1841  int j= 1;
1842  for (int i= 0; i < L.size(); i++, j++)
1843    (*N) (j, M.columns() + 1)= L[i];
1844  int p= getCharacteristic ();
1845  if (fac_NTL_char != p)
1846  {
1847    fac_NTL_char= p;
1848    zz_p::init (p);
1849  }
1850  zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
1851  zz_pE::init (NTLMipo);
1852  mat_zz_pE *NTLN= convertFacCFMatrix2NTLmat_zz_pE(*N);
1853  long rk= gauss (*NTLN);
1854  if (rk != M.columns())
1855  {
1856    delete N;
1857    return CFArray();
1858  }
1859  N= convertNTLmat_zz_pE2FacCFMatrix (*NTLN, alpha);
1860
1861  CFArray A= readOffSolution (*N, rk);
1862
1863  delete N;
1864  return A;
1865}
1866#endif
1867
1868CFArray
1869getMonoms (const CanonicalForm& F)
1870{
1871  if (F.inCoeffDomain())
1872  {
1873    CFArray result= CFArray (1);
1874    result [0]= 1;
1875    return result;
1876  }
1877  if (F.isUnivariate())
1878  {
1879    CFArray result= CFArray (size(F));
1880    int j= 0;
1881    for (CFIterator i= F; i.hasTerms(); i++, j++)
1882      result[j]= power (F.mvar(), i.exp());
1883    return result;
1884  }
1885  int numMon= size (F);
1886  CFArray result= CFArray (numMon);
1887  int j= 0;
1888  CFArray recResult;
1889  Variable x= F.mvar();
1890  CanonicalForm powX;
1891  for (CFIterator i= F; i.hasTerms(); i++)
1892  {
1893    powX= power (x, i.exp());
1894    recResult= getMonoms (i.coeff());
1895    for (int k= 0; k < recResult.size(); k++)
1896      result[j+k]= powX*recResult[k];
1897    j += recResult.size();
1898  }
1899  return result;
1900}
1901
1902#ifdef HAVE_NTL
1903CFArray
1904evaluateMonom (const CanonicalForm& F, const CFList& evalPoints)
1905{
1906  if (F.inCoeffDomain())
1907  {
1908    CFArray result= CFArray (1);
1909    result [0]= F;
1910    return result;
1911  }
1912  if (F.isUnivariate())
1913  {
1914    ASSERT (evalPoints.length() == 1,
1915            "expected an eval point with only one component");
1916    CFArray result= CFArray (size(F));
1917    int j= 0;
1918    CanonicalForm evalPoint= evalPoints.getLast();
1919    for (CFIterator i= F; i.hasTerms(); i++, j++)
1920      result[j]= power (evalPoint, i.exp());
1921    return result;
1922  }
1923  int numMon= size (F);
1924  CFArray result= CFArray (numMon);
1925  int j= 0;
1926  CanonicalForm evalPoint= evalPoints.getLast();
1927  CFList buf= evalPoints;
1928  buf.removeLast();
1929  CFArray recResult;
1930  CanonicalForm powEvalPoint;
1931  for (CFIterator i= F; i.hasTerms(); i++)
1932  {
1933    powEvalPoint= power (evalPoint, i.exp());
1934    recResult= evaluateMonom (i.coeff(), buf);
1935    for (int k= 0; k < recResult.size(); k++)
1936      result[j+k]= powEvalPoint*recResult[k];
1937    j += recResult.size();
1938  }
1939  return result;
1940}
1941
1942CFArray
1943evaluate (const CFArray& A, const CFList& evalPoints)
1944{
1945  CFArray result= A.size();
1946  CanonicalForm tmp;
1947  int k;
1948  for (int i= 0; i < A.size(); i++)
1949  {
1950    tmp= A[i];
1951    k= 1;
1952    for (CFListIterator j= evalPoints; j.hasItem(); j++, k++)
1953      tmp= tmp (j.getItem(), k);
1954    result[i]= tmp;
1955  }
1956  return result;
1957}
1958
1959CFList
1960evaluationPoints (const CanonicalForm& F, const CanonicalForm& G,
1961                  CanonicalForm& Feval, CanonicalForm& Geval,
1962                  const CanonicalForm& LCF, const bool& GF,
1963                  const Variable& alpha, bool& fail, CFList& list
1964                 )
1965{
1966  int k= tmax (F.level(), G.level()) - 1;
1967  Variable x= Variable (1);
1968  CFList result;
1969  FFRandom genFF;
1970  GFRandom genGF;
1971  int p= getCharacteristic ();
1972  double bound;
1973  if (alpha != Variable (1))
1974  {
1975    bound= pow ((double) p, (double) degree (getMipo(alpha)));
1976    bound= pow (bound, (double) k);
1977  }
1978  else if (GF)
1979  {
1980    bound= pow ((double) p, (double) getGFDegree());
1981    bound= pow ((double) bound, (double) k);
1982  }
1983  else
1984    bound= pow ((double) p, (double) k);
1985
1986  CanonicalForm random;
1987  int j;
1988  bool zeroOneOccured= false;
1989  bool allEqual= false;
1990  CanonicalForm buf;
1991  do
1992  {
1993    random= 0;
1994    // possible overflow if list.length() does not fit into a int
1995    if (list.length() >= bound)
1996    {
1997      fail= true;
1998      break;
1999    }
2000    for (int i= 0; i < k; i++)
2001    {
2002      if (GF)
2003      {
2004        result.append (genGF.generate());
2005        random += result.getLast()*power (x, i);
2006      }
2007      else if (alpha.level() != 1)
2008      {
2009        AlgExtRandomF genAlgExt (alpha);
2010        result.append (genAlgExt.generate());
2011        random += result.getLast()*power (x, i);
2012      }
2013      else
2014      {
2015        result.append (genFF.generate());
2016        random += result.getLast()*power (x, i);
2017      }
2018      if (result.getLast().isOne() || result.getLast().isZero())
2019        zeroOneOccured= true;
2020    }
2021    if (find (list, random))
2022    {
2023      zeroOneOccured= false;
2024      allEqual= false;
2025      result= CFList();
2026      continue;
2027    }
2028    if (zeroOneOccured)
2029    {
2030      list.append (random);
2031      zeroOneOccured= false;
2032      allEqual= false;
2033      result= CFList();
2034      continue;
2035    }
2036    // no zero at this point
2037    if (k > 1)
2038    {
2039      allEqual= true;
2040      CFIterator iter= random;
2041      buf= iter.coeff();
2042      iter++;
2043      for (; iter.hasTerms(); iter++)
2044        if (buf != iter.coeff())
2045          allEqual= false;
2046    }
2047    if (allEqual)
2048    {
2049      list.append (random);
2050      allEqual= false;
2051      zeroOneOccured= false;
2052      result= CFList();
2053      continue;
2054    }
2055
2056    Feval= F;
2057    Geval= G;
2058    CanonicalForm LCeval= LCF;
2059    j= 1;
2060    for (CFListIterator i= result; i.hasItem(); i++, j++)
2061    {
2062      Feval= Feval (i.getItem(), j);
2063      Geval= Geval (i.getItem(), j);
2064      LCeval= LCeval (i.getItem(), j);
2065    }
2066
2067    if (LCeval.isZero())
2068    {
2069      if (!find (list, random))
2070        list.append (random);
2071      zeroOneOccured= false;
2072      allEqual= false;
2073      result= CFList();
2074      continue;
2075    }
2076
2077    if (list.length() >= bound)
2078    {
2079      fail= true;
2080      break;
2081    }
2082  } while (find (list, random));
2083
2084  return result;
2085}
2086
2087/// multiply two lists componentwise
2088void mult (CFList& L1, const CFList& L2)
2089{
2090  ASSERT (L1.length() == L2.length(), "lists of the same size expected");
2091
2092  CFListIterator j= L2;
2093  for (CFListIterator i= L1; i.hasItem(); i++, j++)
2094    i.getItem() *= j.getItem();
2095}
2096
2097void eval (const CanonicalForm& A, const CanonicalForm& B, CanonicalForm& Aeval,
2098           CanonicalForm& Beval, const CFList& L)
2099{
2100  Aeval= A;
2101  Beval= B;
2102  int j= 1;
2103  for (CFListIterator i= L; i.hasItem(); i++, j++)
2104  {
2105    Aeval= Aeval (i.getItem(), j);
2106    Beval= Beval (i.getItem(), j);
2107  }
2108}
2109
2110CanonicalForm
2111monicSparseInterpol (const CanonicalForm& F, const CanonicalForm& G,
2112                     const CanonicalForm& skeleton, const Variable& alpha,
2113                     bool& fail, CFArray*& coeffMonoms, CFArray& Monoms
2114                    )
2115{
2116  CanonicalForm A= F;
2117  CanonicalForm B= G;
2118  if (F.isZero() && degree(G) > 0) return G/Lc(G);
2119  else if (G.isZero() && degree (F) > 0) return F/Lc(F);
2120  else if (F.isZero() && G.isZero()) return F.genOne();
2121  if (F.inBaseDomain() || G.inBaseDomain()) return F.genOne();
2122  if (F.isUnivariate() && fdivides(F, G)) return F/Lc(F);
2123  if (G.isUnivariate() && fdivides(G, F)) return G/Lc(G);
2124  if (F == G) return F/Lc(F);
2125
2126  ASSERT (degree (A, 1) == 0, "expected degree (F, 1) == 0");
2127  ASSERT (degree (B, 1) == 0, "expected degree (G, 1) == 0");
2128
2129  CFMap M,N;
2130  int best_level= myCompress (A, B, M, N, false);
2131
2132  if (best_level == 0)
2133    return B.genOne();
2134
2135  A= M(A);
2136  B= M(B);
2137
2138  Variable x= Variable (1);
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
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  delete[] coeffMonoms;
2352
2353  if (alpha.level() != 1 && V_buf != alpha)
2354  {
2355    CFList u, v;
2356    result= mapDown (result, prim_elem, im_prim_elem, alpha, u, v);
2357  }
2358
2359  result= N(result);
2360  if (fdivides (result, F) && fdivides (result, G))
2361    return result;
2362  else
2363  {
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  ASSERT (degree (A, 1) == 0, "expected degree (F, 1) == 0");
2386  ASSERT (degree (B, 1) == 0, "expected degree (G, 1) == 0");
2387
2388  CFMap M,N;
2389  int best_level= myCompress (A, B, M, N, false);
2390
2391  if (best_level == 0)
2392    return B.genOne();
2393
2394  A= M(A);
2395  B= M(B);
2396
2397  Variable x= Variable (1);
2398
2399  //univariate case
2400  if (A.isUnivariate() && B.isUnivariate())
2401    return N (gcd (A, B));
2402
2403  CanonicalForm skel= M(skeleton);
2404
2405  CanonicalForm cA, cB;    // content of A and B
2406  CanonicalForm ppA, ppB;    // primitive part of A and B
2407  CanonicalForm gcdcAcB;
2408  cA = uni_content (A);
2409  cB = uni_content (B);
2410  gcdcAcB= gcd (cA, cB);
2411  ppA= A/cA;
2412  ppB= B/cB;
2413
2414  CanonicalForm lcA, lcB;  // leading coefficients of A and B
2415  CanonicalForm gcdlcAlcB;
2416  lcA= uni_lcoeff (ppA);
2417  lcB= uni_lcoeff (ppB);
2418
2419  if (fdivides (lcA, lcB))
2420  {
2421    if (fdivides (A, B))
2422      return F/Lc(F);
2423  }
2424  if (fdivides (lcB, lcA))
2425  {
2426    if (fdivides (B, A))
2427      return G/Lc(G);
2428  }
2429
2430  gcdlcAlcB= gcd (lcA, lcB);
2431  int skelSize= size (skel, skel.mvar());
2432
2433  int j= 0;
2434  int biggestSize= 0;
2435  int bufSize;
2436  int numberUni= 0;
2437  for (CFIterator i= skel; i.hasTerms(); i++, j++)
2438  {
2439    bufSize= size (i.coeff());
2440    biggestSize= tmax (biggestSize, bufSize);
2441    numberUni += bufSize;
2442  }
2443  numberUni--;
2444  numberUni= (int) ceil ( (double) numberUni / (skelSize - 1));
2445  biggestSize= tmax (biggestSize , numberUni);
2446
2447  numberUni= biggestSize + size (LC(skel)) - 1;
2448  int biggestSize2= tmax (numberUni, biggestSize);
2449
2450  CanonicalForm g, Aeval, Beval;
2451
2452  CFList evalPoints;
2453  CFArray coeffEval;
2454  bool evalFail= false;
2455  CFList list;
2456  bool GF= false;
2457  CanonicalForm LCA= LC (A);
2458  CanonicalForm tmp;
2459  CFArray gcds= CFArray (biggestSize);
2460  CFList * pEvalPoints= new CFList [biggestSize];
2461  Variable V_buf= alpha;
2462  CFList source, dest;
2463  CanonicalForm prim_elem, im_prim_elem;
2464  for (int i= 0; i < biggestSize; i++)
2465  {
2466    if (i == 0)
2467    {
2468      if (getCharacteristic() > 3)
2469        evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf,
2470                                    evalFail, list);
2471      else
2472        evalFail= true;
2473
2474      if (evalFail)
2475      {
2476        if (V_buf.level() != 1)
2477        {
2478          do
2479          {
2480            Variable V_buf2= chooseExtension (V_buf);
2481            source= CFList();
2482            dest= CFList();
2483
2484            bool prim_fail= false;
2485            Variable V_buf3;
2486            prim_elem= primitiveElement (V_buf, V_buf3, prim_fail);
2487
2488            ASSERT (!prim_fail, "failure in integer factorizer");
2489            if (prim_fail)
2490              ; //ERROR
2491            else
2492              im_prim_elem= mapPrimElem (prim_elem, V_buf, V_buf2);
2493
2494            DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf));
2495            DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf2));
2496
2497            for (CFListIterator i= list; i.hasItem(); i++)
2498              i.getItem()= mapUp (i.getItem(), V_buf, V_buf2, prim_elem,
2499                                im_prim_elem, source, dest);
2500            evalFail= false;
2501            evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf,
2502                                          evalFail, list);
2503          } while (evalFail);
2504        }
2505        else
2506        {
2507          CanonicalForm mipo;
2508          int deg= 2;
2509          do {
2510            mipo= randomIrredpoly (deg, x);
2511            V_buf= rootOf (mipo);
2512            evalFail= false;
2513            evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf,
2514                                          evalFail, list);
2515            deg++;
2516          } while (evalFail);
2517        }
2518      }
2519    }
2520    else
2521    {
2522      mult (evalPoints, pEvalPoints[0]);
2523      eval (A, B, Aeval, Beval, evalPoints);
2524    }
2525
2526    g= gcd (Aeval, Beval);
2527    g /= Lc (g);
2528
2529    if (degree (g) != degree (skel) || (skelSize != size (g)))
2530    {
2531      delete[] pEvalPoints;
2532      fail= true;
2533      return 0;
2534    }
2535    CFIterator l= skel;
2536    for (CFIterator k= g; k.hasTerms(); k++, l++)
2537    {
2538      if (k.exp() != l.exp())
2539      {
2540        delete[] pEvalPoints;
2541        fail= true;
2542        return 0;
2543      }
2544    }
2545    pEvalPoints[i]= evalPoints;
2546    gcds[i]= g;
2547
2548    tmp= 0;
2549    int j= 0;
2550    for (CFListIterator k= evalPoints; k.hasItem(); k++, j++)
2551      tmp += k.getItem()*power (x, j);
2552    list.append (tmp);
2553  }
2554
2555  if (Monoms.size() == 0)
2556    Monoms= getMonoms (skel);
2557
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    delete[] pEvalPoints;
2870    delete[] pMat;
2871    delete[] pL;
2872    delete[] coeffMonoms;
2873    delete[] pM;
2874
2875    if (bufpEvalPoints != NULL)
2876      delete [] bufpEvalPoints;
2877    if (fdivides (result, F) && fdivides (result, G))
2878      return result;
2879    else
2880    {
2881      fail= true;
2882      return 0;
2883    }
2884  } // end of deKleine, Monagan & Wittkopf
2885
2886  result += Monoms[0];
2887  int ind2= 0, ind3= 2;
2888  ind= 0;
2889  for (int l= 0; l < minimalColumnsIndex; l++)
2890    ind += coeffMonoms[l].size();
2891  for (int l= solution.size() - pMat[0].columns() + 1; l < solution.size();
2892       l++, ind2++, ind3++)
2893  {
2894    result += solution[l]*Monoms [1 + ind2];
2895    for (int i= 0; i < pMat[0].rows(); i++)
2896      firstColumn[i] += solution[l]*pMat[0] (i + 1, ind3);
2897  }
2898  for (int l= 0; l < solution.size() + 1 - pMat[0].columns(); l++)
2899    result += solution[l]*Monoms [ind + l];
2900  ind= coeffMonoms[0].size();
2901  for (int k= 1; k < skelSize; k++)
2902  {
2903    if (k == minimalColumnsIndex)
2904    {
2905      ind += coeffMonoms[k].size();
2906      continue;
2907    }
2908    if (k != minimalColumnsIndex)
2909    {
2910      for (int i= 0; i < biggestSize; i++)
2911        pL[k] [i] *= firstColumn [i];
2912    }
2913
2914    if (biggestSize != coeffMonoms[k].size() && k != minimalColumnsIndex)
2915    {
2916      bufArray= CFArray (coeffMonoms[k].size());
2917      for (int i= 0; i < bufArray.size(); i++)
2918        bufArray [i]= pL[k] [i];
2919      solution= solveGeneralVandermonde (pM[k], bufArray);
2920    }
2921    else
2922      solution= solveGeneralVandermonde (pM[k], pL[k]);
2923
2924    if (solution.size() == 0)
2925    {
2926      delete[] pEvalPoints;
2927      delete[] pMat;
2928      delete[] pL;
2929      delete[] coeffMonoms;
2930      delete[] pM;
2931      fail= true;
2932      return 0;
2933    }
2934    if (k != minimalColumnsIndex)
2935    {
2936      for (int l= 0; l < solution.size(); l++)
2937        result += solution[l]*Monoms [ind + l];
2938      ind += solution.size();
2939    }
2940  }
2941
2942  delete[] pEvalPoints;
2943  delete[] pMat;
2944  delete[] pL;
2945  delete[] pM;
2946  delete[] coeffMonoms;
2947
2948  if (alpha.level() != 1 && V_buf != alpha)
2949  {
2950    CFList u, v;
2951    result= mapDown (result, prim_elem, im_prim_elem, alpha, u, v);
2952  }
2953  result= N(result);
2954
2955  if (fdivides (result, F) && fdivides (result, G))
2956    return result;
2957  else
2958  {
2959    fail= true;
2960    return 0;
2961  }
2962}
2963
2964CanonicalForm sparseGCDFq (const CanonicalForm& F, const CanonicalForm& G,
2965                           const Variable & alpha, CFList& l, bool& topLevel)
2966{
2967  CanonicalForm A= F;
2968  CanonicalForm B= G;
2969  if (F.isZero() && degree(G) > 0) return G/Lc(G);
2970  else if (G.isZero() && degree (F) > 0) return F/Lc(F);
2971  else if (F.isZero() && G.isZero()) return F.genOne();
2972  if (F.inBaseDomain() || G.inBaseDomain()) return F.genOne();
2973  if (F.isUnivariate() && fdivides(F, G)) return F/Lc(F);
2974  if (G.isUnivariate() && fdivides(G, F)) return G/Lc(G);
2975  if (F == G) return F/Lc(F);
2976
2977  CFMap M,N;
2978  int best_level= myCompress (A, B, M, N, topLevel);
2979
2980  if (best_level == 0) return B.genOne();
2981
2982  A= M(A);
2983  B= M(B);
2984
2985  Variable x= Variable (1);
2986
2987  //univariate case
2988  if (A.isUnivariate() && B.isUnivariate())
2989    return N (gcd (A, B));
2990
2991  CanonicalForm cA, cB;    // content of A and B
2992  CanonicalForm ppA, ppB;    // primitive part of A and B
2993  CanonicalForm gcdcAcB;
2994
2995  cA = uni_content (A);
2996  cB = uni_content (B);
2997  gcdcAcB= gcd (cA, cB);
2998  ppA= A/cA;
2999  ppB= B/cB;
3000
3001  CanonicalForm lcA, lcB;  // leading coefficients of A and B
3002  CanonicalForm gcdlcAlcB;
3003  lcA= uni_lcoeff (ppA);
3004  lcB= uni_lcoeff (ppB);
3005
3006  if (fdivides (lcA, lcB))
3007  {
3008    if (fdivides (A, B))
3009      return F/Lc(F);
3010  }
3011  if (fdivides (lcB, lcA))
3012  {
3013    if (fdivides (B, A))
3014      return G/Lc(G);
3015  }
3016
3017  gcdlcAlcB= gcd (lcA, lcB);
3018
3019  int d= totaldegree (ppA, Variable (2), Variable (ppA.level()));
3020  int d0;
3021
3022  if (d == 0)
3023    return N(gcdcAcB);
3024  d0= totaldegree (ppB, Variable (2), Variable (ppB.level()));
3025
3026  if (d0 < d)
3027    d= d0;
3028
3029  if (d == 0)
3030    return N(gcdcAcB);
3031
3032  CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH, skeleton;
3033  CanonicalForm newtonPoly= 1;
3034  m= gcdlcAlcB;
3035  G_m= 0;
3036  H= 0;
3037  bool fail= false;
3038  topLevel= false;
3039  bool inextension= false;
3040  Variable V_buf= alpha;
3041  CanonicalForm prim_elem, im_prim_elem;
3042  CFList source, dest;
3043  do // first do
3044  {
3045    random_element= randomElement (m, V_buf, l, fail);
3046    if (random_element == 0 && !fail)
3047    {
3048      if (!find (l, random_element))
3049        l.append (random_element);
3050      continue;
3051    }
3052    if (fail)
3053    {
3054      source= CFList();
3055      dest= CFList();
3056
3057      Variable V_buf3= V_buf;
3058      V_buf= chooseExtension (V_buf);
3059      bool prim_fail= false;
3060      Variable V_buf2;
3061      prim_elem= primitiveElement (alpha, V_buf2, prim_fail);
3062
3063      if (V_buf3 != alpha)
3064      {
3065        m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
3066        G_m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
3067        newtonPoly= mapDown (newtonPoly, prim_elem, im_prim_elem, alpha,
3068                             source, dest);
3069        ppA= mapDown (ppA, prim_elem, im_prim_elem, alpha, source, dest);
3070        ppB= mapDown (ppB, prim_elem, im_prim_elem, alpha, source, dest);
3071        gcdlcAlcB= mapDown (gcdlcAlcB, prim_elem, im_prim_elem, alpha, source,
3072                            dest);
3073        for (CFListIterator i= l; i.hasItem(); i++)
3074          i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, alpha,
3075                                source, dest);
3076      }
3077
3078      ASSERT (!prim_fail, "failure in integer factorizer");
3079      if (prim_fail)
3080        ; //ERROR
3081      else
3082        im_prim_elem= mapPrimElem (prim_elem, alpha, V_buf);
3083
3084      DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (alpha));
3085      DEBOUTLN (cerr, "getMipo (V_buf2)= " << getMipo (V_buf2));
3086      inextension= true;
3087      for (CFListIterator i= l; i.hasItem(); i++)
3088        i.getItem()= mapUp (i.getItem(), alpha, V_buf, prim_elem,
3089                             im_prim_elem, source, dest);
3090      m= mapUp (m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3091      G_m= mapUp (G_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3092      newtonPoly= mapUp (newtonPoly, alpha, V_buf, prim_elem, im_prim_elem,
3093                          source, dest);
3094      ppA= mapUp (ppA, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3095      ppB= mapUp (ppB, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3096      gcdlcAlcB= mapUp (gcdlcAlcB, alpha, V_buf, prim_elem, im_prim_elem,
3097                         source, dest);
3098
3099      fail= false;
3100      random_element= randomElement (m, V_buf, l, fail );
3101      DEBOUTLN (cerr, "fail= " << fail);
3102      CFList list;
3103      TIMING_START (gcd_recursion);
3104      G_random_element=
3105      sparseGCDFq (ppA (random_element, x), ppB (random_element, x), V_buf,
3106                        list, topLevel);
3107      TIMING_END_AND_PRINT (gcd_recursion,
3108                            "time for recursive call: ");
3109      DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3110    }
3111    else
3112    {
3113      CFList list;
3114      TIMING_START (gcd_recursion);
3115      G_random_element=
3116      sparseGCDFq (ppA(random_element, x), ppB(random_element, x), V_buf,
3117                        list, topLevel);
3118      TIMING_END_AND_PRINT (gcd_recursion,
3119                            "time for recursive call: ");
3120      DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3121    }
3122
3123    if (!G_random_element.inCoeffDomain())
3124      d0= totaldegree (G_random_element, Variable(2),
3125                       Variable (G_random_element.level()));
3126    else
3127      d0= 0;
3128
3129    if (d0 == 0)
3130      return N(gcdcAcB);
3131    if (d0 >  d)
3132    {
3133      if (!find (l, random_element))
3134        l.append (random_element);
3135      continue;
3136    }
3137
3138    G_random_element=
3139    (gcdlcAlcB(random_element, x)/uni_lcoeff (G_random_element))
3140    * G_random_element;
3141
3142    skeleton= G_random_element;
3143    if (!G_random_element.inCoeffDomain())
3144      d0= totaldegree (G_random_element, Variable(2),
3145                       Variable (G_random_element.level()));
3146    else
3147      d0= 0;
3148
3149    if (d0 <  d)
3150    {
3151      m= gcdlcAlcB;
3152      newtonPoly= 1;
3153      G_m= 0;
3154      d= d0;
3155    }
3156
3157    H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
3158    if (uni_lcoeff (H) == gcdlcAlcB)
3159    {
3160      cH= uni_content (H);
3161      ppH= H/cH;
3162      if (inextension)
3163      {
3164        CFList u, v;
3165        //maybe it's better to test if ppH is an element of F(\alpha) before
3166        //mapping down
3167        if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
3168        {
3169          DEBOUTLN (cerr, "ppH before mapDown= " << ppH);
3170          ppH= mapDown (ppH, prim_elem, im_prim_elem, alpha, u, v);
3171          ppH /= Lc(ppH);
3172          DEBOUTLN (cerr, "ppH after mapDown= " << ppH);
3173          return N(gcdcAcB*ppH);
3174        }
3175      }
3176      else if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
3177        return N(gcdcAcB*ppH);
3178    }
3179    G_m= H;
3180    newtonPoly= newtonPoly*(x - random_element);
3181    m= m*(x - random_element);
3182    if (!find (l, random_element))
3183      l.append (random_element);
3184
3185    if (getCharacteristic () > 3 && size (skeleton) < 100)
3186    {
3187      CFArray Monoms;
3188      CFArray *coeffMonoms;
3189      do //second do
3190      {
3191        random_element= randomElement (m, V_buf, l, fail);
3192        if (random_element == 0 && !fail)
3193        {
3194          if (!find (l, random_element))
3195            l.append (random_element);
3196          continue;
3197        }
3198        if (fail)
3199        {
3200          source= CFList();
3201          dest= CFList();
3202
3203          Variable V_buf3= V_buf;
3204          V_buf= chooseExtension (V_buf);
3205          bool prim_fail= false;
3206          Variable V_buf2;
3207          prim_elem= primitiveElement (alpha, V_buf2, prim_fail);
3208
3209          if (V_buf3 != alpha)
3210          {
3211            m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
3212            G_m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
3213            newtonPoly= mapDown (newtonPoly, prim_elem, im_prim_elem, alpha,
3214                                 source, dest);
3215            ppA= mapDown (ppA, prim_elem, im_prim_elem, alpha, source, dest);
3216            ppB= mapDown (ppB, prim_elem, im_prim_elem, alpha, source, dest);
3217            gcdlcAlcB= mapDown (gcdlcAlcB, prim_elem, im_prim_elem, alpha,
3218                                source, dest);
3219            for (CFListIterator i= l; i.hasItem(); i++)
3220              i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, alpha,
3221                                    source, dest);
3222          }
3223
3224          ASSERT (!prim_fail, "failure in integer factorizer");
3225          if (prim_fail)
3226            ; //ERROR
3227          else
3228            im_prim_elem= mapPrimElem (prim_elem, alpha, V_buf);
3229
3230          DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (alpha));
3231          DEBOUTLN (cerr, "getMipo (V_buf2)= " << getMipo (V_buf2));
3232          inextension= true;
3233          for (CFListIterator i= l; i.hasItem(); i++)
3234            i.getItem()= mapUp (i.getItem(), alpha, V_buf, prim_elem,
3235                                im_prim_elem, source, dest);
3236          m= mapUp (m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3237          G_m= mapUp (G_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3238          newtonPoly= mapUp (newtonPoly, alpha, V_buf, prim_elem, im_prim_elem,
3239                              source, dest);
3240          ppA= mapUp (ppA, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3241          ppB= mapUp (ppB, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3242
3243          gcdlcAlcB= mapUp (gcdlcAlcB, alpha, V_buf, prim_elem, im_prim_elem,
3244                            source, dest);
3245
3246          fail= false;
3247          random_element= randomElement (m, V_buf, l, fail);
3248          DEBOUTLN (cerr, "fail= " << fail);
3249          CFList list;
3250          TIMING_START (gcd_recursion);
3251
3252          //sparseInterpolation
3253          bool sparseFail= false;
3254          if (LC (skeleton).inCoeffDomain())
3255            G_random_element=
3256            monicSparseInterpol (ppA (random_element, x), ppB(random_element,x),
3257                                skeleton,V_buf, sparseFail, coeffMonoms,Monoms);
3258          else
3259            G_random_element=
3260            nonMonicSparseInterpol (ppA(random_element,x),ppB(random_element,x),
3261                                    skeleton, V_buf, sparseFail, coeffMonoms,
3262                                    Monoms);
3263          if (sparseFail)
3264            break;
3265
3266          TIMING_END_AND_PRINT (gcd_recursion,
3267                                "time for recursive call: ");
3268          DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3269        }
3270        else
3271        {
3272          CFList list;
3273          TIMING_START (gcd_recursion);
3274          bool sparseFail= false;
3275          if (LC (skeleton).inCoeffDomain())
3276            G_random_element=
3277            monicSparseInterpol (ppA (random_element, x),ppB(random_element, x),
3278                                skeleton,V_buf, sparseFail,coeffMonoms, Monoms);
3279          else
3280            G_random_element=
3281            nonMonicSparseInterpol (ppA(random_element,x),ppB(random_element,x),
3282                                    skeleton, V_buf, sparseFail, coeffMonoms,
3283                                    Monoms);
3284          if (sparseFail)
3285            break;
3286
3287          TIMING_END_AND_PRINT (gcd_recursion,
3288                                "time for recursive call: ");
3289          DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3290        }
3291
3292        if (!G_random_element.inCoeffDomain())
3293          d0= totaldegree (G_random_element, Variable(2),
3294                           Variable (G_random_element.level()));
3295        else
3296          d0= 0;
3297
3298        if (d0 == 0)
3299          return N(gcdcAcB);
3300        if (d0 >  d)
3301        {
3302          if (!find (l, random_element))
3303            l.append (random_element);
3304          continue;
3305        }
3306
3307        G_random_element=
3308        (gcdlcAlcB(random_element, x)/uni_lcoeff (G_random_element))
3309        * G_random_element;
3310
3311        if (!G_random_element.inCoeffDomain())
3312          d0= totaldegree (G_random_element, Variable(2),
3313                          Variable (G_random_element.level()));
3314        else
3315          d0= 0;
3316
3317        if (d0 <  d)
3318        {
3319          m= gcdlcAlcB;
3320          newtonPoly= 1;
3321          G_m= 0;
3322          d= d0;
3323        }
3324
3325        TIMING_START (newton_interpolation);
3326        H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
3327        TIMING_END_AND_PRINT (newton_interpolation,
3328                              "time for newton interpolation: ");
3329
3330        //termination test
3331        if (uni_lcoeff (H) == gcdlcAlcB)
3332        {
3333          cH= uni_content (H);
3334          ppH= H/cH;
3335          if (inextension)
3336          {
3337            CFList u, v;
3338            //maybe it's better to test if ppH is an element of F(\alpha) before
3339            //mapping down
3340            if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
3341            {
3342              DEBOUTLN (cerr, "ppH before mapDown= " << ppH);
3343              ppH= mapDown (ppH, prim_elem, im_prim_elem, alpha, u, v);
3344              ppH /= Lc(ppH);
3345              DEBOUTLN (cerr, "ppH after mapDown= " << ppH);
3346              return N(gcdcAcB*ppH);
3347            }
3348          }
3349          else if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
3350          {
3351            return N(gcdcAcB*ppH);
3352          }
3353        }
3354
3355        G_m= H;
3356        newtonPoly= newtonPoly*(x - random_element);
3357        m= m*(x - random_element);
3358        if (!find (l, random_element))
3359          l.append (random_element);
3360
3361      } while (1);
3362    }
3363  } while (1);
3364}
3365
3366CanonicalForm sparseGCDFp (const CanonicalForm& F, const CanonicalForm& G,
3367                           bool& topLevel, CFList& l)
3368{
3369  CanonicalForm A= F;
3370  CanonicalForm B= G;
3371  if (F.isZero() && degree(G) > 0) return G/Lc(G);
3372  else if (G.isZero() && degree (F) > 0) return F/Lc(F);
3373  else if (F.isZero() && G.isZero()) return F.genOne();
3374  if (F.inBaseDomain() || G.inBaseDomain()) return F.genOne();
3375  if (F.isUnivariate() && fdivides(F, G)) return F/Lc(F);
3376  if (G.isUnivariate() && fdivides(G, F)) return G/Lc(G);
3377  if (F == G) return F/Lc(F);
3378
3379  CFMap M,N;
3380  int best_level= myCompress (A, B, M, N, topLevel);
3381
3382  if (best_level == 0) return B.genOne();
3383
3384  A= M(A);
3385  B= M(B);
3386
3387  Variable x= Variable (1);
3388
3389  //univariate case
3390  if (A.isUnivariate() && B.isUnivariate())
3391    return N (gcd (A, B));
3392
3393  CanonicalForm cA, cB;    // content of A and B
3394  CanonicalForm ppA, ppB;    // primitive part of A and B
3395  CanonicalForm gcdcAcB;
3396
3397  cA = uni_content (A);
3398  cB = uni_content (B);
3399  gcdcAcB= gcd (cA, cB);
3400  ppA= A/cA;
3401  ppB= B/cB;
3402
3403  CanonicalForm lcA, lcB;  // leading coefficients of A and B
3404  CanonicalForm gcdlcAlcB;
3405  lcA= uni_lcoeff (ppA);
3406  lcB= uni_lcoeff (ppB);
3407
3408  if (fdivides (lcA, lcB))
3409  {
3410    if (fdivides (A, B))
3411      return F/Lc(F);
3412  }
3413  if (fdivides (lcB, lcA))
3414  {
3415    if (fdivides (B, A))
3416      return G/Lc(G);
3417  }
3418
3419  gcdlcAlcB= gcd (lcA, lcB);
3420
3421  int d= totaldegree (ppA, Variable (2), Variable (ppA.level()));
3422  int d0;
3423
3424  if (d == 0)
3425    return N(gcdcAcB);
3426
3427  d0= totaldegree (ppB, Variable (2), Variable (ppB.level()));
3428
3429  if (d0 < d)
3430    d= d0;
3431
3432  if (d == 0)
3433    return N(gcdcAcB);
3434
3435  CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH, skeleton;
3436  CanonicalForm newtonPoly= 1;
3437  m= gcdlcAlcB;
3438  G_m= 0;
3439  H= 0;
3440  bool fail= false;
3441  topLevel= false;
3442  bool inextension= false;
3443  Variable V_buf, alpha;
3444  CanonicalForm prim_elem, im_prim_elem;
3445  CFList source, dest;
3446  do //first do
3447  {
3448    if (inextension)
3449      random_element= randomElement (m, V_buf, l, fail);
3450    else
3451      random_element= FpRandomElement (m, l, fail);
3452    if (random_element == 0 && !fail)
3453    {
3454      if (!find (l, random_element))
3455        l.append (random_element);
3456      continue;
3457    }
3458
3459    if (!fail && !inextension)
3460    {
3461      CFList list;
3462      TIMING_START (gcd_recursion);
3463      G_random_element=
3464      sparseGCDFp (ppA (random_element,x), ppB (random_element,x), topLevel,
3465                   list);
3466      TIMING_END_AND_PRINT (gcd_recursion,
3467                            "time for recursive call: ");
3468      DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3469    }
3470    else if (!fail && inextension)
3471    {
3472      CFList list;
3473      TIMING_START (gcd_recursion);
3474      G_random_element=
3475      sparseGCDFq (ppA (random_element, x), ppB (random_element, x), alpha,
3476                        list, topLevel);
3477      TIMING_END_AND_PRINT (gcd_recursion,
3478                            "time for recursive call: ");
3479      DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3480    }
3481    else if (fail && !inextension)
3482    {
3483      source= CFList();
3484      dest= CFList();
3485      CFList list;
3486      CanonicalForm mipo;
3487      int deg= 2;
3488      do
3489      {
3490        mipo= randomIrredpoly (deg, x);
3491        alpha= rootOf (mipo);
3492        inextension= true;
3493        fail= false;
3494        random_element= randomElement (m, alpha, l, fail);
3495        deg++;
3496      } while (fail);
3497      V_buf= alpha;
3498      list= CFList();
3499      TIMING_START (gcd_recursion);
3500      G_random_element=
3501      sparseGCDFq (ppA (random_element, x), ppB (random_element, x), alpha,
3502                        list, topLevel);
3503      TIMING_END_AND_PRINT (gcd_recursion,
3504                            "time for recursive call: ");
3505      DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3506    }
3507    else if (fail && inextension)
3508    {
3509      source= CFList();
3510      dest= CFList();
3511
3512      Variable V_buf3= V_buf;
3513      V_buf= chooseExtension (V_buf);
3514      bool prim_fail= false;
3515      Variable V_buf2;
3516      CanonicalForm prim_elem, im_prim_elem;
3517      prim_elem= primitiveElement (alpha, V_buf2, prim_fail);
3518
3519      if (V_buf3 != alpha)
3520      {
3521        m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
3522        G_m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
3523        newtonPoly= mapDown (newtonPoly, prim_elem, im_prim_elem, alpha, source,
3524                             dest);
3525        ppA= mapDown (ppA, prim_elem, im_prim_elem, alpha, source, dest);
3526        ppB= mapDown (ppB, prim_elem, im_prim_elem, alpha, source, dest);
3527        gcdlcAlcB= mapDown (gcdlcAlcB, prim_elem, im_prim_elem, alpha, source,
3528                            dest);
3529        for (CFListIterator i= l; i.hasItem(); i++)
3530          i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, alpha,
3531                                source, dest);
3532      }
3533
3534      ASSERT (!prim_fail, "failure in integer factorizer");
3535      if (prim_fail)
3536        ; //ERROR
3537      else
3538        im_prim_elem= mapPrimElem (prim_elem, alpha, V_buf);
3539
3540      DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (alpha));
3541      DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf2));
3542
3543      for (CFListIterator i= l; i.hasItem(); i++)
3544        i.getItem()= mapUp (i.getItem(), alpha, V_buf, prim_elem,
3545                             im_prim_elem, source, dest);
3546      m= mapUp (m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3547      G_m= mapUp (G_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3548      newtonPoly= mapUp (newtonPoly, alpha, V_buf, prim_elem, im_prim_elem,
3549                          source, dest);
3550      ppA= mapUp (ppA, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3551      ppB= mapUp (ppB, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3552      gcdlcAlcB= mapUp (gcdlcAlcB, alpha, V_buf, prim_elem, im_prim_elem,
3553                         source, dest);
3554      fail= false;
3555      random_element= randomElement (m, V_buf, l, fail );
3556      DEBOUTLN (cerr, "fail= " << fail);
3557      CFList list;
3558      TIMING_START (gcd_recursion);
3559      G_random_element=
3560      sparseGCDFq (ppA (random_element, x), ppB (random_element, x), V_buf,
3561                        list, topLevel);
3562      TIMING_END_AND_PRINT (gcd_recursion,
3563                            "time for recursive call: ");
3564      DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3565    }
3566
3567    if (!G_random_element.inCoeffDomain())
3568      d0= totaldegree (G_random_element, Variable(2),
3569                       Variable (G_random_element.level()));
3570    else
3571      d0= 0;
3572
3573    if (d0 == 0)
3574      return N(gcdcAcB);
3575    if (d0 >  d)
3576    {
3577      if (!find (l, random_element))
3578        l.append (random_element);
3579      continue;
3580    }
3581
3582    G_random_element=
3583    (gcdlcAlcB(random_element, x)/uni_lcoeff (G_random_element))
3584    * G_random_element;
3585
3586    skeleton= G_random_element;
3587
3588    if (!G_random_element.inCoeffDomain())
3589      d0= totaldegree (G_random_element, Variable(2),
3590                       Variable (G_random_element.level()));
3591    else
3592      d0= 0;
3593
3594    if (d0 <  d)
3595    {
3596      m= gcdlcAlcB;
3597      newtonPoly= 1;
3598      G_m= 0;
3599      d= d0;
3600    }
3601
3602    H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
3603
3604    if (uni_lcoeff (H) == gcdlcAlcB)
3605    {
3606      cH= uni_content (H);
3607      ppH= H/cH;
3608      ppH /= Lc (ppH);
3609      DEBOUTLN (cerr, "ppH= " << ppH);
3610
3611      if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
3612        return N(gcdcAcB*ppH);
3613    }
3614    G_m= H;
3615    newtonPoly= newtonPoly*(x - random_element);
3616    m= m*(x - random_element);
3617    if (!find (l, random_element))
3618      l.append (random_element);
3619
3620    if ((getCharacteristic() > 3 && size (skeleton) < 200))
3621    {
3622      CFArray Monoms;
3623      CFArray* coeffMonoms;
3624
3625      do //second do
3626      {
3627        if (inextension)
3628          random_element= randomElement (m, alpha, l, fail);
3629        else
3630          random_element= FpRandomElement (m, l, fail);
3631        if (random_element == 0 && !fail)
3632        {
3633          if (!find (l, random_element))
3634            l.append (random_element);
3635          continue;
3636        }
3637
3638        bool sparseFail= false;
3639        if (!fail && !inextension)
3640        {
3641          CFList list;
3642          TIMING_START (gcd_recursion);
3643          if (LC (skeleton).inCoeffDomain())
3644            G_random_element=
3645            monicSparseInterpol(ppA(random_element, x), ppB (random_element, x),
3646                                skeleton, Variable(1), sparseFail, coeffMonoms,
3647                                Monoms);
3648          else
3649            G_random_element=
3650            nonMonicSparseInterpol(ppA(random_element,x), ppB(random_element,x),
3651                                    skeleton, Variable (1), sparseFail,
3652                                    coeffMonoms, Monoms);
3653          TIMING_END_AND_PRINT (gcd_recursion,
3654                                "time for recursive call: ");
3655          DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3656        }
3657        else if (!fail && inextension)
3658        {
3659          CFList list;
3660          TIMING_START (gcd_recursion);
3661          if (LC (skeleton).inCoeffDomain())
3662            G_random_element=
3663            monicSparseInterpol(ppA (random_element,x), ppB (random_element, x),
3664                                skeleton, alpha, sparseFail, coeffMonoms,
3665                                Monoms);
3666          else
3667            G_random_element=
3668            nonMonicSparseInterpol(ppA(random_element,x), ppB(random_element,x),
3669                                   skeleton, alpha, sparseFail, coeffMonoms,
3670                                   Monoms);
3671          TIMING_END_AND_PRINT (gcd_recursion,
3672                                "time for recursive call: ");
3673          DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3674        }
3675        else if (fail && !inextension)
3676        {
3677          source= CFList();
3678          dest= CFList();
3679          CFList list;
3680          CanonicalForm mipo;
3681          int deg= 2;
3682          do
3683          {
3684            mipo= randomIrredpoly (deg, x);
3685            alpha= rootOf (mipo);
3686            inextension= true;
3687            fail= false;
3688            random_element= randomElement (m, alpha, l, fail);
3689            deg++;
3690          } while (fail);
3691          V_buf= alpha;
3692          list= CFList();
3693          TIMING_START (gcd_recursion);
3694          if (LC (skeleton).inCoeffDomain())
3695            G_random_element=
3696            monicSparseInterpol (ppA (random_element,x), ppB (random_element,x),
3697                                skeleton, alpha, sparseFail, coeffMonoms,
3698                                Monoms);
3699          else
3700            G_random_element=
3701            nonMonicSparseInterpol(ppA(random_element,x), ppB(random_element,x),
3702                                   skeleton, alpha, sparseFail, coeffMonoms,
3703                                   Monoms);
3704          TIMING_END_AND_PRINT (gcd_recursion,
3705                                "time for recursive call: ");
3706          DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3707        }
3708        else if (fail && inextension)
3709        {
3710          source= CFList();
3711          dest= CFList();
3712
3713          Variable V_buf3= V_buf;
3714          V_buf= chooseExtension (V_buf);
3715          bool prim_fail= false;
3716          Variable V_buf2;
3717          CanonicalForm prim_elem, im_prim_elem;
3718          prim_elem= primitiveElement (alpha, V_buf2, prim_fail);
3719
3720          if (V_buf3 != alpha)
3721          {
3722            m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
3723            G_m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
3724            newtonPoly= mapDown (newtonPoly, prim_elem, im_prim_elem, alpha,
3725                                 source, dest);
3726            ppA= mapDown (ppA, prim_elem, im_prim_elem, alpha, source, dest);
3727            ppB= mapDown (ppB, prim_elem, im_prim_elem, alpha, source, dest);
3728            gcdlcAlcB= mapDown (gcdlcAlcB, prim_elem, im_prim_elem, alpha,
3729                                source, dest);
3730            for (CFListIterator i= l; i.hasItem(); i++)
3731              i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, alpha,
3732                                    source, dest);
3733          }
3734
3735          ASSERT (!prim_fail, "failure in integer factorizer");
3736          if (prim_fail)
3737            ; //ERROR
3738          else
3739            im_prim_elem= mapPrimElem (prim_elem, alpha, V_buf);
3740
3741          DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (alpha));
3742          DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf2));
3743
3744          for (CFListIterator i= l; i.hasItem(); i++)
3745            i.getItem()= mapUp (i.getItem(), alpha, V_buf, prim_elem,
3746                                im_prim_elem, source, dest);
3747          m= mapUp (m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3748          G_m= mapUp (G_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3749          newtonPoly= mapUp (newtonPoly, alpha, V_buf, prim_elem, im_prim_elem,
3750                              source, dest);
3751          ppA= mapUp (ppA, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3752          ppB= mapUp (ppB, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3753          gcdlcAlcB= mapUp (gcdlcAlcB, alpha, V_buf, prim_elem, im_prim_elem,
3754                            source, dest);
3755          fail= false;
3756          random_element= randomElement (m, V_buf, l, fail );
3757          DEBOUTLN (cerr, "fail= " << fail);
3758          CFList list;
3759          TIMING_START (gcd_recursion);
3760          if (LC (skeleton).inCoeffDomain())
3761            G_random_element=
3762            monicSparseInterpol (ppA (random_element, x), ppB (random_element, x),
3763                                skeleton, V_buf, sparseFail, coeffMonoms,
3764                                Monoms);
3765          else
3766            G_random_element=
3767            nonMonicSparseInterpol (ppA(random_element,x), ppB(random_element,x),
3768                                    skeleton, V_buf, sparseFail,
3769                                    coeffMonoms, Monoms);
3770          TIMING_END_AND_PRINT (gcd_recursion,
3771                                "time for recursive call: ");
3772          DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3773        }
3774
3775        if (sparseFail)
3776          break;
3777
3778        if (!G_random_element.inCoeffDomain())
3779          d0= totaldegree (G_random_element, Variable(2),
3780                           Variable (G_random_element.level()));
3781        else
3782          d0= 0;
3783
3784        if (d0 == 0)
3785          return N(gcdcAcB);
3786        if (d0 >  d)
3787        {
3788          if (!find (l, random_element))
3789            l.append (random_element);
3790          continue;
3791        }
3792
3793        G_random_element=
3794        (gcdlcAlcB(random_element, x)/uni_lcoeff (G_random_element))
3795        * G_random_element;
3796
3797        if (!G_random_element.inCoeffDomain())
3798          d0= totaldegree (G_random_element, Variable(2),
3799                           Variable (G_random_element.level()));
3800        else
3801          d0= 0;
3802
3803        if (d0 <  d)
3804        {
3805          m= gcdlcAlcB;
3806          newtonPoly= 1;
3807          G_m= 0;
3808          d= d0;
3809        }
3810
3811        TIMING_START (newton_interpolation);
3812        H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
3813        TIMING_END_AND_PRINT (newton_interpolation,
3814                              "time for newton interpolation: ");
3815
3816        //termination test
3817        if (uni_lcoeff (H) == gcdlcAlcB)
3818        {
3819          cH= uni_content (H);
3820          ppH= H/cH;
3821          ppH /= Lc (ppH);
3822          DEBOUTLN (cerr, "ppH= " << ppH);
3823          if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
3824            return N(gcdcAcB*ppH);
3825        }
3826
3827        G_m= H;
3828        newtonPoly= newtonPoly*(x - random_element);
3829        m= m*(x - random_element);
3830        if (!find (l, random_element))
3831          l.append (random_element);
3832
3833      } while (1); //end of second do
3834    }
3835    else
3836      return N(gcdcAcB*GCD_small_p (ppA, ppB));
3837  } while (1); //end of first do
3838}
3839
3840static inline
3841int compress4EZGCD (const CanonicalForm& F, const CanonicalForm& G, CFMap & M,
3842                    CFMap & N, int& both_non_zero)
3843{
3844  int n= tmax (F.level(), G.level());
3845  int * degsf= new int [n + 1];
3846  int * degsg= new int [n + 1];
3847
3848  for (int i = 0; i <= n; i++)
3849    degsf[i]= degsg[i]= 0;
3850
3851  degsf= degrees (F, degsf);
3852  degsg= degrees (G, degsg);
3853
3854  both_non_zero= 0;
3855  int f_zero= 0;
3856  int g_zero= 0;
3857
3858  for (int i= 1; i <= n; i++)
3859  {
3860    if (degsf[i] != 0 && degsg[i] != 0)
3861    {
3862      both_non_zero++;
3863      continue;
3864    }
3865    if (degsf[i] == 0 && degsg[i] != 0 && i <= G.level())
3866    {
3867      f_zero++;
3868      continue;
3869    }
3870    if (degsg[i] == 0 && degsf[i] && i <= F.level())
3871    {
3872      g_zero++;
3873      continue;
3874    }
3875  }
3876
3877  if (both_non_zero == 0)
3878  {
3879    delete [] degsf;
3880    delete [] degsg;
3881    return 0;
3882  }
3883
3884  // map Variables which do not occur in both polynomials to higher levels
3885  int k= 1;
3886  int l= 1;
3887  for (int i= 1; i <= n; i++)
3888  {
3889    if (degsf[i] != 0 && degsg[i] == 0 && i <= F.level())
3890    {
3891      if (k + both_non_zero != i)
3892      {
3893        M.newpair (Variable (i), Variable (k + both_non_zero));
3894        N.newpair (Variable (k + both_non_zero), Variable (i));
3895      }
3896      k++;
3897    }
3898    if (degsf[i] == 0 && degsg[i] != 0 && i <= G.level())
3899    {
3900      if (l + g_zero + both_non_zero != i)
3901      {
3902        M.newpair (Variable (i), Variable (l + g_zero + both_non_zero));
3903        N.newpair (Variable (l + g_zero + both_non_zero), Variable (i));
3904      }
3905      l++;
3906    }
3907  }
3908
3909  // sort Variables x_{i} in decreasing order of
3910  // min(deg_{x_{i}}(f),deg_{x_{i}}(g))
3911  int m= tmin (F.level(), G.level());
3912  int max_min_deg;
3913  k= both_non_zero;
3914  l= 0;
3915  int i= 1;
3916  while (k > 0)
3917  {
3918    max_min_deg= tmin (degsf[i], degsg[i]);
3919    while (max_min_deg == 0)
3920    {
3921      i++;
3922      max_min_deg= tmin (degsf[i], degsg[i]);
3923    }
3924    for (int j= i + 1; j <= m; j++)
3925    {
3926      if ((tmin (degsf[j],degsg[j]) < max_min_deg) &&
3927          (tmin (degsf[j], degsg[j]) != 0))
3928      {
3929        max_min_deg= tmin (degsf[j], degsg[j]);
3930        l= j;
3931      }
3932    }
3933
3934    if (l != 0)
3935    {
3936      if (l != k)
3937      {
3938        M.newpair (Variable (l), Variable(k));
3939        N.newpair (Variable (k), Variable(l));
3940        degsf[l]= 0;
3941        degsg[l]= 0;
3942        l= 0;
3943      }
3944      else
3945      {
3946        degsf[l]= 0;
3947        degsg[l]= 0;
3948        l= 0;
3949      }
3950    }
3951    else if (l == 0)
3952    {
3953      if (i != k)
3954      {
3955        M.newpair (Variable (i), Variable (k));
3956        N.newpair (Variable (k), Variable (i));
3957        degsf[i]= 0;
3958        degsg[i]= 0;
3959      }
3960      else
3961      {
3962        degsf[i]= 0;
3963        degsg[i]= 0;
3964      }
3965      i++;
3966    }
3967    k--;
3968  }
3969
3970  delete [] degsf;
3971  delete [] degsg;
3972
3973  return both_non_zero;
3974}
3975
3976static inline
3977CanonicalForm myShift2Zero (const CanonicalForm& F, CFList& Feval,
3978                            const CFList& evaluation)
3979{
3980  CanonicalForm A= F;
3981  int k= 2;
3982  for (CFListIterator i= evaluation; i.hasItem(); i++, k++)
3983    A= A (Variable (k) + i.getItem(), k);
3984
3985  CanonicalForm buf= A;
3986  Feval= CFList();
3987  Feval.append (buf);
3988  for (k= evaluation.length() + 1; k > 2; k--)
3989  {
3990    buf= mod (buf, Variable (k));
3991    Feval.insert (buf);
3992  }
3993  return A;
3994}
3995
3996static inline
3997CanonicalForm myReverseShift (const CanonicalForm& F, const CFList& evaluation)
3998{
3999  int l= evaluation.length() + 1;
4000  CanonicalForm result= F;
4001  CFListIterator j= evaluation;
4002  for (int i= 2; i < l + 1; i++, j++)
4003  {
4004    if (F.level() < i)
4005      continue;
4006    result= result (Variable (i) - j.getItem(), i);
4007  }
4008  return result;
4009}
4010
4011static inline
4012Evaluation optimize4Lift (const CanonicalForm& F, CFMap & M,
4013                    CFMap & N, const Evaluation& A)
4014{
4015  int n= F.level();
4016  int * degsf= new int [n + 1];
4017
4018  for (int i = 0; i <= n; i++)
4019    degsf[i]= 0;
4020
4021  degsf= degrees (F, degsf);
4022
4023  Evaluation result= Evaluation (A.min(), A.max());
4024  ASSERT (A.min() == 2, "expected A.min() == 2");
4025  int max_deg;
4026  int k= n;
4027  int l= 1;
4028  int i= 2;
4029  int pos= 2;
4030  while (k > 1)
4031  {
4032    max_deg= degsf [i];
4033    while (max_deg == 0)
4034    {
4035      i++;
4036      max_deg= degsf [i];
4037    }
4038    l= i;
4039    for (int j= i + 1; j <=  n; j++)
4040    {
4041      if (degsf[j] > max_deg)
4042      {
4043        max_deg= degsf[j];
4044        l= j;
4045      }
4046    }
4047
4048    if (l <= n)
4049    {
4050      if (l != pos)
4051      {
4052        result.setValue (pos, A [l]);
4053        M.newpair (Variable (l), Variable (pos));
4054        N.newpair (Variable (pos), Variable (l));
4055        degsf[l]= 0;
4056        l= 2;
4057        if (k == 2 && n == 3)
4058        {
4059          result.setValue (l, A [pos]);
4060          M.newpair (Variable (pos), Variable (l));
4061          N.newpair (Variable (l), Variable (pos));
4062          degsf[pos]= 0;
4063        }
4064      }
4065      else
4066      {
4067        result.setValue (l, A [l]);
4068        degsf [l]= 0;
4069      }
4070    }
4071    pos++;
4072    k--;
4073    l= 2;
4074  }
4075
4076  delete [] degsf;
4077
4078  return result;
4079}
4080
4081static inline
4082int Hensel_P (const CanonicalForm & UU, CFArray & G, const Evaluation & AA,
4083              const CFArray& LeadCoeffs )
4084{
4085  CFList factors;
4086  factors.append (G[1]);
4087  factors.append (G[2]);
4088
4089  CFMap NN, MM;
4090  Evaluation A= optimize4Lift (UU, MM, NN, AA);
4091
4092  CanonicalForm U= MM (UU);
4093  CFArray LCs= CFArray (1,2);
4094  LCs [1]= MM (LeadCoeffs [1]);
4095  LCs [2]= MM (LeadCoeffs [2]);
4096
4097  CFList evaluation;
4098  long termEstimate= size (U);
4099  for (int i= A.min(); i <= A.max(); i++)
4100  {
4101    if (!A[i].isZero() && (getCharacteristic() > degree (U,i))) //TODO find a good estimate for getCharacteristic() <= degree (U,i)
4102    {
4103      termEstimate *= degree (U,i)*2;
4104      termEstimate /= 3;
4105    }
4106    evaluation.append (A [i]);
4107  }
4108  if (termEstimate/getNumVars(U) > 500)
4109    return -1;
4110  CFList UEval;
4111  CanonicalForm shiftedU= myShift2Zero (U, UEval, evaluation);
4112
4113  if (size (shiftedU)/getNumVars (U) > 500)
4114    return -1;
4115
4116  CFArray shiftedLCs= CFArray (2);
4117  CFList shiftedLCsEval1, shiftedLCsEval2;
4118  shiftedLCs[0]= myShift2Zero (LCs[1], shiftedLCsEval1, evaluation);
4119  shiftedLCs[1]= myShift2Zero (LCs[2], shiftedLCsEval2, evaluation);
4120  factors.insert (1);
4121  int liftBound= degree (UEval.getLast(), 2) + 1;
4122  CFArray Pi;
4123  CFMatrix M= CFMatrix (liftBound, factors.length() - 1);
4124  CFList diophant;
4125  CFArray lcs= CFArray (2);
4126  lcs [0]= shiftedLCsEval1.getFirst();
4127  lcs [1]= shiftedLCsEval2.getFirst();
4128  nonMonicHenselLift12 (UEval.getFirst(), factors, liftBound, Pi, diophant, M,
4129                        lcs, false);
4130
4131  for (CFListIterator i= factors; i.hasItem(); i++)
4132  {
4133    if (!fdivides (i.getItem(), UEval.getFirst()))
4134      return 0;
4135  }
4136
4137  int * liftBounds;
4138  bool noOneToOne= false;
4139  if (U.level() > 2)
4140  {
4141    liftBounds= new int [U.level() - 1]; /* index: 0.. U.level()-2 */
4142    liftBounds[0]= liftBound;
4143    for (int i= 1; i < U.level() - 1; i++)
4144      liftBounds[i]= degree (shiftedU, Variable (i + 2)) + 1;
4145    factors= nonMonicHenselLift2 (UEval, factors, liftBounds, U.level() - 1,
4146                                  false, shiftedLCsEval1, shiftedLCsEval2, Pi,
4147                                  diophant, noOneToOne);
4148    delete [] liftBounds;
4149    if (noOneToOne)
4150      return 0;
4151  }
4152  G[1]= factors.getFirst();
4153  G[2]= factors.getLast();
4154  G[1]= myReverseShift (G[1], evaluation);
4155  G[2]= myReverseShift (G[2], evaluation);
4156  G[1]= NN (G[1]);
4157  G[2]= NN (G[2]);
4158  return 1;
4159}
4160
4161static inline
4162bool findeval_P (const CanonicalForm & F, const CanonicalForm & G,
4163                 CanonicalForm & Fb, CanonicalForm & Gb, CanonicalForm & Db,
4164                 REvaluation & b, int delta, int degF, int degG, int maxeval,
4165                 int & count, int& k, int bound, int& l)
4166{
4167  if( count == 0 && delta != 0)
4168  {
4169    if( count++ > maxeval )
4170      return false;
4171  }
4172  if (count > 0)
4173  {
4174    b.nextpoint(k);
4175    if (k == 0)
4176      k++;
4177    l++;
4178    if (l > bound)
4179    {
4180      l= 1;
4181      k++;
4182      if (k > tmax (F.level(), G.level()) - 1)
4183        return false;
4184      b.nextpoint (k);
4185    }
4186    if (count++ > maxeval)
4187      return false;
4188  }
4189  while( true )
4190  {
4191    Fb = b( F );
4192    if( degree( Fb, 1 ) == degF )
4193    {
4194      Gb = b( G );
4195      if( degree( Gb, 1 ) == degG )
4196      {
4197        Db = gcd( Fb, Gb );
4198        if( delta > 0 )
4199        {
4200          if( degree( Db, 1 ) <= delta )
4201            return true;
4202        }
4203        else
4204          return true;
4205      }
4206    }
4207    if (k == 0)
4208      k++;
4209    b.nextpoint(k);
4210    l++;
4211    if (l > bound)
4212    {
4213      l= 1;
4214      k++;
4215      if (k > tmax (F.level(), G.level()) - 1)
4216        return false;
4217      b.nextpoint (k);
4218    }
4219    if( count++ > maxeval )
4220      return false;
4221  }
4222}
4223
4224// parameters for heuristic
4225static int maxNumEval= 200;
4226static int sizePerVars1= 500; //try dense gcd if size/#variables is bigger
4227
4228/// Extended Zassenhaus GCD for finite fields
4229CanonicalForm EZGCD_P( const CanonicalForm & FF, const CanonicalForm & GG )
4230{
4231  if (FF.isZero() && degree(GG) > 0) return GG/Lc(GG);
4232  else if (GG.isZero() && degree (FF) > 0) return FF/Lc(FF);
4233  else if (FF.isZero() && GG.isZero()) return FF.genOne();
4234  if (FF.inBaseDomain() || GG.inBaseDomain()) return FF.genOne();
4235  if (FF.isUnivariate() && fdivides(FF, GG)) return FF/Lc(FF);
4236  if (GG.isUnivariate() && fdivides(GG, FF)) return GG/Lc(GG);
4237  if (FF == GG) return FF/Lc(FF);
4238
4239  int maxNumVars= tmax (getNumVars (FF), getNumVars (GG));
4240  Variable a, oldA;
4241  int sizeF= size (FF);
4242  int sizeG= size (GG);
4243
4244  if (sizeF/maxNumVars > sizePerVars1 && sizeG/maxNumVars > sizePerVars1)
4245  {
4246    if (hasFirstAlgVar (FF, a) || hasFirstAlgVar (GG, a))
4247      return GCD_Fp_extension (FF, GG, a);
4248    else if (CFFactory::gettype() == GaloisFieldDomain)
4249      return GCD_GF (FF, GG);
4250    else
4251      return GCD_small_p (FF, GG);
4252  }
4253
4254  CanonicalForm F, G, f, g, d, Fb, Gb, Db, Fbt, Gbt, Dbt, B0, B, D0, lcF, lcG,
4255                lcD;
4256  CFArray DD( 1, 2 ), lcDD( 1, 2 );
4257  int degF, degG, delta, count;
4258  int maxeval;
4259  maxeval= tmin((getCharacteristic()/
4260                (int)(ilog2(getCharacteristic())*log2exp))*2, maxNumEval);
4261  count= 0; // number of eval. used
4262  REvaluation b, bt;
4263  int gcdfound = 0;
4264  Variable x = Variable(1);
4265
4266  F= FF;
4267  G= GG;
4268
4269  CFMap M,N;
4270  int smallestDegLev;
4271  TIMING_START (ez_p_compress)
4272  int best_level= compress4EZGCD (F, G, M, N, smallestDegLev);
4273
4274  if (best_level == 0) return G.genOne();
4275
4276  F= M (F);
4277  G= M (G);
4278  TIMING_END_AND_PRINT (ez_p_compress, "time for compression in EZ_P: ")
4279
4280  TIMING_START (ez_p_content)
4281  f = content( F, x ); g = content( G, x );
4282  d = gcd( f, g );
4283  F /= f; G /= g;
4284  TIMING_END_AND_PRINT (ez_p_content, "time to extract content in EZ_P: ")
4285
4286  if( F.isUnivariate() && G.isUnivariate() )
4287  {
4288    if( F.mvar() == G.mvar() )
4289      d *= gcd( F, G );
4290    else
4291      return N (d);
4292    return N (d);
4293  }
4294  if ( F.isUnivariate())
4295  {
4296    g= content (G,G.mvar());
4297    return N(d*gcd(F,g));
4298  }
4299  if ( G.isUnivariate())
4300  {
4301    f= content (F,F.mvar());
4302    return N(d*gcd(G,f));
4303  }
4304
4305  maxNumVars= tmax (getNumVars (F), getNumVars (G));
4306  sizeF= size (F);
4307  sizeG= size (G);
4308
4309  if (sizeF/maxNumVars > sizePerVars1 && sizeG/maxNumVars > sizePerVars1)
4310  {
4311    if (hasFirstAlgVar (F, a) || hasFirstAlgVar (G, a))
4312      return N (d*GCD_Fp_extension (F, G, a));
4313    else if (CFFactory::gettype() == GaloisFieldDomain)
4314      return N (d*GCD_GF (F, G));
4315    else
4316      return N (d*GCD_small_p (F, G));
4317  }
4318
4319  int dummy= 0;
4320  if( gcd_test_one( F, G, false, dummy ) )
4321  {
4322    return N (d);
4323  }
4324
4325  bool passToGF= false;
4326  bool extOfExt= false;
4327  int p= getCharacteristic();
4328  bool algExtension= (hasFirstAlgVar(F,a) || hasFirstAlgVar(G,a));
4329  int k= 1;
4330  CanonicalForm primElem, imPrimElem;
4331  CFList source, dest;
4332  if (p < 50 && CFFactory::gettype() != GaloisFieldDomain && !algExtension)
4333  {
4334    if (p == 2)
4335      setCharacteristic (2, 12, 'Z');
4336    else if (p == 3)
4337      setCharacteristic (3, 4, 'Z');
4338    else if (p == 5 || p == 7)
4339      setCharacteristic (p, 3, 'Z');
4340    else
4341      setCharacteristic (p, 2, 'Z');
4342    passToGF= true;
4343    F= F.mapinto();
4344    G= G.mapinto();
4345    maxeval= 2*ipower (p, getGFDegree());
4346  }
4347  else if (CFFactory::gettype() == GaloisFieldDomain &&
4348           ipower (p , getGFDegree()) < 50)
4349  {
4350    k= getGFDegree();
4351    if (ipower (p, 2*k) > 50)
4352      setCharacteristic (p, 2*k, gf_name);
4353    else
4354      setCharacteristic (p, 3*k, gf_name);
4355    F= GFMapUp (F, k);
4356    G= GFMapUp (G, k);
4357    maxeval= tmin (2*ipower (p, getGFDegree()), maxNumEval);
4358  }
4359  else if (p < 50 && algExtension && CFFactory::gettype() != GaloisFieldDomain)
4360  {
4361    int d= degree (getMipo (a));
4362    oldA= a;
4363    Variable v2;
4364    if (p == 2 && d < 6)
4365    {
4366      if (fac_NTL_char != p)
4367      {
4368        fac_NTL_char= p;
4369        zz_p::init (p);
4370      }
4371      bool primFail= false;
4372      Variable vBuf;
4373      primElem= primitiveElement (a, vBuf, primFail);
4374      ASSERT (!primFail, "failure in integer factorizer");
4375      if (d < 3)
4376      {
4377        zz_pX NTLIrredpoly;
4378        BuildIrred (NTLIrredpoly, d*3);
4379        CanonicalForm newMipo= convertNTLzzpX2CF (NTLIrredpoly, Variable (1));
4380        v2= rootOf (newMipo);
4381      }
4382      else
4383      {
4384        zz_pX NTLIrredpoly;
4385        BuildIrred (NTLIrredpoly, d*2);
4386        CanonicalForm newMipo= convertNTLzzpX2CF (NTLIrredpoly, Variable (1));
4387        v2= rootOf (newMipo);
4388      }
4389      imPrimElem= mapPrimElem (primElem, a, v2);
4390      extOfExt= true;
4391    }
4392    else if ((p == 3 && d < 4) || ((p == 5 || p == 7) && d < 3))
4393    {
4394      if (fac_NTL_char != p)
4395      {
4396        fac_NTL_char= p;
4397        zz_p::init (p);
4398      }
4399      bool primFail= false;
4400      Variable vBuf;
4401      primElem= primitiveElement (a, vBuf, primFail);
4402      ASSERT (!primFail, "failure in integer factorizer");
4403      zz_pX NTLIrredpoly;
4404      BuildIrred (NTLIrredpoly, d*2);
4405      CanonicalForm newMipo= convertNTLzzpX2CF (NTLIrredpoly, Variable (1));
4406      v2= rootOf (newMipo);
4407      imPrimElem= mapPrimElem (primElem, a, v2);
4408      extOfExt= true;
4409    }
4410    if (extOfExt)
4411    {
4412      maxeval= tmin (2*ipower (p, degree (getMipo (v2))), maxNumEval);
4413      F= mapUp (F, a, v2, primElem, imPrimElem, source, dest);
4414      G= mapUp (G, a, v2, primElem, imPrimElem, source, dest);
4415      a= v2;
4416    }
4417  }
4418
4419  lcF = LC( F, x ); lcG = LC( G, x );
4420  lcD = gcd( lcF, lcG );
4421
4422  delta = 0;
4423  degF = degree( F, x ); degG = degree( G, x );
4424
4425  if(hasFirstAlgVar(G,a))
4426    b = REvaluation( 2, tmax(F.level(), G.level()), AlgExtRandomF( a ) );
4427  else
4428  { // both not in extension given by algebraic variable
4429    if (CFFactory::gettype() != GaloisFieldDomain)
4430      b = REvaluation( 2, tmax(F.level(), G.level()), FFRandom() );
4431    else
4432      b = REvaluation( 2, tmax(F.level(), G.level()), GFRandom() );
4433  }
4434
4435  CanonicalForm cand, contcand;
4436  CanonicalForm result;
4437  int o, t;
4438  o= 0;
4439  t= 1;
4440  int goodPointCount= 0;
4441  while( !gcdfound )
4442  {
4443    TIMING_START (ez_p_eval);
4444    if( !findeval_P( F, G, Fb, Gb, Db, b, delta, degF, degG, maxeval, count, o,
4445         maxeval/maxNumVars, t ))
4446    { // too many eval. used --> try another method
4447      Off (SW_USE_EZGCD_P);
4448      result= gcd (F,G);
4449      On (SW_USE_EZGCD_P);
4450      if (passToGF)
4451      {
4452        CanonicalForm mipo= gf_mipo;
4453        setCharacteristic (p);
4454        Variable alpha= rootOf (mipo.mapinto());
4455        result= GF2FalphaRep (result, alpha);
4456      }
4457      if (k > 1)
4458      {
4459        result= GFMapDown (result, k);
4460        setCharacteristic (p, k, gf_name);
4461      }
4462      if (extOfExt)
4463        result= mapDown (result, primElem, imPrimElem, oldA, dest, source);
4464      return N (d*result);
4465    }
4466    TIMING_END_AND_PRINT (ez_p_eval, "time for eval point search in EZ_P1: ");
4467    delta = degree( Db );
4468    if( delta == 0 )
4469    {
4470      if (passToGF)
4471        setCharacteristic (p);
4472      if (k > 1)
4473        setCharacteristic (p, k, gf_name);
4474      return N (d);
4475    }
4476    while( true )
4477    {
4478      bt = b;
4479      TIMING_START (ez_p_eval);
4480      if( !findeval_P(F,G,Fbt,Gbt,Dbt, bt, delta, degF, degG, maxeval, count, o,
4481           maxeval/maxNumVars, t ))
4482      { // too many eval. used --> try another method
4483        Off (SW_USE_EZGCD_P);
4484        result= gcd (F,G);
4485        On (SW_USE_EZGCD_P);
4486        if (passToGF)
4487        {
4488          CanonicalForm mipo= gf_mipo;
4489          setCharacteristic (p);
4490          Variable alpha= rootOf (mipo.mapinto());
4491          result= GF2FalphaRep (result, alpha);
4492        }
4493        if (k > 1)
4494        {
4495          result= GFMapDown (result, k);
4496          setCharacteristic (p, k, gf_name);
4497        }
4498        if (extOfExt)
4499          result= mapDown (result, primElem, imPrimElem, oldA, dest, source);
4500        return N (d*result);
4501      }
4502      TIMING_END_AND_PRINT (ez_p_eval, "time for eval point search in EZ_P2: ");
4503      int dd = degree( Dbt );
4504      if( dd == 0 )
4505      {
4506        if (passToGF)
4507          setCharacteristic (p);
4508        if (k > 1)
4509          setCharacteristic (p, k, gf_name);
4510        return N (d);
4511      }
4512      if( dd == delta )
4513      {
4514        goodPointCount++;
4515        if (goodPointCount == 5)
4516          break;
4517      }
4518      if( dd < delta )
4519      {
4520        goodPointCount= 0;
4521        delta = dd;
4522        b = bt;
4523        Db = Dbt; Fb = Fbt; Gb = Gbt;
4524      }
4525      if (delta == degF)
4526      {
4527        if (degF <= degG && fdivides (F, G))
4528        {
4529          if (passToGF)
4530          {
4531            CanonicalForm mipo= gf_mipo;
4532            setCharacteristic (p);
4533            Variable alpha= rootOf (mipo.mapinto());
4534            F= GF2FalphaRep (F, alpha);
4535          }
4536          if (k > 1)
4537          {
4538            F= GFMapDown (F, k);
4539            setCharacteristic (p, k, gf_name);
4540          }
4541          if (extOfExt)
4542            F= mapDown (F, primElem, imPrimElem, oldA, dest, source);
4543          return N (d*F);
4544        }
4545        else
4546          delta--;
4547      }
4548      else if (delta == degG)
4549      {
4550        if (degG <= degF && fdivides (G, F))
4551        {
4552          if (passToGF)
4553          {
4554            CanonicalForm mipo= gf_mipo;
4555            setCharacteristic (p);
4556            Variable alpha= rootOf (mipo.mapinto());
4557            G= GF2FalphaRep (G, alpha);
4558          }
4559          if (k > 1)
4560          {
4561            G= GFMapDown (G, k);
4562            setCharacteristic (p, k, gf_name);
4563          }
4564          if (extOfExt)
4565            G= mapDown (G, primElem, imPrimElem, oldA, dest, source);
4566          return N (d*G);
4567        }
4568        else
4569          delta--;
4570      }
4571    }
4572    if( delta != degF && delta != degG )
4573    {
4574      bool B_is_F;
4575      CanonicalForm xxx1, xxx2;
4576      CanonicalForm buf;
4577      DD[1] = Fb / Db;
4578      buf= Gb/Db;
4579      xxx1 = gcd( DD[1], Db );
4580      xxx2 = gcd( buf, Db );
4581      if (((xxx1.inCoeffDomain() && xxx2.inCoeffDomain()) &&
4582          (size (F) <= size (G)))
4583          || (xxx1.inCoeffDomain() && !xxx2.inCoeffDomain()))
4584      {
4585        B = F;
4586        DD[2] = Db;
4587        lcDD[1] = lcF;
4588        lcDD[2] = lcD;
4589        B_is_F = true;
4590      }
4591      else if (((xxx1.inCoeffDomain() && xxx2.inCoeffDomain()) &&
4592               (size (G) < size (F)))
4593               || (!xxx1.inCoeffDomain() && xxx2.inCoeffDomain()))
4594      {
4595        DD[1] = buf;
4596        B = G;
4597        DD[2] = Db;
4598        lcDD[1] = lcG;
4599        lcDD[2] = lcD;
4600        B_is_F = false;
4601      }
4602      else // special case handling
4603      {
4604        Off (SW_USE_EZGCD_P);
4605        result= gcd (F,G);
4606        On (SW_USE_EZGCD_P);
4607        if (passToGF)
4608        {
4609          CanonicalForm mipo= gf_mipo;
4610          setCharacteristic (p);
4611          Variable alpha= rootOf (mipo.mapinto());
4612          result= GF2FalphaRep (result, alpha);
4613        }
4614        if (k > 1)
4615        {
4616          result= GFMapDown (result, k);
4617          setCharacteristic (p, k, gf_name);
4618        }
4619        if (extOfExt)
4620          result= mapDown (result, primElem, imPrimElem, oldA, dest, source);
4621        return N (d*result);
4622      }
4623      DD[2] = DD[2] * ( b( lcDD[2] ) / lc( DD[2] ) );
4624      DD[1] = DD[1] * ( b( lcDD[1] ) / lc( DD[1] ) );
4625
4626      if (size (B*lcDD[2])/maxNumVars > sizePerVars1)
4627      {
4628        if (algExtension)
4629        {
4630          result= GCD_Fp_extension (F, G, a);
4631          if (extOfExt)
4632            result= mapDown (result, primElem, imPrimElem, oldA, dest, source);
4633          return N (d*result);
4634        }
4635        if (CFFactory::gettype() == GaloisFieldDomain)
4636        {
4637          result= GCD_GF (F, G);
4638          if (passToGF)
4639          {
4640            CanonicalForm mipo= gf_mipo;
4641            setCharacteristic (p);
4642            Variable alpha= rootOf (mipo.mapinto());
4643            result= GF2FalphaRep (result, alpha);
4644          }
4645          if (k > 1)
4646          {
4647            result= GFMapDown (result, k);
4648            setCharacteristic (p, k, gf_name);
4649          }
4650          return N (d*result);
4651        }
4652        else
4653          return N (d*GCD_small_p (F,G));
4654      }
4655
4656      TIMING_START (ez_p_hensel_lift);
4657      gcdfound= Hensel_P (B*lcD, DD, b, lcDD);
4658      TIMING_END_AND_PRINT (ez_p_hensel_lift, "time for Hensel lift in EZ_P: ");
4659
4660      if (gcdfound == -1) //things became dense
4661      {
4662        if (algExtension)
4663        {
4664          result= GCD_Fp_extension (F, G, a);
4665          if (extOfExt)
4666            result= mapDown (result, primElem, imPrimElem, oldA, dest, source);
4667          return N (d*result);
4668        }
4669        if (CFFactory::gettype() == GaloisFieldDomain)
4670        {
4671          result= GCD_GF (F, G);
4672          if (passToGF)
4673          {
4674            CanonicalForm mipo= gf_mipo;
4675            setCharacteristic (p);
4676            Variable alpha= rootOf (mipo.mapinto());
4677            result= GF2FalphaRep (result, alpha);
4678          }
4679          if (k > 1)
4680          {
4681            result= GFMapDown (result, k);
4682            setCharacteristic (p, k, gf_name);
4683          }
4684          return N (d*result);
4685        }
4686        else
4687        {
4688          if (p >= cf_getBigPrime(0))
4689            return N (d*sparseGCDFp (F,G));
4690          else
4691            return N (d*GCD_small_p (F,G));
4692        }
4693      }
4694
4695      if (gcdfound == 1)
4696      {
4697        TIMING_START (termination_test);
4698        contcand= content (DD[2], Variable (1));
4699        cand = DD[2] / contcand;
4700        if (B_is_F)
4701          gcdfound = fdivides( cand, G ) && cand*(DD[1]/(lcD/contcand)) == F;
4702        else
4703          gcdfound = fdivides( cand, F ) && cand*(DD[1]/(lcD/contcand)) == G;
4704        TIMING_END_AND_PRINT (termination_test,
4705                              "time for termination test EZ_P: ");
4706
4707        if (passToGF && gcdfound)
4708        {
4709          CanonicalForm mipo= gf_mipo;
4710          setCharacteristic (p);
4711          Variable alpha= rootOf (mipo.mapinto());
4712          cand= GF2FalphaRep (cand, alpha);
4713        }
4714        if (k > 1 && gcdfound)
4715        {
4716          cand= GFMapDown (cand, k);
4717          setCharacteristic (p, k, gf_name);
4718        }
4719        if (extOfExt && gcdfound)
4720          cand= mapDown (cand, primElem, imPrimElem, oldA, dest, source);
4721      }
4722    }
4723    delta--;
4724    goodPointCount= 0;
4725  }
4726  return N (d*cand);
4727}
4728
4729
4730#endif
Note: See TracBrowser for help on using the repository browser.