source: git/factory/cf_gcd_smallp.cc @ ffd883

fieker-DuValspielwiese
Last change on this file since ffd883 was babf57a, checked in by Martin Lee <martinlee84@…>, 11 years ago
fix: memory leaks
  • Property mode set to 100644
File size: 127.1 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  delete N;
1724  N= convertNmod_mat_t2FacCFMatrix (FLINTN);
1725  nmod_mat_clear (FLINTN);
1726#else
1727  int p= getCharacteristic ();
1728  if (fac_NTL_char != p)
1729  {
1730    fac_NTL_char= p;
1731    zz_p::init (p);
1732  }
1733  mat_zz_p *NTLN= convertFacCFMatrix2NTLmat_zz_p(*N);
1734  delete N;
1735  long rk= gauss (*NTLN);
1736
1737  N= convertNTLmat_zz_p2FacCFMatrix (*NTLN);
1738  delete NTLN;
1739#endif
1740
1741  L= CFArray (M.rows());
1742  for (int i= 0; i < M.rows(); i++)
1743    L[i]= (*N) (i + 1, M.columns() + 1);
1744  M= (*N) (1, M.rows(), 1, M.columns());
1745  delete N;
1746  return rk;
1747}
1748
1749long
1750gaussianElimFq (CFMatrix& M, CFArray& L, const Variable& alpha)
1751{
1752  ASSERT (L.size() <= M.rows(), "dimension exceeded");
1753  CFMatrix *N;
1754  N= new CFMatrix (M.rows(), M.columns() + 1);
1755
1756  for (int i= 1; i <= M.rows(); i++)
1757    for (int j= 1; j <= M.columns(); j++)
1758      (*N) (i, j)= M (i, j);
1759
1760  int j= 1;
1761  for (int i= 0; i < L.size(); i++, j++)
1762    (*N) (j, M.columns() + 1)= L[i];
1763  int p= getCharacteristic ();
1764  if (fac_NTL_char != p)
1765  {
1766    fac_NTL_char= p;
1767    zz_p::init (p);
1768  }
1769  zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
1770  zz_pE::init (NTLMipo);
1771  mat_zz_pE *NTLN= convertFacCFMatrix2NTLmat_zz_pE(*N);
1772  long rk= gauss (*NTLN);
1773
1774  delete N;
1775  N= convertNTLmat_zz_pE2FacCFMatrix (*NTLN, alpha);
1776
1777  delete NTLN;
1778
1779  M= (*N) (1, M.rows(), 1, M.columns());
1780  L= CFArray (M.rows());
1781  for (int i= 0; i < M.rows(); i++)
1782    L[i]= (*N) (i + 1, M.columns() + 1);
1783
1784  delete N;
1785  return rk;
1786}
1787
1788CFArray
1789solveSystemFp (const CFMatrix& M, const CFArray& L)
1790{
1791  ASSERT (L.size() <= M.rows(), "dimension exceeded");
1792  CFMatrix *N;
1793  N= new CFMatrix (M.rows(), M.columns() + 1);
1794
1795  for (int i= 1; i <= M.rows(); i++)
1796    for (int j= 1; j <= M.columns(); j++)
1797      (*N) (i, j)= M (i, j);
1798
1799  int j= 1;
1800  for (int i= 0; i < L.size(); i++, j++)
1801    (*N) (j, M.columns() + 1)= L[i];
1802
1803#ifdef HAVE_FLINT
1804  nmod_mat_t FLINTN;
1805  convertFacCFMatrix2nmod_mat_t (FLINTN, *N);
1806  long rk= nmod_mat_rref (FLINTN);
1807#else
1808  int p= getCharacteristic ();
1809  if (fac_NTL_char != p)
1810  {
1811    fac_NTL_char= p;
1812    zz_p::init (p);
1813  }
1814  mat_zz_p *NTLN= convertFacCFMatrix2NTLmat_zz_p(*N);
1815  long rk= gauss (*NTLN);
1816#endif
1817  delete N;
1818  if (rk != M.columns())
1819  {
1820#ifdef HAVE_FLINT
1821    nmod_mat_clear (FLINTN);
1822#else
1823    delete NTLN;
1824#endif
1825    return CFArray();
1826  }
1827#ifdef HAVE_FLINT
1828  N= convertNmod_mat_t2FacCFMatrix (FLINTN);
1829  nmod_mat_clear (FLINTN);
1830#else
1831  N= convertNTLmat_zz_p2FacCFMatrix (*NTLN);
1832  delete NTLN;
1833#endif
1834  CFArray A= readOffSolution (*N, rk);
1835
1836  delete N;
1837  return A;
1838}
1839
1840CFArray
1841solveSystemFq (const CFMatrix& M, const CFArray& L, const Variable& alpha)
1842{
1843  ASSERT (L.size() <= M.rows(), "dimension exceeded");
1844  CFMatrix *N;
1845  N= new CFMatrix (M.rows(), M.columns() + 1);
1846
1847  for (int i= 1; i <= M.rows(); i++)
1848    for (int j= 1; j <= M.columns(); j++)
1849      (*N) (i, j)= M (i, j);
1850  int j= 1;
1851  for (int i= 0; i < L.size(); i++, j++)
1852    (*N) (j, M.columns() + 1)= L[i];
1853  int p= getCharacteristic ();
1854  if (fac_NTL_char != p)
1855  {
1856    fac_NTL_char= p;
1857    zz_p::init (p);
1858  }
1859  zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
1860  zz_pE::init (NTLMipo);
1861  mat_zz_pE *NTLN= convertFacCFMatrix2NTLmat_zz_pE(*N);
1862  long rk= gauss (*NTLN);
1863
1864  delete N;
1865  if (rk != M.columns())
1866  {
1867    delete NTLN;
1868    return CFArray();
1869  }
1870  N= convertNTLmat_zz_pE2FacCFMatrix (*NTLN, alpha);
1871
1872  delete NTLN;
1873
1874  CFArray A= readOffSolution (*N, rk);
1875
1876  delete N;
1877  return A;
1878}
1879#endif
1880
1881CFArray
1882getMonoms (const CanonicalForm& F)
1883{
1884  if (F.inCoeffDomain())
1885  {
1886    CFArray result= CFArray (1);
1887    result [0]= 1;
1888    return result;
1889  }
1890  if (F.isUnivariate())
1891  {
1892    CFArray result= CFArray (size(F));
1893    int j= 0;
1894    for (CFIterator i= F; i.hasTerms(); i++, j++)
1895      result[j]= power (F.mvar(), i.exp());
1896    return result;
1897  }
1898  int numMon= size (F);
1899  CFArray result= CFArray (numMon);
1900  int j= 0;
1901  CFArray recResult;
1902  Variable x= F.mvar();
1903  CanonicalForm powX;
1904  for (CFIterator i= F; i.hasTerms(); i++)
1905  {
1906    powX= power (x, i.exp());
1907    recResult= getMonoms (i.coeff());
1908    for (int k= 0; k < recResult.size(); k++)
1909      result[j+k]= powX*recResult[k];
1910    j += recResult.size();
1911  }
1912  return result;
1913}
1914
1915#ifdef HAVE_NTL
1916CFArray
1917evaluateMonom (const CanonicalForm& F, const CFList& evalPoints)
1918{
1919  if (F.inCoeffDomain())
1920  {
1921    CFArray result= CFArray (1);
1922    result [0]= F;
1923    return result;
1924  }
1925  if (F.isUnivariate())
1926  {
1927    ASSERT (evalPoints.length() == 1,
1928            "expected an eval point with only one component");
1929    CFArray result= CFArray (size(F));
1930    int j= 0;
1931    CanonicalForm evalPoint= evalPoints.getLast();
1932    for (CFIterator i= F; i.hasTerms(); i++, j++)
1933      result[j]= power (evalPoint, i.exp());
1934    return result;
1935  }
1936  int numMon= size (F);
1937  CFArray result= CFArray (numMon);
1938  int j= 0;
1939  CanonicalForm evalPoint= evalPoints.getLast();
1940  CFList buf= evalPoints;
1941  buf.removeLast();
1942  CFArray recResult;
1943  CanonicalForm powEvalPoint;
1944  for (CFIterator i= F; i.hasTerms(); i++)
1945  {
1946    powEvalPoint= power (evalPoint, i.exp());
1947    recResult= evaluateMonom (i.coeff(), buf);
1948    for (int k= 0; k < recResult.size(); k++)
1949      result[j+k]= powEvalPoint*recResult[k];
1950    j += recResult.size();
1951  }
1952  return result;
1953}
1954
1955CFArray
1956evaluate (const CFArray& A, const CFList& evalPoints)
1957{
1958  CFArray result= A.size();
1959  CanonicalForm tmp;
1960  int k;
1961  for (int i= 0; i < A.size(); i++)
1962  {
1963    tmp= A[i];
1964    k= 1;
1965    for (CFListIterator j= evalPoints; j.hasItem(); j++, k++)
1966      tmp= tmp (j.getItem(), k);
1967    result[i]= tmp;
1968  }
1969  return result;
1970}
1971
1972CFList
1973evaluationPoints (const CanonicalForm& F, const CanonicalForm& G,
1974                  CanonicalForm& Feval, CanonicalForm& Geval,
1975                  const CanonicalForm& LCF, const bool& GF,
1976                  const Variable& alpha, bool& fail, CFList& list
1977                 )
1978{
1979  int k= tmax (F.level(), G.level()) - 1;
1980  Variable x= Variable (1);
1981  CFList result;
1982  FFRandom genFF;
1983  GFRandom genGF;
1984  int p= getCharacteristic ();
1985  double bound;
1986  if (alpha != Variable (1))
1987  {
1988    bound= pow ((double) p, (double) degree (getMipo(alpha)));
1989    bound= pow (bound, (double) k);
1990  }
1991  else if (GF)
1992  {
1993    bound= pow ((double) p, (double) getGFDegree());
1994    bound= pow ((double) bound, (double) k);
1995  }
1996  else
1997    bound= pow ((double) p, (double) k);
1998
1999  CanonicalForm random;
2000  int j;
2001  bool zeroOneOccured= false;
2002  bool allEqual= false;
2003  CanonicalForm buf;
2004  do
2005  {
2006    random= 0;
2007    // possible overflow if list.length() does not fit into a int
2008    if (list.length() >= bound)
2009    {
2010      fail= true;
2011      break;
2012    }
2013    for (int i= 0; i < k; i++)
2014    {
2015      if (GF)
2016      {
2017        result.append (genGF.generate());
2018        random += result.getLast()*power (x, i);
2019      }
2020      else if (alpha.level() != 1)
2021      {
2022        AlgExtRandomF genAlgExt (alpha);
2023        result.append (genAlgExt.generate());
2024        random += result.getLast()*power (x, i);
2025      }
2026      else
2027      {
2028        result.append (genFF.generate());
2029        random += result.getLast()*power (x, i);
2030      }
2031      if (result.getLast().isOne() || result.getLast().isZero())
2032        zeroOneOccured= true;
2033    }
2034    if (find (list, random))
2035    {
2036      zeroOneOccured= false;
2037      allEqual= false;
2038      result= CFList();
2039      continue;
2040    }
2041    if (zeroOneOccured)
2042    {
2043      list.append (random);
2044      zeroOneOccured= false;
2045      allEqual= false;
2046      result= CFList();
2047      continue;
2048    }
2049    // no zero at this point
2050    if (k > 1)
2051    {
2052      allEqual= true;
2053      CFIterator iter= random;
2054      buf= iter.coeff();
2055      iter++;
2056      for (; iter.hasTerms(); iter++)
2057        if (buf != iter.coeff())
2058          allEqual= false;
2059    }
2060    if (allEqual)
2061    {
2062      list.append (random);
2063      allEqual= false;
2064      zeroOneOccured= false;
2065      result= CFList();
2066      continue;
2067    }
2068
2069    Feval= F;
2070    Geval= G;
2071    CanonicalForm LCeval= LCF;
2072    j= 1;
2073    for (CFListIterator i= result; i.hasItem(); i++, j++)
2074    {
2075      Feval= Feval (i.getItem(), j);
2076      Geval= Geval (i.getItem(), j);
2077      LCeval= LCeval (i.getItem(), j);
2078    }
2079
2080    if (LCeval.isZero())
2081    {
2082      if (!find (list, random))
2083        list.append (random);
2084      zeroOneOccured= false;
2085      allEqual= false;
2086      result= CFList();
2087      continue;
2088    }
2089
2090    if (list.length() >= bound)
2091    {
2092      fail= true;
2093      break;
2094    }
2095  } while (find (list, random));
2096
2097  return result;
2098}
2099
2100/// multiply two lists componentwise
2101void mult (CFList& L1, const CFList& L2)
2102{
2103  ASSERT (L1.length() == L2.length(), "lists of the same size expected");
2104
2105  CFListIterator j= L2;
2106  for (CFListIterator i= L1; i.hasItem(); i++, j++)
2107    i.getItem() *= j.getItem();
2108}
2109
2110void eval (const CanonicalForm& A, const CanonicalForm& B, CanonicalForm& Aeval,
2111           CanonicalForm& Beval, const CFList& L)
2112{
2113  Aeval= A;
2114  Beval= B;
2115  int j= 1;
2116  for (CFListIterator i= L; i.hasItem(); i++, j++)
2117  {
2118    Aeval= Aeval (i.getItem(), j);
2119    Beval= Beval (i.getItem(), j);
2120  }
2121}
2122
2123CanonicalForm
2124monicSparseInterpol (const CanonicalForm& F, const CanonicalForm& G,
2125                     const CanonicalForm& skeleton, const Variable& alpha,
2126                     bool& fail, CFArray*& coeffMonoms, CFArray& Monoms
2127                    )
2128{
2129  CanonicalForm A= F;
2130  CanonicalForm B= G;
2131  if (F.isZero() && degree(G) > 0) return G/Lc(G);
2132  else if (G.isZero() && degree (F) > 0) return F/Lc(F);
2133  else if (F.isZero() && G.isZero()) return F.genOne();
2134  if (F.inBaseDomain() || G.inBaseDomain()) return F.genOne();
2135  if (F.isUnivariate() && fdivides(F, G)) return F/Lc(F);
2136  if (G.isUnivariate() && fdivides(G, F)) return G/Lc(G);
2137  if (F == G) return F/Lc(F);
2138
2139  ASSERT (degree (A, 1) == 0, "expected degree (F, 1) == 0");
2140  ASSERT (degree (B, 1) == 0, "expected degree (G, 1) == 0");
2141
2142  CFMap M,N;
2143  int best_level= myCompress (A, B, M, N, false);
2144
2145  if (best_level == 0)
2146    return B.genOne();
2147
2148  A= M(A);
2149  B= M(B);
2150
2151  Variable x= Variable (1);
2152
2153  //univariate case
2154  if (A.isUnivariate() && B.isUnivariate())
2155    return N (gcd (A, B));
2156
2157  CanonicalForm skel= M(skeleton);
2158  CanonicalForm cA, cB;    // content of A and B
2159  CanonicalForm ppA, ppB;    // primitive part of A and B
2160  CanonicalForm gcdcAcB;
2161  cA = uni_content (A);
2162  cB = uni_content (B);
2163  gcdcAcB= gcd (cA, cB);
2164  ppA= A/cA;
2165  ppB= B/cB;
2166
2167  CanonicalForm lcA, lcB;  // leading coefficients of A and B
2168  CanonicalForm gcdlcAlcB;
2169  lcA= uni_lcoeff (ppA);
2170  lcB= uni_lcoeff (ppB);
2171
2172  if (fdivides (lcA, lcB))
2173  {
2174    if (fdivides (A, B))
2175      return F/Lc(F);
2176  }
2177  if (fdivides (lcB, lcA))
2178  {
2179    if (fdivides (B, A))
2180      return G/Lc(G);
2181  }
2182
2183  gcdlcAlcB= gcd (lcA, lcB);
2184  int skelSize= size (skel, skel.mvar());
2185
2186  int j= 0;
2187  int biggestSize= 0;
2188
2189  for (CFIterator i= skel; i.hasTerms(); i++, j++)
2190    biggestSize= tmax (biggestSize, size (i.coeff()));
2191
2192  CanonicalForm g, Aeval, Beval;
2193
2194  CFList evalPoints;
2195  bool evalFail= false;
2196  CFList list;
2197  bool GF= false;
2198  CanonicalForm LCA= LC (A);
2199  CanonicalForm tmp;
2200  CFArray gcds= CFArray (biggestSize);
2201  CFList * pEvalPoints= new CFList [biggestSize];
2202  Variable V_buf= alpha;
2203  CFList source, dest;
2204  CanonicalForm prim_elem, im_prim_elem;
2205  for (int i= 0; i < biggestSize; i++)
2206  {
2207    if (i == 0)
2208      evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf, evalFail,
2209                                    list);
2210    else
2211    {
2212      mult (evalPoints, pEvalPoints [0]);
2213      eval (A, B, Aeval, Beval, evalPoints);
2214    }
2215
2216    if (evalFail)
2217    {
2218      if (V_buf.level() != 1)
2219      {
2220        do
2221        {
2222          Variable V_buf2= chooseExtension (V_buf);
2223          source= CFList();
2224          dest= CFList();
2225
2226          bool prim_fail= false;
2227          Variable V_buf3;
2228          prim_elem= primitiveElement (V_buf, V_buf3, prim_fail);
2229
2230          ASSERT (!prim_fail, "failure in integer factorizer");
2231          if (prim_fail)
2232            ; //ERROR
2233          else
2234            im_prim_elem= mapPrimElem (prim_elem, V_buf, V_buf2);
2235
2236          DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf));
2237          DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf2));
2238
2239          for (CFListIterator j= list; j.hasItem(); j++)
2240            j.getItem()= mapUp (j.getItem(), V_buf, V_buf2, prim_elem,
2241                                im_prim_elem, source, dest);
2242          for (int k= 0; k < i; k++)
2243          {
2244            for (CFListIterator j= pEvalPoints[k]; j.hasItem(); j++)
2245              j.getItem()= mapUp (j.getItem(), V_buf, V_buf2, prim_elem,
2246                                  im_prim_elem, source, dest);
2247            gcds[k]= mapUp (gcds[k], V_buf, V_buf2, prim_elem, im_prim_elem,
2248                            source, dest);
2249          }
2250
2251          if (alpha.level() != 1)
2252          {
2253            A= mapUp (A, V_buf, V_buf2, prim_elem, im_prim_elem, source,dest);
2254            B= mapUp (B, V_buf, V_buf2, prim_elem, im_prim_elem, source,dest);
2255          }
2256          evalFail= false;
2257          evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf,
2258                                        evalFail, list);
2259        } while (evalFail);
2260      }
2261      else
2262      {
2263        CanonicalForm mipo;
2264        int deg= 2;
2265        do {
2266          mipo= randomIrredpoly (deg, x);
2267          V_buf= rootOf (mipo);
2268          evalFail= false;
2269          evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf,
2270                                        evalFail, list);
2271          deg++;
2272        } while (evalFail);
2273      }
2274    }
2275
2276    g= gcd (Aeval, Beval);
2277    g /= Lc (g);
2278
2279    if (degree (g) != degree (skel) || (skelSize != size (g)))
2280    {
2281      delete[] pEvalPoints;
2282      fail= true;
2283      return 0;
2284    }
2285    CFIterator l= skel;
2286    for (CFIterator k= g; k.hasTerms(); k++, l++)
2287    {
2288      if (k.exp() != l.exp())
2289      {
2290        delete[] pEvalPoints;
2291        fail= true;
2292        return 0;
2293      }
2294    }
2295    pEvalPoints[i]= evalPoints;
2296    gcds[i]= g;
2297
2298    tmp= 0;
2299    int j= 0;
2300    for (CFListIterator k= evalPoints; k.hasItem(); k++, j++)
2301      tmp += k.getItem()*power (x, j);
2302    list.append (tmp);
2303  }
2304
2305  if (Monoms.size() == 0)
2306    Monoms= getMonoms (skel);
2307
2308  coeffMonoms= new CFArray [skelSize];
2309  j= 0;
2310  for (CFIterator i= skel; i.hasTerms(); i++, j++)
2311    coeffMonoms[j]= getMonoms (i.coeff());
2312
2313  CFArray* pL= new CFArray [skelSize];
2314  CFArray* pM= new CFArray [skelSize];
2315  for (int i= 0; i < biggestSize; i++)
2316  {
2317    CFIterator l= gcds [i];
2318    evalPoints= pEvalPoints [i];
2319    for (int k= 0; k < skelSize; k++, l++)
2320    {
2321      if (i == 0)
2322        pL[k]= CFArray (biggestSize);
2323      pL[k] [i]= l.coeff();
2324
2325      if (i == 0)
2326        pM[k]= evaluate (coeffMonoms [k], evalPoints);
2327    }
2328  }
2329
2330  CFArray solution;
2331  CanonicalForm result= 0;
2332  int ind= 0;
2333  CFArray bufArray;
2334  CFMatrix Mat;
2335  for (int k= 0; k < skelSize; k++)
2336  {
2337    if (biggestSize != coeffMonoms[k].size())
2338    {
2339      bufArray= CFArray (coeffMonoms[k].size());
2340      for (int i= 0; i < coeffMonoms[k].size(); i++)
2341        bufArray [i]= pL[k] [i];
2342      solution= solveGeneralVandermonde (pM [k], bufArray);
2343    }
2344    else
2345      solution= solveGeneralVandermonde (pM [k], pL[k]);
2346
2347    if (solution.size() == 0)
2348    {
2349      delete[] pEvalPoints;
2350      delete[] pM;
2351      delete[] pL;
2352      delete[] coeffMonoms;
2353      fail= true;
2354      return 0;
2355    }
2356    for (int l= 0; l < solution.size(); l++)
2357      result += solution[l]*Monoms [ind + l];
2358    ind += solution.size();
2359  }
2360
2361  delete[] pEvalPoints;
2362  delete[] pM;
2363  delete[] pL;
2364  delete[] coeffMonoms;
2365
2366  if (alpha.level() != 1 && V_buf != alpha)
2367  {
2368    CFList u, v;
2369    result= mapDown (result, prim_elem, im_prim_elem, alpha, u, v);
2370  }
2371
2372  result= N(result);
2373  if (fdivides (result, F) && fdivides (result, G))
2374    return result;
2375  else
2376  {
2377    fail= true;
2378    return 0;
2379  }
2380}
2381
2382CanonicalForm
2383nonMonicSparseInterpol (const CanonicalForm& F, const CanonicalForm& G,
2384                        const CanonicalForm& skeleton, const Variable& alpha,
2385                        bool& fail, CFArray*& coeffMonoms, CFArray& Monoms
2386                       )
2387{
2388  CanonicalForm A= F;
2389  CanonicalForm B= G;
2390  if (F.isZero() && degree(G) > 0) return G/Lc(G);
2391  else if (G.isZero() && degree (F) > 0) return F/Lc(F);
2392  else if (F.isZero() && G.isZero()) return F.genOne();
2393  if (F.inBaseDomain() || G.inBaseDomain()) return F.genOne();
2394  if (F.isUnivariate() && fdivides(F, G)) return F/Lc(F);
2395  if (G.isUnivariate() && fdivides(G, F)) return G/Lc(G);
2396  if (F == G) return F/Lc(F);
2397
2398  ASSERT (degree (A, 1) == 0, "expected degree (F, 1) == 0");
2399  ASSERT (degree (B, 1) == 0, "expected degree (G, 1) == 0");
2400
2401  CFMap M,N;
2402  int best_level= myCompress (A, B, M, N, false);
2403
2404  if (best_level == 0)
2405    return B.genOne();
2406
2407  A= M(A);
2408  B= M(B);
2409
2410  Variable x= Variable (1);
2411
2412  //univariate case
2413  if (A.isUnivariate() && B.isUnivariate())
2414    return N (gcd (A, B));
2415
2416  CanonicalForm skel= M(skeleton);
2417
2418  CanonicalForm cA, cB;    // content of A and B
2419  CanonicalForm ppA, ppB;    // primitive part of A and B
2420  CanonicalForm gcdcAcB;
2421  cA = uni_content (A);
2422  cB = uni_content (B);
2423  gcdcAcB= gcd (cA, cB);
2424  ppA= A/cA;
2425  ppB= B/cB;
2426
2427  CanonicalForm lcA, lcB;  // leading coefficients of A and B
2428  CanonicalForm gcdlcAlcB;
2429  lcA= uni_lcoeff (ppA);
2430  lcB= uni_lcoeff (ppB);
2431
2432  if (fdivides (lcA, lcB))
2433  {
2434    if (fdivides (A, B))
2435      return F/Lc(F);
2436  }
2437  if (fdivides (lcB, lcA))
2438  {
2439    if (fdivides (B, A))
2440      return G/Lc(G);
2441  }
2442
2443  gcdlcAlcB= gcd (lcA, lcB);
2444  int skelSize= size (skel, skel.mvar());
2445
2446  int j= 0;
2447  int biggestSize= 0;
2448  int bufSize;
2449  int numberUni= 0;
2450  for (CFIterator i= skel; i.hasTerms(); i++, j++)
2451  {
2452    bufSize= size (i.coeff());
2453    biggestSize= tmax (biggestSize, bufSize);
2454    numberUni += bufSize;
2455  }
2456  numberUni--;
2457  numberUni= (int) ceil ( (double) numberUni / (skelSize - 1));
2458  biggestSize= tmax (biggestSize , numberUni);
2459
2460  numberUni= biggestSize + size (LC(skel)) - 1;
2461  int biggestSize2= tmax (numberUni, biggestSize);
2462
2463  CanonicalForm g, Aeval, Beval;
2464
2465  CFList evalPoints;
2466  CFArray coeffEval;
2467  bool evalFail= false;
2468  CFList list;
2469  bool GF= false;
2470  CanonicalForm LCA= LC (A);
2471  CanonicalForm tmp;
2472  CFArray gcds= CFArray (biggestSize);
2473  CFList * pEvalPoints= new CFList [biggestSize];
2474  Variable V_buf= alpha;
2475  CFList source, dest;
2476  CanonicalForm prim_elem, im_prim_elem;
2477  for (int i= 0; i < biggestSize; i++)
2478  {
2479    if (i == 0)
2480    {
2481      if (getCharacteristic() > 3)
2482        evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf,
2483                                    evalFail, list);
2484      else
2485        evalFail= true;
2486
2487      if (evalFail)
2488      {
2489        if (V_buf.level() != 1)
2490        {
2491          do
2492          {
2493            Variable V_buf2= chooseExtension (V_buf);
2494            source= CFList();
2495            dest= CFList();
2496
2497            bool prim_fail= false;
2498            Variable V_buf3;
2499            prim_elem= primitiveElement (V_buf, V_buf3, prim_fail);
2500
2501            ASSERT (!prim_fail, "failure in integer factorizer");
2502            if (prim_fail)
2503              ; //ERROR
2504            else
2505              im_prim_elem= mapPrimElem (prim_elem, V_buf, V_buf2);
2506
2507            DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf));
2508            DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf2));
2509
2510            for (CFListIterator i= list; i.hasItem(); i++)
2511              i.getItem()= mapUp (i.getItem(), V_buf, V_buf2, prim_elem,
2512                                im_prim_elem, source, dest);
2513            evalFail= false;
2514            evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf,
2515                                          evalFail, list);
2516          } while (evalFail);
2517        }
2518        else
2519        {
2520          CanonicalForm mipo;
2521          int deg= 2;
2522          do {
2523            mipo= randomIrredpoly (deg, x);
2524            V_buf= rootOf (mipo);
2525            evalFail= false;
2526            evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf,
2527                                          evalFail, list);
2528            deg++;
2529          } while (evalFail);
2530        }
2531      }
2532    }
2533    else
2534    {
2535      mult (evalPoints, pEvalPoints[0]);
2536      eval (A, B, Aeval, Beval, evalPoints);
2537    }
2538
2539    g= gcd (Aeval, Beval);
2540    g /= Lc (g);
2541
2542    if (degree (g) != degree (skel) || (skelSize != size (g)))
2543    {
2544      delete[] pEvalPoints;
2545      fail= true;
2546      return 0;
2547    }
2548    CFIterator l= skel;
2549    for (CFIterator k= g; k.hasTerms(); k++, l++)
2550    {
2551      if (k.exp() != l.exp())
2552      {
2553        delete[] pEvalPoints;
2554        fail= true;
2555        return 0;
2556      }
2557    }
2558    pEvalPoints[i]= evalPoints;
2559    gcds[i]= g;
2560
2561    tmp= 0;
2562    int j= 0;
2563    for (CFListIterator k= evalPoints; k.hasItem(); k++, j++)
2564      tmp += k.getItem()*power (x, j);
2565    list.append (tmp);
2566  }
2567
2568  if (Monoms.size() == 0)
2569    Monoms= getMonoms (skel);
2570
2571  coeffMonoms= new CFArray [skelSize];
2572
2573  j= 0;
2574  for (CFIterator i= skel; i.hasTerms(); i++, j++)
2575    coeffMonoms[j]= getMonoms (i.coeff());
2576
2577  int minimalColumnsIndex;
2578  if (skelSize > 1)
2579    minimalColumnsIndex= 1;
2580  else
2581    minimalColumnsIndex= 0;
2582  int minimalColumns=-1;
2583
2584  CFArray* pM= new CFArray [skelSize];
2585  CFMatrix Mat;
2586  // find the Matrix with minimal number of columns
2587  for (int i= 0; i < skelSize; i++)
2588  {
2589    pM[i]= CFArray (coeffMonoms[i].size());
2590    if (i == 1)
2591      minimalColumns= coeffMonoms[i].size();
2592    if (i > 1)
2593    {
2594      if (minimalColumns > coeffMonoms[i].size())
2595      {
2596        minimalColumns= coeffMonoms[i].size();
2597        minimalColumnsIndex= i;
2598      }
2599    }
2600  }
2601  CFMatrix* pMat= new CFMatrix [2];
2602  pMat[0]= CFMatrix (biggestSize, coeffMonoms[0].size());
2603  pMat[1]= CFMatrix (biggestSize, coeffMonoms[minimalColumnsIndex].size());
2604  CFArray* pL= new CFArray [skelSize];
2605  for (int i= 0; i < biggestSize; i++)
2606  {
2607    CFIterator l= gcds [i];
2608    evalPoints= pEvalPoints [i];
2609    for (int k= 0; k < skelSize; k++, l++)
2610    {
2611      if (i == 0)
2612        pL[k]= CFArray (biggestSize);
2613      pL[k] [i]= l.coeff();
2614
2615      if (i == 0 && k != 0 && k != minimalColumnsIndex)
2616        pM[k]= evaluate (coeffMonoms[k], evalPoints);
2617      else if (i == 0 && (k == 0 || k == minimalColumnsIndex))
2618        coeffEval= evaluate (coeffMonoms[k], evalPoints);
2619      if (i > 0 && (k == 0 || k == minimalColumnsIndex))
2620        coeffEval= evaluate (coeffMonoms[k], evalPoints);
2621
2622      if (k == 0)
2623      {
2624        if (pMat[k].rows() >= i + 1)
2625        {
2626          for (int ind= 1; ind < coeffEval.size() + 1; ind++)
2627            pMat[k] (i + 1, ind)= coeffEval[ind - 1];
2628        }
2629      }
2630      if (k == minimalColumnsIndex)
2631      {
2632        if (pMat[1].rows() >= i + 1)
2633        {
2634          for (int ind= 1; ind < coeffEval.size() + 1; ind++)
2635            pMat[1] (i + 1, ind)= coeffEval[ind - 1];
2636        }
2637      }
2638    }
2639  }
2640
2641  CFArray solution;
2642  CanonicalForm result= 0;
2643  int ind= 1;
2644  int matRows, matColumns;
2645  matRows= pMat[1].rows();
2646  matColumns= pMat[0].columns() - 1;
2647  matColumns += pMat[1].columns();
2648
2649  Mat= CFMatrix (matRows, matColumns);
2650  for (int i= 1; i <= matRows; i++)
2651    for (int j= 1; j <= pMat[1].columns(); j++)
2652      Mat (i, j)= pMat[1] (i, j);
2653
2654  for (int j= pMat[1].columns() + 1; j <= matColumns;
2655       j++, ind++)
2656  {
2657    for (int i= 1; i <= matRows; i++)
2658      Mat (i, j)= (-pMat [0] (i, ind + 1))*pL[minimalColumnsIndex] [i - 1];
2659  }
2660
2661  CFArray firstColumn= CFArray (pMat[0].rows());
2662  for (int i= 0; i < pMat[0].rows(); i++)
2663    firstColumn [i]= pMat[0] (i + 1, 1);
2664
2665  CFArray bufArray= pL[minimalColumnsIndex];
2666
2667  for (int i= 0; i < biggestSize; i++)
2668    pL[minimalColumnsIndex] [i] *= firstColumn[i];
2669
2670  CFMatrix bufMat= pMat[1];
2671  pMat[1]= Mat;
2672
2673  if (V_buf.level() != 1)
2674    solution= solveSystemFq (pMat[1],
2675                             pL[minimalColumnsIndex], V_buf);
2676  else
2677    solution= solveSystemFp (pMat[1],
2678                             pL[minimalColumnsIndex]);
2679
2680  if (solution.size() == 0)
2681  { //Javadi's method failed, try deKleine, Monagan, Wittkopf instead
2682    CFMatrix bufMat0= pMat[0];
2683    delete [] pMat;
2684    pMat= new CFMatrix [skelSize];
2685    pL[minimalColumnsIndex]= bufArray;
2686    CFList* bufpEvalPoints= NULL;
2687    CFArray bufGcds;
2688    if (biggestSize != biggestSize2)
2689    {
2690      bufpEvalPoints= pEvalPoints;
2691      pEvalPoints= new CFList [biggestSize2];
2692      bufGcds= gcds;
2693      gcds= CFArray (biggestSize2);
2694      for (int i= 0; i < biggestSize; i++)
2695      {
2696        pEvalPoints[i]= bufpEvalPoints [i];
2697        gcds[i]= bufGcds[i];
2698      }
2699      for (int i= 0; i < biggestSize2 - biggestSize; i++)
2700      {
2701        mult (evalPoints, pEvalPoints[0]);
2702        eval (A, B, Aeval, Beval, evalPoints);
2703        g= gcd (Aeval, Beval);
2704        g /= Lc (g);
2705        if (degree (g) != degree (skel) || (skelSize != size (g)))
2706        {
2707          delete[] pEvalPoints;
2708          delete[] pMat;
2709          delete[] pL;
2710          delete[] coeffMonoms;
2711          delete[] pM;
2712          if (bufpEvalPoints != NULL)
2713            delete [] bufpEvalPoints;
2714          fail= true;
2715          return 0;
2716        }
2717        CFIterator l= skel;
2718        for (CFIterator k= g; k.hasTerms(); k++, l++)
2719        {
2720          if (k.exp() != l.exp())
2721          {
2722            delete[] pEvalPoints;
2723            delete[] pMat;
2724            delete[] pL;
2725            delete[] coeffMonoms;
2726            delete[] pM;
2727            if (bufpEvalPoints != NULL)
2728              delete [] bufpEvalPoints;
2729            fail= true;
2730            return 0;
2731          }
2732        }
2733        pEvalPoints[i + biggestSize]= evalPoints;
2734        gcds[i + biggestSize]= g;
2735      }
2736    }
2737    for (int i= 0; i < biggestSize; i++)
2738    {
2739      CFIterator l= gcds [i];
2740      evalPoints= pEvalPoints [i];
2741      for (int k= 1; k < skelSize; k++, l++)
2742      {
2743        if (i == 0)
2744          pMat[k]= CFMatrix (biggestSize2,coeffMonoms[k].size()+biggestSize2-1);
2745        if (k == minimalColumnsIndex)
2746          continue;
2747        coeffEval= evaluate (coeffMonoms[k], evalPoints);
2748        if (pMat[k].rows() >= i + 1)
2749        {
2750          for (int ind= 1; ind < coeffEval.size() + 1; ind++)
2751            pMat[k] (i + 1, ind)= coeffEval[ind - 1];
2752        }
2753      }
2754    }
2755    Mat= bufMat0;
2756    pMat[0]= CFMatrix (biggestSize2, coeffMonoms[0].size() + biggestSize2 - 1);
2757    for (int j= 1; j <= Mat.rows(); j++)
2758      for (int k= 1; k <= Mat.columns(); k++)
2759        pMat [0] (j,k)= Mat (j,k);
2760    Mat= bufMat;
2761    for (int j= 1; j <= Mat.rows(); j++)
2762      for (int k= 1; k <= Mat.columns(); k++)
2763        pMat [minimalColumnsIndex] (j,k)= Mat (j,k);
2764    // write old matrix entries into new matrices
2765    for (int i= 0; i < skelSize; i++)
2766    {
2767      bufArray= pL[i];
2768      pL[i]= CFArray (biggestSize2);
2769      for (int j= 0; j < bufArray.size(); j++)
2770        pL[i] [j]= bufArray [j];
2771    }
2772    //write old vector entries into new and add new entries to old matrices
2773    for (int i= 0; i < biggestSize2 - biggestSize; i++)
2774    {
2775      CFIterator l= gcds [i + biggestSize];
2776      evalPoints= pEvalPoints [i + biggestSize];
2777      for (int k= 0; k < skelSize; k++, l++)
2778      {
2779        pL[k] [i + biggestSize]= l.coeff();
2780        coeffEval= evaluate (coeffMonoms[k], evalPoints);
2781        if (pMat[k].rows() >= i + biggestSize + 1)
2782        {
2783          for (int ind= 1; ind < coeffEval.size() + 1; ind++)
2784            pMat[k] (i + biggestSize + 1, ind)= coeffEval[ind - 1];
2785        }
2786      }
2787    }
2788    // begin new
2789    for (int i= 0; i < skelSize; i++)
2790    {
2791      if (pL[i].size() > 1)
2792      {
2793        for (int j= 2; j <= pMat[i].rows(); j++)
2794          pMat[i] (j, coeffMonoms[i].size() + j - 1)=
2795              -pL[i] [j - 1];
2796      }
2797    }
2798
2799    matColumns= biggestSize2 - 1;
2800    matRows= 0;
2801    for (int i= 0; i < skelSize; i++)
2802    {
2803      if (V_buf.level() == 1)
2804        (void) gaussianElimFp (pMat[i], pL[i]);
2805      else
2806        (void) gaussianElimFq (pMat[i], pL[i], V_buf);
2807
2808      if (pMat[i] (coeffMonoms[i].size(), coeffMonoms[i].size()) == 0)
2809      {
2810        delete[] pEvalPoints;
2811        delete[] pMat;
2812        delete[] pL;
2813        delete[] coeffMonoms;
2814        delete[] pM;
2815        if (bufpEvalPoints != NULL)
2816          delete [] bufpEvalPoints;
2817        fail= true;
2818        return 0;
2819      }
2820      matRows += pMat[i].rows() - coeffMonoms[i].size();
2821    }
2822
2823    CFMatrix bufMat;
2824    Mat= CFMatrix (matRows, matColumns);
2825    ind= 0;
2826    bufArray= CFArray (matRows);
2827    CFArray bufArray2;
2828    for (int i= 0; i < skelSize; i++)
2829    {
2830      bufMat= pMat[i] (coeffMonoms[i].size() + 1, pMat[i].rows(),
2831                       coeffMonoms[i].size() + 1, pMat[i].columns());
2832
2833      for (int j= 1; j <= bufMat.rows(); j++)
2834        for (int k= 1; k <= bufMat.columns(); k++)
2835          Mat (j + ind, k)= bufMat(j, k);
2836      bufArray2= coeffMonoms[i].size();
2837      for (int j= 1; j <= pMat[i].rows(); j++)
2838      {
2839        if (j > coeffMonoms[i].size())
2840          bufArray [j-coeffMonoms[i].size() + ind - 1]= pL[i] [j - 1];
2841        else
2842          bufArray2 [j - 1]= pL[i] [j - 1];
2843      }
2844      pL[i]= bufArray2;
2845      ind += bufMat.rows();
2846      pMat[i]= pMat[i] (1, coeffMonoms[i].size(), 1, pMat[i].columns());
2847    }
2848
2849    if (V_buf.level() != 1)
2850      solution= solveSystemFq (Mat, bufArray, V_buf);
2851    else
2852      solution= solveSystemFp (Mat, bufArray);
2853
2854    if (solution.size() == 0)
2855    {
2856      delete[] pEvalPoints;
2857      delete[] pMat;
2858      delete[] pL;
2859      delete[] coeffMonoms;
2860      delete[] pM;
2861      if (bufpEvalPoints != NULL)
2862        delete [] bufpEvalPoints;
2863      fail= true;
2864      return 0;
2865    }
2866
2867    ind= 0;
2868    result= 0;
2869    CFArray bufSolution;
2870    for (int i= 0; i < skelSize; i++)
2871    {
2872      bufSolution= readOffSolution (pMat[i], pL[i], solution);
2873      for (int i= 0; i < bufSolution.size(); i++, ind++)
2874        result += Monoms [ind]*bufSolution[i];
2875    }
2876    if (alpha.level() != 1 && V_buf != alpha)
2877    {
2878      CFList u, v;
2879      result= mapDown (result, prim_elem, im_prim_elem, alpha, u, v);
2880    }
2881    result= N(result);
2882    delete[] pEvalPoints;
2883    delete[] pMat;
2884    delete[] pL;
2885    delete[] coeffMonoms;
2886    delete[] pM;
2887
2888    if (bufpEvalPoints != NULL)
2889      delete [] bufpEvalPoints;
2890    if (fdivides (result, F) && fdivides (result, G))
2891      return result;
2892    else
2893    {
2894      fail= true;
2895      return 0;
2896    }
2897  } // end of deKleine, Monagan & Wittkopf
2898
2899  result += Monoms[0];
2900  int ind2= 0, ind3= 2;
2901  ind= 0;
2902  for (int l= 0; l < minimalColumnsIndex; l++)
2903    ind += coeffMonoms[l].size();
2904  for (int l= solution.size() - pMat[0].columns() + 1; l < solution.size();
2905       l++, ind2++, ind3++)
2906  {
2907    result += solution[l]*Monoms [1 + ind2];
2908    for (int i= 0; i < pMat[0].rows(); i++)
2909      firstColumn[i] += solution[l]*pMat[0] (i + 1, ind3);
2910  }
2911  for (int l= 0; l < solution.size() + 1 - pMat[0].columns(); l++)
2912    result += solution[l]*Monoms [ind + l];
2913  ind= coeffMonoms[0].size();
2914  for (int k= 1; k < skelSize; k++)
2915  {
2916    if (k == minimalColumnsIndex)
2917    {
2918      ind += coeffMonoms[k].size();
2919      continue;
2920    }
2921    if (k != minimalColumnsIndex)
2922    {
2923      for (int i= 0; i < biggestSize; i++)
2924        pL[k] [i] *= firstColumn [i];
2925    }
2926
2927    if (biggestSize != coeffMonoms[k].size() && k != minimalColumnsIndex)
2928    {
2929      bufArray= CFArray (coeffMonoms[k].size());
2930      for (int i= 0; i < bufArray.size(); i++)
2931        bufArray [i]= pL[k] [i];
2932      solution= solveGeneralVandermonde (pM[k], bufArray);
2933    }
2934    else
2935      solution= solveGeneralVandermonde (pM[k], pL[k]);
2936
2937    if (solution.size() == 0)
2938    {
2939      delete[] pEvalPoints;
2940      delete[] pMat;
2941      delete[] pL;
2942      delete[] coeffMonoms;
2943      delete[] pM;
2944      fail= true;
2945      return 0;
2946    }
2947    if (k != minimalColumnsIndex)
2948    {
2949      for (int l= 0; l < solution.size(); l++)
2950        result += solution[l]*Monoms [ind + l];
2951      ind += solution.size();
2952    }
2953  }
2954
2955  delete[] pEvalPoints;
2956  delete[] pMat;
2957  delete[] pL;
2958  delete[] pM;
2959  delete[] coeffMonoms;
2960
2961  if (alpha.level() != 1 && V_buf != alpha)
2962  {
2963    CFList u, v;
2964    result= mapDown (result, prim_elem, im_prim_elem, alpha, u, v);
2965  }
2966  result= N(result);
2967
2968  if (fdivides (result, F) && fdivides (result, G))
2969    return result;
2970  else
2971  {
2972    fail= true;
2973    return 0;
2974  }
2975}
2976
2977CanonicalForm sparseGCDFq (const CanonicalForm& F, const CanonicalForm& G,
2978                           const Variable & alpha, CFList& l, bool& topLevel)
2979{
2980  CanonicalForm A= F;
2981  CanonicalForm B= G;
2982  if (F.isZero() && degree(G) > 0) return G/Lc(G);
2983  else if (G.isZero() && degree (F) > 0) return F/Lc(F);
2984  else if (F.isZero() && G.isZero()) return F.genOne();
2985  if (F.inBaseDomain() || G.inBaseDomain()) return F.genOne();
2986  if (F.isUnivariate() && fdivides(F, G)) return F/Lc(F);
2987  if (G.isUnivariate() && fdivides(G, F)) return G/Lc(G);
2988  if (F == G) return F/Lc(F);
2989
2990  CFMap M,N;
2991  int best_level= myCompress (A, B, M, N, topLevel);
2992
2993  if (best_level == 0) return B.genOne();
2994
2995  A= M(A);
2996  B= M(B);
2997
2998  Variable x= Variable (1);
2999
3000  //univariate case
3001  if (A.isUnivariate() && B.isUnivariate())
3002    return N (gcd (A, B));
3003
3004  CanonicalForm cA, cB;    // content of A and B
3005  CanonicalForm ppA, ppB;    // primitive part of A and B
3006  CanonicalForm gcdcAcB;
3007
3008  cA = uni_content (A);
3009  cB = uni_content (B);
3010  gcdcAcB= gcd (cA, cB);
3011  ppA= A/cA;
3012  ppB= B/cB;
3013
3014  CanonicalForm lcA, lcB;  // leading coefficients of A and B
3015  CanonicalForm gcdlcAlcB;
3016  lcA= uni_lcoeff (ppA);
3017  lcB= uni_lcoeff (ppB);
3018
3019  if (fdivides (lcA, lcB))
3020  {
3021    if (fdivides (A, B))
3022      return F/Lc(F);
3023  }
3024  if (fdivides (lcB, lcA))
3025  {
3026    if (fdivides (B, A))
3027      return G/Lc(G);
3028  }
3029
3030  gcdlcAlcB= gcd (lcA, lcB);
3031
3032  int d= totaldegree (ppA, Variable (2), Variable (ppA.level()));
3033  int d0;
3034
3035  if (d == 0)
3036    return N(gcdcAcB);
3037  d0= totaldegree (ppB, Variable (2), Variable (ppB.level()));
3038
3039  if (d0 < d)
3040    d= d0;
3041
3042  if (d == 0)
3043    return N(gcdcAcB);
3044
3045  CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH, skeleton;
3046  CanonicalForm newtonPoly= 1;
3047  m= gcdlcAlcB;
3048  G_m= 0;
3049  H= 0;
3050  bool fail= false;
3051  topLevel= false;
3052  bool inextension= false;
3053  Variable V_buf= alpha;
3054  CanonicalForm prim_elem, im_prim_elem;
3055  CFList source, dest;
3056  do // first do
3057  {
3058    random_element= randomElement (m, V_buf, l, fail);
3059    if (random_element == 0 && !fail)
3060    {
3061      if (!find (l, random_element))
3062        l.append (random_element);
3063      continue;
3064    }
3065    if (fail)
3066    {
3067      source= CFList();
3068      dest= CFList();
3069
3070      Variable V_buf3= V_buf;
3071      V_buf= chooseExtension (V_buf);
3072      bool prim_fail= false;
3073      Variable V_buf2;
3074      prim_elem= primitiveElement (alpha, V_buf2, prim_fail);
3075
3076      if (V_buf3 != alpha)
3077      {
3078        m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
3079        G_m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
3080        newtonPoly= mapDown (newtonPoly, prim_elem, im_prim_elem, alpha,
3081                             source, dest);
3082        ppA= mapDown (ppA, prim_elem, im_prim_elem, alpha, source, dest);
3083        ppB= mapDown (ppB, prim_elem, im_prim_elem, alpha, source, dest);
3084        gcdlcAlcB= mapDown (gcdlcAlcB, prim_elem, im_prim_elem, alpha, source,
3085                            dest);
3086        for (CFListIterator i= l; i.hasItem(); i++)
3087          i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, alpha,
3088                                source, dest);
3089      }
3090
3091      ASSERT (!prim_fail, "failure in integer factorizer");
3092      if (prim_fail)
3093        ; //ERROR
3094      else
3095        im_prim_elem= mapPrimElem (prim_elem, alpha, V_buf);
3096
3097      DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (alpha));
3098      DEBOUTLN (cerr, "getMipo (V_buf2)= " << getMipo (V_buf2));
3099      inextension= true;
3100      for (CFListIterator i= l; i.hasItem(); i++)
3101        i.getItem()= mapUp (i.getItem(), alpha, V_buf, prim_elem,
3102                             im_prim_elem, source, dest);
3103      m= mapUp (m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3104      G_m= mapUp (G_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3105      newtonPoly= mapUp (newtonPoly, alpha, V_buf, prim_elem, im_prim_elem,
3106                          source, dest);
3107      ppA= mapUp (ppA, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3108      ppB= mapUp (ppB, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3109      gcdlcAlcB= mapUp (gcdlcAlcB, alpha, V_buf, prim_elem, im_prim_elem,
3110                         source, dest);
3111
3112      fail= false;
3113      random_element= randomElement (m, V_buf, l, fail );
3114      DEBOUTLN (cerr, "fail= " << fail);
3115      CFList list;
3116      TIMING_START (gcd_recursion);
3117      G_random_element=
3118      sparseGCDFq (ppA (random_element, x), ppB (random_element, x), V_buf,
3119                        list, topLevel);
3120      TIMING_END_AND_PRINT (gcd_recursion,
3121                            "time for recursive call: ");
3122      DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3123    }
3124    else
3125    {
3126      CFList list;
3127      TIMING_START (gcd_recursion);
3128      G_random_element=
3129      sparseGCDFq (ppA(random_element, x), ppB(random_element, x), V_buf,
3130                        list, topLevel);
3131      TIMING_END_AND_PRINT (gcd_recursion,
3132                            "time for recursive call: ");
3133      DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3134    }
3135
3136    if (!G_random_element.inCoeffDomain())
3137      d0= totaldegree (G_random_element, Variable(2),
3138                       Variable (G_random_element.level()));
3139    else
3140      d0= 0;
3141
3142    if (d0 == 0)
3143      return N(gcdcAcB);
3144    if (d0 >  d)
3145    {
3146      if (!find (l, random_element))
3147        l.append (random_element);
3148      continue;
3149    }
3150
3151    G_random_element=
3152    (gcdlcAlcB(random_element, x)/uni_lcoeff (G_random_element))
3153    * G_random_element;
3154
3155    skeleton= G_random_element;
3156    if (!G_random_element.inCoeffDomain())
3157      d0= totaldegree (G_random_element, Variable(2),
3158                       Variable (G_random_element.level()));
3159    else
3160      d0= 0;
3161
3162    if (d0 <  d)
3163    {
3164      m= gcdlcAlcB;
3165      newtonPoly= 1;
3166      G_m= 0;
3167      d= d0;
3168    }
3169
3170    H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
3171    if (uni_lcoeff (H) == gcdlcAlcB)
3172    {
3173      cH= uni_content (H);
3174      ppH= H/cH;
3175      if (inextension)
3176      {
3177        CFList u, v;
3178        //maybe it's better to test if ppH is an element of F(\alpha) before
3179        //mapping down
3180        if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
3181        {
3182          DEBOUTLN (cerr, "ppH before mapDown= " << ppH);
3183          ppH= mapDown (ppH, prim_elem, im_prim_elem, alpha, u, v);
3184          ppH /= Lc(ppH);
3185          DEBOUTLN (cerr, "ppH after mapDown= " << ppH);
3186          return N(gcdcAcB*ppH);
3187        }
3188      }
3189      else if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
3190        return N(gcdcAcB*ppH);
3191    }
3192    G_m= H;
3193    newtonPoly= newtonPoly*(x - random_element);
3194    m= m*(x - random_element);
3195    if (!find (l, random_element))
3196      l.append (random_element);
3197
3198    if (getCharacteristic () > 3 && size (skeleton) < 100)
3199    {
3200      CFArray Monoms;
3201      CFArray *coeffMonoms;
3202      do //second do
3203      {
3204        random_element= randomElement (m, V_buf, l, fail);
3205        if (random_element == 0 && !fail)
3206        {
3207          if (!find (l, random_element))
3208            l.append (random_element);
3209          continue;
3210        }
3211        if (fail)
3212        {
3213          source= CFList();
3214          dest= CFList();
3215
3216          Variable V_buf3= V_buf;
3217          V_buf= chooseExtension (V_buf);
3218          bool prim_fail= false;
3219          Variable V_buf2;
3220          prim_elem= primitiveElement (alpha, V_buf2, prim_fail);
3221
3222          if (V_buf3 != alpha)
3223          {
3224            m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
3225            G_m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
3226            newtonPoly= mapDown (newtonPoly, prim_elem, im_prim_elem, alpha,
3227                                 source, dest);
3228            ppA= mapDown (ppA, prim_elem, im_prim_elem, alpha, source, dest);
3229            ppB= mapDown (ppB, prim_elem, im_prim_elem, alpha, source, dest);
3230            gcdlcAlcB= mapDown (gcdlcAlcB, prim_elem, im_prim_elem, alpha,
3231                                source, dest);
3232            for (CFListIterator i= l; i.hasItem(); i++)
3233              i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, alpha,
3234                                    source, dest);
3235          }
3236
3237          ASSERT (!prim_fail, "failure in integer factorizer");
3238          if (prim_fail)
3239            ; //ERROR
3240          else
3241            im_prim_elem= mapPrimElem (prim_elem, alpha, V_buf);
3242
3243          DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (alpha));
3244          DEBOUTLN (cerr, "getMipo (V_buf2)= " << getMipo (V_buf2));
3245          inextension= true;
3246          for (CFListIterator i= l; i.hasItem(); i++)
3247            i.getItem()= mapUp (i.getItem(), alpha, V_buf, prim_elem,
3248                                im_prim_elem, source, dest);
3249          m= mapUp (m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3250          G_m= mapUp (G_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3251          newtonPoly= mapUp (newtonPoly, alpha, V_buf, prim_elem, im_prim_elem,
3252                              source, dest);
3253          ppA= mapUp (ppA, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3254          ppB= mapUp (ppB, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3255
3256          gcdlcAlcB= mapUp (gcdlcAlcB, alpha, V_buf, prim_elem, im_prim_elem,
3257                            source, dest);
3258
3259          fail= false;
3260          random_element= randomElement (m, V_buf, l, fail);
3261          DEBOUTLN (cerr, "fail= " << fail);
3262          CFList list;
3263          TIMING_START (gcd_recursion);
3264
3265          //sparseInterpolation
3266          bool sparseFail= false;
3267          if (LC (skeleton).inCoeffDomain())
3268            G_random_element=
3269            monicSparseInterpol (ppA (random_element, x), ppB(random_element,x),
3270                                skeleton,V_buf, sparseFail, coeffMonoms,Monoms);
3271          else
3272            G_random_element=
3273            nonMonicSparseInterpol (ppA(random_element,x),ppB(random_element,x),
3274                                    skeleton, V_buf, sparseFail, coeffMonoms,
3275                                    Monoms);
3276          if (sparseFail)
3277            break;
3278
3279          TIMING_END_AND_PRINT (gcd_recursion,
3280                                "time for recursive call: ");
3281          DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3282        }
3283        else
3284        {
3285          CFList list;
3286          TIMING_START (gcd_recursion);
3287          bool sparseFail= false;
3288          if (LC (skeleton).inCoeffDomain())
3289            G_random_element=
3290            monicSparseInterpol (ppA (random_element, x),ppB(random_element, x),
3291                                skeleton,V_buf, sparseFail,coeffMonoms, Monoms);
3292          else
3293            G_random_element=
3294            nonMonicSparseInterpol (ppA(random_element,x),ppB(random_element,x),
3295                                    skeleton, V_buf, sparseFail, coeffMonoms,
3296                                    Monoms);
3297          if (sparseFail)
3298            break;
3299
3300          TIMING_END_AND_PRINT (gcd_recursion,
3301                                "time for recursive call: ");
3302          DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3303        }
3304
3305        if (!G_random_element.inCoeffDomain())
3306          d0= totaldegree (G_random_element, Variable(2),
3307                           Variable (G_random_element.level()));
3308        else
3309          d0= 0;
3310
3311        if (d0 == 0)
3312          return N(gcdcAcB);
3313        if (d0 >  d)
3314        {
3315          if (!find (l, random_element))
3316            l.append (random_element);
3317          continue;
3318        }
3319
3320        G_random_element=
3321        (gcdlcAlcB(random_element, x)/uni_lcoeff (G_random_element))
3322        * G_random_element;
3323
3324        if (!G_random_element.inCoeffDomain())
3325          d0= totaldegree (G_random_element, Variable(2),
3326                          Variable (G_random_element.level()));
3327        else
3328          d0= 0;
3329
3330        if (d0 <  d)
3331        {
3332          m= gcdlcAlcB;
3333          newtonPoly= 1;
3334          G_m= 0;
3335          d= d0;
3336        }
3337
3338        TIMING_START (newton_interpolation);
3339        H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
3340        TIMING_END_AND_PRINT (newton_interpolation,
3341                              "time for newton interpolation: ");
3342
3343        //termination test
3344        if (uni_lcoeff (H) == gcdlcAlcB)
3345        {
3346          cH= uni_content (H);
3347          ppH= H/cH;
3348          if (inextension)
3349          {
3350            CFList u, v;
3351            //maybe it's better to test if ppH is an element of F(\alpha) before
3352            //mapping down
3353            if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
3354            {
3355              DEBOUTLN (cerr, "ppH before mapDown= " << ppH);
3356              ppH= mapDown (ppH, prim_elem, im_prim_elem, alpha, u, v);
3357              ppH /= Lc(ppH);
3358              DEBOUTLN (cerr, "ppH after mapDown= " << ppH);
3359              return N(gcdcAcB*ppH);
3360            }
3361          }
3362          else if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
3363          {
3364            return N(gcdcAcB*ppH);
3365          }
3366        }
3367
3368        G_m= H;
3369        newtonPoly= newtonPoly*(x - random_element);
3370        m= m*(x - random_element);
3371        if (!find (l, random_element))
3372          l.append (random_element);
3373
3374      } while (1);
3375    }
3376  } while (1);
3377}
3378
3379CanonicalForm sparseGCDFp (const CanonicalForm& F, const CanonicalForm& G,
3380                           bool& topLevel, CFList& l)
3381{
3382  CanonicalForm A= F;
3383  CanonicalForm B= G;
3384  if (F.isZero() && degree(G) > 0) return G/Lc(G);
3385  else if (G.isZero() && degree (F) > 0) return F/Lc(F);
3386  else if (F.isZero() && G.isZero()) return F.genOne();
3387  if (F.inBaseDomain() || G.inBaseDomain()) return F.genOne();
3388  if (F.isUnivariate() && fdivides(F, G)) return F/Lc(F);
3389  if (G.isUnivariate() && fdivides(G, F)) return G/Lc(G);
3390  if (F == G) return F/Lc(F);
3391
3392  CFMap M,N;
3393  int best_level= myCompress (A, B, M, N, topLevel);
3394
3395  if (best_level == 0) return B.genOne();
3396
3397  A= M(A);
3398  B= M(B);
3399
3400  Variable x= Variable (1);
3401
3402  //univariate case
3403  if (A.isUnivariate() && B.isUnivariate())
3404    return N (gcd (A, B));
3405
3406  CanonicalForm cA, cB;    // content of A and B
3407  CanonicalForm ppA, ppB;    // primitive part of A and B
3408  CanonicalForm gcdcAcB;
3409
3410  cA = uni_content (A);
3411  cB = uni_content (B);
3412  gcdcAcB= gcd (cA, cB);
3413  ppA= A/cA;
3414  ppB= B/cB;
3415
3416  CanonicalForm lcA, lcB;  // leading coefficients of A and B
3417  CanonicalForm gcdlcAlcB;
3418  lcA= uni_lcoeff (ppA);
3419  lcB= uni_lcoeff (ppB);
3420
3421  if (fdivides (lcA, lcB))
3422  {
3423    if (fdivides (A, B))
3424      return F/Lc(F);
3425  }
3426  if (fdivides (lcB, lcA))
3427  {
3428    if (fdivides (B, A))
3429      return G/Lc(G);
3430  }
3431
3432  gcdlcAlcB= gcd (lcA, lcB);
3433
3434  int d= totaldegree (ppA, Variable (2), Variable (ppA.level()));
3435  int d0;
3436
3437  if (d == 0)
3438    return N(gcdcAcB);
3439
3440  d0= totaldegree (ppB, Variable (2), Variable (ppB.level()));
3441
3442  if (d0 < d)
3443    d= d0;
3444
3445  if (d == 0)
3446    return N(gcdcAcB);
3447
3448  CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH, skeleton;
3449  CanonicalForm newtonPoly= 1;
3450  m= gcdlcAlcB;
3451  G_m= 0;
3452  H= 0;
3453  bool fail= false;
3454  topLevel= false;
3455  bool inextension= false;
3456  Variable V_buf, alpha;
3457  CanonicalForm prim_elem, im_prim_elem;
3458  CFList source, dest;
3459  do //first do
3460  {
3461    if (inextension)
3462      random_element= randomElement (m, V_buf, l, fail);
3463    else
3464      random_element= FpRandomElement (m, l, fail);
3465    if (random_element == 0 && !fail)
3466    {
3467      if (!find (l, random_element))
3468        l.append (random_element);
3469      continue;
3470    }
3471
3472    if (!fail && !inextension)
3473    {
3474      CFList list;
3475      TIMING_START (gcd_recursion);
3476      G_random_element=
3477      sparseGCDFp (ppA (random_element,x), ppB (random_element,x), topLevel,
3478                   list);
3479      TIMING_END_AND_PRINT (gcd_recursion,
3480                            "time for recursive call: ");
3481      DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3482    }
3483    else if (!fail && inextension)
3484    {
3485      CFList list;
3486      TIMING_START (gcd_recursion);
3487      G_random_element=
3488      sparseGCDFq (ppA (random_element, x), ppB (random_element, x), alpha,
3489                        list, topLevel);
3490      TIMING_END_AND_PRINT (gcd_recursion,
3491                            "time for recursive call: ");
3492      DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3493    }
3494    else if (fail && !inextension)
3495    {
3496      source= CFList();
3497      dest= CFList();
3498      CFList list;
3499      CanonicalForm mipo;
3500      int deg= 2;
3501      do
3502      {
3503        mipo= randomIrredpoly (deg, x);
3504        alpha= rootOf (mipo);
3505        inextension= true;
3506        fail= false;
3507        random_element= randomElement (m, alpha, l, fail);
3508        deg++;
3509      } while (fail);
3510      V_buf= alpha;
3511      list= CFList();
3512      TIMING_START (gcd_recursion);
3513      G_random_element=
3514      sparseGCDFq (ppA (random_element, x), ppB (random_element, x), alpha,
3515                        list, topLevel);
3516      TIMING_END_AND_PRINT (gcd_recursion,
3517                            "time for recursive call: ");
3518      DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3519    }
3520    else if (fail && inextension)
3521    {
3522      source= CFList();
3523      dest= CFList();
3524
3525      Variable V_buf3= V_buf;
3526      V_buf= chooseExtension (V_buf);
3527      bool prim_fail= false;
3528      Variable V_buf2;
3529      CanonicalForm prim_elem, im_prim_elem;
3530      prim_elem= primitiveElement (alpha, V_buf2, prim_fail);
3531
3532      if (V_buf3 != alpha)
3533      {
3534        m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
3535        G_m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
3536        newtonPoly= mapDown (newtonPoly, prim_elem, im_prim_elem, alpha, source,
3537                             dest);
3538        ppA= mapDown (ppA, prim_elem, im_prim_elem, alpha, source, dest);
3539        ppB= mapDown (ppB, prim_elem, im_prim_elem, alpha, source, dest);
3540        gcdlcAlcB= mapDown (gcdlcAlcB, prim_elem, im_prim_elem, alpha, source,
3541                            dest);
3542        for (CFListIterator i= l; i.hasItem(); i++)
3543          i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, alpha,
3544                                source, dest);
3545      }
3546
3547      ASSERT (!prim_fail, "failure in integer factorizer");
3548      if (prim_fail)
3549        ; //ERROR
3550      else
3551        im_prim_elem= mapPrimElem (prim_elem, alpha, V_buf);
3552
3553      DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (alpha));
3554      DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf2));
3555
3556      for (CFListIterator i= l; i.hasItem(); i++)
3557        i.getItem()= mapUp (i.getItem(), alpha, V_buf, prim_elem,
3558                             im_prim_elem, source, dest);
3559      m= mapUp (m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3560      G_m= mapUp (G_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3561      newtonPoly= mapUp (newtonPoly, alpha, V_buf, prim_elem, im_prim_elem,
3562                          source, dest);
3563      ppA= mapUp (ppA, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3564      ppB= mapUp (ppB, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3565      gcdlcAlcB= mapUp (gcdlcAlcB, alpha, V_buf, prim_elem, im_prim_elem,
3566                         source, dest);
3567      fail= false;
3568      random_element= randomElement (m, V_buf, l, fail );
3569      DEBOUTLN (cerr, "fail= " << fail);
3570      CFList list;
3571      TIMING_START (gcd_recursion);
3572      G_random_element=
3573      sparseGCDFq (ppA (random_element, x), ppB (random_element, x), V_buf,
3574                        list, topLevel);
3575      TIMING_END_AND_PRINT (gcd_recursion,
3576                            "time for recursive call: ");
3577      DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3578    }
3579
3580    if (!G_random_element.inCoeffDomain())
3581      d0= totaldegree (G_random_element, Variable(2),
3582                       Variable (G_random_element.level()));
3583    else
3584      d0= 0;
3585
3586    if (d0 == 0)
3587      return N(gcdcAcB);
3588    if (d0 >  d)
3589    {
3590      if (!find (l, random_element))
3591        l.append (random_element);
3592      continue;
3593    }
3594
3595    G_random_element=
3596    (gcdlcAlcB(random_element, x)/uni_lcoeff (G_random_element))
3597    * G_random_element;
3598
3599    skeleton= G_random_element;
3600
3601    if (!G_random_element.inCoeffDomain())
3602      d0= totaldegree (G_random_element, Variable(2),
3603                       Variable (G_random_element.level()));
3604    else
3605      d0= 0;
3606
3607    if (d0 <  d)
3608    {
3609      m= gcdlcAlcB;
3610      newtonPoly= 1;
3611      G_m= 0;
3612      d= d0;
3613    }
3614
3615    H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
3616
3617    if (uni_lcoeff (H) == gcdlcAlcB)
3618    {
3619      cH= uni_content (H);
3620      ppH= H/cH;
3621      ppH /= Lc (ppH);
3622      DEBOUTLN (cerr, "ppH= " << ppH);
3623
3624      if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
3625        return N(gcdcAcB*ppH);
3626    }
3627    G_m= H;
3628    newtonPoly= newtonPoly*(x - random_element);
3629    m= m*(x - random_element);
3630    if (!find (l, random_element))
3631      l.append (random_element);
3632
3633    if ((getCharacteristic() > 3 && size (skeleton) < 200))
3634    {
3635      CFArray Monoms;
3636      CFArray* coeffMonoms;
3637
3638      do //second do
3639      {
3640        if (inextension)
3641          random_element= randomElement (m, alpha, l, fail);
3642        else
3643          random_element= FpRandomElement (m, l, fail);
3644        if (random_element == 0 && !fail)
3645        {
3646          if (!find (l, random_element))
3647            l.append (random_element);
3648          continue;
3649        }
3650
3651        bool sparseFail= false;
3652        if (!fail && !inextension)
3653        {
3654          CFList list;
3655          TIMING_START (gcd_recursion);
3656          if (LC (skeleton).inCoeffDomain())
3657            G_random_element=
3658            monicSparseInterpol(ppA(random_element, x), ppB (random_element, x),
3659                                skeleton, Variable(1), sparseFail, coeffMonoms,
3660                                Monoms);
3661          else
3662            G_random_element=
3663            nonMonicSparseInterpol(ppA(random_element,x), ppB(random_element,x),
3664                                    skeleton, Variable (1), sparseFail,
3665                                    coeffMonoms, Monoms);
3666          TIMING_END_AND_PRINT (gcd_recursion,
3667                                "time for recursive call: ");
3668          DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3669        }
3670        else if (!fail && inextension)
3671        {
3672          CFList list;
3673          TIMING_START (gcd_recursion);
3674          if (LC (skeleton).inCoeffDomain())
3675            G_random_element=
3676            monicSparseInterpol(ppA (random_element,x), ppB (random_element, x),
3677                                skeleton, alpha, sparseFail, coeffMonoms,
3678                                Monoms);
3679          else
3680            G_random_element=
3681            nonMonicSparseInterpol(ppA(random_element,x), ppB(random_element,x),
3682                                   skeleton, alpha, sparseFail, coeffMonoms,
3683                                   Monoms);
3684          TIMING_END_AND_PRINT (gcd_recursion,
3685                                "time for recursive call: ");
3686          DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3687        }
3688        else if (fail && !inextension)
3689        {
3690          source= CFList();
3691          dest= CFList();
3692          CFList list;
3693          CanonicalForm mipo;
3694          int deg= 2;
3695          do
3696          {
3697            mipo= randomIrredpoly (deg, x);
3698            alpha= rootOf (mipo);
3699            inextension= true;
3700            fail= false;
3701            random_element= randomElement (m, alpha, l, fail);
3702            deg++;
3703          } while (fail);
3704          V_buf= alpha;
3705          list= CFList();
3706          TIMING_START (gcd_recursion);
3707          if (LC (skeleton).inCoeffDomain())
3708            G_random_element=
3709            monicSparseInterpol (ppA (random_element,x), ppB (random_element,x),
3710                                skeleton, alpha, sparseFail, coeffMonoms,
3711                                Monoms);
3712          else
3713            G_random_element=
3714            nonMonicSparseInterpol(ppA(random_element,x), ppB(random_element,x),
3715                                   skeleton, alpha, sparseFail, coeffMonoms,
3716                                   Monoms);
3717          TIMING_END_AND_PRINT (gcd_recursion,
3718                                "time for recursive call: ");
3719          DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3720        }
3721        else if (fail && inextension)
3722        {
3723          source= CFList();
3724          dest= CFList();
3725
3726          Variable V_buf3= V_buf;
3727          V_buf= chooseExtension (V_buf);
3728          bool prim_fail= false;
3729          Variable V_buf2;
3730          CanonicalForm prim_elem, im_prim_elem;
3731          prim_elem= primitiveElement (alpha, V_buf2, prim_fail);
3732
3733          if (V_buf3 != alpha)
3734          {
3735            m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
3736            G_m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
3737            newtonPoly= mapDown (newtonPoly, prim_elem, im_prim_elem, alpha,
3738                                 source, dest);
3739            ppA= mapDown (ppA, prim_elem, im_prim_elem, alpha, source, dest);
3740            ppB= mapDown (ppB, prim_elem, im_prim_elem, alpha, source, dest);
3741            gcdlcAlcB= mapDown (gcdlcAlcB, prim_elem, im_prim_elem, alpha,
3742                                source, dest);
3743            for (CFListIterator i= l; i.hasItem(); i++)
3744              i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, alpha,
3745                                    source, dest);
3746          }
3747
3748          ASSERT (!prim_fail, "failure in integer factorizer");
3749          if (prim_fail)
3750            ; //ERROR
3751          else
3752            im_prim_elem= mapPrimElem (prim_elem, alpha, V_buf);
3753
3754          DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (alpha));
3755          DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf2));
3756
3757          for (CFListIterator i= l; i.hasItem(); i++)
3758            i.getItem()= mapUp (i.getItem(), alpha, V_buf, prim_elem,
3759                                im_prim_elem, source, dest);
3760          m= mapUp (m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3761          G_m= mapUp (G_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3762          newtonPoly= mapUp (newtonPoly, alpha, V_buf, prim_elem, im_prim_elem,
3763                              source, dest);
3764          ppA= mapUp (ppA, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3765          ppB= mapUp (ppB, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
3766          gcdlcAlcB= mapUp (gcdlcAlcB, alpha, V_buf, prim_elem, im_prim_elem,
3767                            source, dest);
3768          fail= false;
3769          random_element= randomElement (m, V_buf, l, fail );
3770          DEBOUTLN (cerr, "fail= " << fail);
3771          CFList list;
3772          TIMING_START (gcd_recursion);
3773          if (LC (skeleton).inCoeffDomain())
3774            G_random_element=
3775            monicSparseInterpol (ppA (random_element, x), ppB (random_element, x),
3776                                skeleton, V_buf, sparseFail, coeffMonoms,
3777                                Monoms);
3778          else
3779            G_random_element=
3780            nonMonicSparseInterpol (ppA(random_element,x), ppB(random_element,x),
3781                                    skeleton, V_buf, sparseFail,
3782                                    coeffMonoms, Monoms);
3783          TIMING_END_AND_PRINT (gcd_recursion,
3784                                "time for recursive call: ");
3785          DEBOUTLN (cerr, "G_random_element= " << G_random_element);
3786        }
3787
3788        if (sparseFail)
3789          break;
3790
3791        if (!G_random_element.inCoeffDomain())
3792          d0= totaldegree (G_random_element, Variable(2),
3793                           Variable (G_random_element.level()));
3794        else
3795          d0= 0;
3796
3797        if (d0 == 0)
3798          return N(gcdcAcB);
3799        if (d0 >  d)
3800        {
3801          if (!find (l, random_element))
3802            l.append (random_element);
3803          continue;
3804        }
3805
3806        G_random_element=
3807        (gcdlcAlcB(random_element, x)/uni_lcoeff (G_random_element))
3808        * G_random_element;
3809
3810        if (!G_random_element.inCoeffDomain())
3811          d0= totaldegree (G_random_element, Variable(2),
3812                           Variable (G_random_element.level()));
3813        else
3814          d0= 0;
3815
3816        if (d0 <  d)
3817        {
3818          m= gcdlcAlcB;
3819          newtonPoly= 1;
3820          G_m= 0;
3821          d= d0;
3822        }
3823
3824        TIMING_START (newton_interpolation);
3825        H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
3826        TIMING_END_AND_PRINT (newton_interpolation,
3827                              "time for newton interpolation: ");
3828
3829        //termination test
3830        if (uni_lcoeff (H) == gcdlcAlcB)
3831        {
3832          cH= uni_content (H);
3833          ppH= H/cH;
3834          ppH /= Lc (ppH);
3835          DEBOUTLN (cerr, "ppH= " << ppH);
3836          if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
3837            return N(gcdcAcB*ppH);
3838        }
3839
3840        G_m= H;
3841        newtonPoly= newtonPoly*(x - random_element);
3842        m= m*(x - random_element);
3843        if (!find (l, random_element))
3844          l.append (random_element);
3845
3846      } while (1); //end of second do
3847    }
3848    else
3849      return N(gcdcAcB*GCD_small_p (ppA, ppB));
3850  } while (1); //end of first do
3851}
3852
3853static inline
3854int compress4EZGCD (const CanonicalForm& F, const CanonicalForm& G, CFMap & M,
3855                    CFMap & N, int& both_non_zero)
3856{
3857  int n= tmax (F.level(), G.level());
3858  int * degsf= new int [n + 1];
3859  int * degsg= new int [n + 1];
3860
3861  for (int i = 0; i <= n; i++)
3862    degsf[i]= degsg[i]= 0;
3863
3864  degsf= degrees (F, degsf);
3865  degsg= degrees (G, degsg);
3866
3867  both_non_zero= 0;
3868  int f_zero= 0;
3869  int g_zero= 0;
3870
3871  for (int i= 1; i <= n; i++)
3872  {
3873    if (degsf[i] != 0 && degsg[i] != 0)
3874    {
3875      both_non_zero++;
3876      continue;
3877    }
3878    if (degsf[i] == 0 && degsg[i] != 0 && i <= G.level())
3879    {
3880      f_zero++;
3881      continue;
3882    }
3883    if (degsg[i] == 0 && degsf[i] && i <= F.level())
3884    {
3885      g_zero++;
3886      continue;
3887    }
3888  }
3889
3890  if (both_non_zero == 0)
3891  {
3892    delete [] degsf;
3893    delete [] degsg;
3894    return 0;
3895  }
3896
3897  // map Variables which do not occur in both polynomials to higher levels
3898  int k= 1;
3899  int l= 1;
3900  for (int i= 1; i <= n; i++)
3901  {
3902    if (degsf[i] != 0 && degsg[i] == 0 && i <= F.level())
3903    {
3904      if (k + both_non_zero != i)
3905      {
3906        M.newpair (Variable (i), Variable (k + both_non_zero));
3907        N.newpair (Variable (k + both_non_zero), Variable (i));
3908      }
3909      k++;
3910    }
3911    if (degsf[i] == 0 && degsg[i] != 0 && i <= G.level())
3912    {
3913      if (l + g_zero + both_non_zero != i)
3914      {
3915        M.newpair (Variable (i), Variable (l + g_zero + both_non_zero));
3916        N.newpair (Variable (l + g_zero + both_non_zero), Variable (i));
3917      }
3918      l++;
3919    }
3920  }
3921
3922  // sort Variables x_{i} in decreasing order of
3923  // min(deg_{x_{i}}(f),deg_{x_{i}}(g))
3924  int m= tmin (F.level(), G.level());
3925  int max_min_deg;
3926  k= both_non_zero;
3927  l= 0;
3928  int i= 1;
3929  while (k > 0)
3930  {
3931    max_min_deg= tmin (degsf[i], degsg[i]);
3932    while (max_min_deg == 0)
3933    {
3934      i++;
3935      max_min_deg= tmin (degsf[i], degsg[i]);
3936    }
3937    for (int j= i + 1; j <= m; j++)
3938    {
3939      if ((tmin (degsf[j],degsg[j]) < max_min_deg) &&
3940          (tmin (degsf[j], degsg[j]) != 0))
3941      {
3942        max_min_deg= tmin (degsf[j], degsg[j]);
3943        l= j;
3944      }
3945    }
3946
3947    if (l != 0)
3948    {
3949      if (l != k)
3950      {
3951        M.newpair (Variable (l), Variable(k));
3952        N.newpair (Variable (k), Variable(l));
3953        degsf[l]= 0;
3954        degsg[l]= 0;
3955        l= 0;
3956      }
3957      else
3958      {
3959        degsf[l]= 0;
3960        degsg[l]= 0;
3961        l= 0;
3962      }
3963    }
3964    else if (l == 0)
3965    {
3966      if (i != k)
3967      {
3968        M.newpair (Variable (i), Variable (k));
3969        N.newpair (Variable (k), Variable (i));
3970        degsf[i]= 0;
3971        degsg[i]= 0;
3972      }
3973      else
3974      {
3975        degsf[i]= 0;
3976        degsg[i]= 0;
3977      }
3978      i++;
3979    }
3980    k--;
3981  }
3982
3983  delete [] degsf;
3984  delete [] degsg;
3985
3986  return both_non_zero;
3987}
3988
3989static inline
3990CanonicalForm myShift2Zero (const CanonicalForm& F, CFList& Feval,
3991                            const CFList& evaluation)
3992{
3993  CanonicalForm A= F;
3994  int k= 2;
3995  for (CFListIterator i= evaluation; i.hasItem(); i++, k++)
3996    A= A (Variable (k) + i.getItem(), k);
3997
3998  CanonicalForm buf= A;
3999  Feval= CFList();
4000  Feval.append (buf);
4001  for (k= evaluation.length() + 1; k > 2; k--)
4002  {
4003    buf= mod (buf, Variable (k));
4004    Feval.insert (buf);
4005  }
4006  return A;
4007}
4008
4009static inline
4010CanonicalForm myReverseShift (const CanonicalForm& F, const CFList& evaluation)
4011{
4012  int l= evaluation.length() + 1;
4013  CanonicalForm result= F;
4014  CFListIterator j= evaluation;
4015  for (int i= 2; i < l + 1; i++, j++)
4016  {
4017    if (F.level() < i)
4018      continue;
4019    result= result (Variable (i) - j.getItem(), i);
4020  }
4021  return result;
4022}
4023
4024static inline
4025Evaluation optimize4Lift (const CanonicalForm& F, CFMap & M,
4026                    CFMap & N, const Evaluation& A)
4027{
4028  int n= F.level();
4029  int * degsf= new int [n + 1];
4030
4031  for (int i = 0; i <= n; i++)
4032    degsf[i]= 0;
4033
4034  degsf= degrees (F, degsf);
4035
4036  Evaluation result= Evaluation (A.min(), A.max());
4037  ASSERT (A.min() == 2, "expected A.min() == 2");
4038  int max_deg;
4039  int k= n;
4040  int l= 1;
4041  int i= 2;
4042  int pos= 2;
4043  while (k > 1)
4044  {
4045    max_deg= degsf [i];
4046    while (max_deg == 0)
4047    {
4048      i++;
4049      max_deg= degsf [i];
4050    }
4051    l= i;
4052    for (int j= i + 1; j <=  n; j++)
4053    {
4054      if (degsf[j] > max_deg)
4055      {
4056        max_deg= degsf[j];
4057        l= j;
4058      }
4059    }
4060
4061    if (l <= n)
4062    {
4063      if (l != pos)
4064      {
4065        result.setValue (pos, A [l]);
4066        M.newpair (Variable (l), Variable (pos));
4067        N.newpair (Variable (pos), Variable (l));
4068        degsf[l]= 0;
4069        l= 2;
4070        if (k == 2 && n == 3)
4071        {
4072          result.setValue (l, A [pos]);
4073          M.newpair (Variable (pos), Variable (l));
4074          N.newpair (Variable (l), Variable (pos));
4075          degsf[pos]= 0;
4076        }
4077      }
4078      else
4079      {
4080        result.setValue (l, A [l]);
4081        degsf [l]= 0;
4082      }
4083    }
4084    pos++;
4085    k--;
4086    l= 2;
4087  }
4088
4089  delete [] degsf;
4090
4091  return result;
4092}
4093
4094static inline
4095int Hensel_P (const CanonicalForm & UU, CFArray & G, const Evaluation & AA,
4096              const CFArray& LeadCoeffs )
4097{
4098  CFList factors;
4099  factors.append (G[1]);
4100  factors.append (G[2]);
4101
4102  CFMap NN, MM;
4103  Evaluation A= optimize4Lift (UU, MM, NN, AA);
4104
4105  CanonicalForm U= MM (UU);
4106  CFArray LCs= CFArray (1,2);
4107  LCs [1]= MM (LeadCoeffs [1]);
4108  LCs [2]= MM (LeadCoeffs [2]);
4109
4110  CFList evaluation;
4111  long termEstimate= size (U);
4112  for (int i= A.min(); i <= A.max(); i++)
4113  {
4114    if (!A[i].isZero() && (getCharacteristic() > degree (U,i))) //TODO find a good estimate for getCharacteristic() <= degree (U,i)
4115    {
4116      termEstimate *= degree (U,i)*2;
4117      termEstimate /= 3;
4118    }
4119    evaluation.append (A [i]);
4120  }
4121  if (termEstimate/getNumVars(U) > 500)
4122    return -1;
4123  CFList UEval;
4124  CanonicalForm shiftedU= myShift2Zero (U, UEval, evaluation);
4125
4126  if (size (shiftedU)/getNumVars (U) > 500)
4127    return -1;
4128
4129  CFArray shiftedLCs= CFArray (2);
4130  CFList shiftedLCsEval1, shiftedLCsEval2;
4131  shiftedLCs[0]= myShift2Zero (LCs[1], shiftedLCsEval1, evaluation);
4132  shiftedLCs[1]= myShift2Zero (LCs[2], shiftedLCsEval2, evaluation);
4133  factors.insert (1);
4134  int liftBound= degree (UEval.getLast(), 2) + 1;
4135  CFArray Pi;
4136  CFMatrix M= CFMatrix (liftBound, factors.length() - 1);
4137  CFList diophant;
4138  CFArray lcs= CFArray (2);
4139  lcs [0]= shiftedLCsEval1.getFirst();
4140  lcs [1]= shiftedLCsEval2.getFirst();
4141  nonMonicHenselLift12 (UEval.getFirst(), factors, liftBound, Pi, diophant, M,
4142                        lcs, false);
4143
4144  for (CFListIterator i= factors; i.hasItem(); i++)
4145  {
4146    if (!fdivides (i.getItem(), UEval.getFirst()))
4147      return 0;
4148  }
4149
4150  int * liftBounds;
4151  bool noOneToOne= false;
4152  if (U.level() > 2)
4153  {
4154    liftBounds= new int [U.level() - 1]; /* index: 0.. U.level()-2 */
4155    liftBounds[0]= liftBound;
4156    for (int i= 1; i < U.level() - 1; i++)
4157      liftBounds[i]= degree (shiftedU, Variable (i + 2)) + 1;
4158    factors= nonMonicHenselLift2 (UEval, factors, liftBounds, U.level() - 1,
4159                                  false, shiftedLCsEval1, shiftedLCsEval2, Pi,
4160                                  diophant, noOneToOne);
4161    delete [] liftBounds;
4162    if (noOneToOne)
4163      return 0;
4164  }
4165  G[1]= factors.getFirst();
4166  G[2]= factors.getLast();
4167  G[1]= myReverseShift (G[1], evaluation);
4168  G[2]= myReverseShift (G[2], evaluation);
4169  G[1]= NN (G[1]);
4170  G[2]= NN (G[2]);
4171  return 1;
4172}
4173
4174static inline
4175bool findeval_P (const CanonicalForm & F, const CanonicalForm & G,
4176                 CanonicalForm & Fb, CanonicalForm & Gb, CanonicalForm & Db,
4177                 REvaluation & b, int delta, int degF, int degG, int maxeval,
4178                 int & count, int& k, int bound, int& l)
4179{
4180  if( count == 0 && delta != 0)
4181  {
4182    if( count++ > maxeval )
4183      return false;
4184  }
4185  if (count > 0)
4186  {
4187    b.nextpoint(k);
4188    if (k == 0)
4189      k++;
4190    l++;
4191    if (l > bound)
4192    {
4193      l= 1;
4194      k++;
4195      if (k > tmax (F.level(), G.level()) - 1)
4196        return false;
4197      b.nextpoint (k);
4198    }
4199    if (count++ > maxeval)
4200      return false;
4201  }
4202  while( true )
4203  {
4204    Fb = b( F );
4205    if( degree( Fb, 1 ) == degF )
4206    {
4207      Gb = b( G );
4208      if( degree( Gb, 1 ) == degG )
4209      {
4210        Db = gcd( Fb, Gb );
4211        if( delta > 0 )
4212        {
4213          if( degree( Db, 1 ) <= delta )
4214            return true;
4215        }
4216        else
4217          return true;
4218      }
4219    }
4220    if (k == 0)
4221      k++;
4222    b.nextpoint(k);
4223    l++;
4224    if (l > bound)
4225    {
4226      l= 1;
4227      k++;
4228      if (k > tmax (F.level(), G.level()) - 1)
4229        return false;
4230      b.nextpoint (k);
4231    }
4232    if( count++ > maxeval )
4233      return false;
4234  }
4235}
4236
4237// parameters for heuristic
4238static int maxNumEval= 200;
4239static int sizePerVars1= 500; //try dense gcd if size/#variables is bigger
4240
4241/// Extended Zassenhaus GCD for finite fields
4242CanonicalForm EZGCD_P( const CanonicalForm & FF, const CanonicalForm & GG )
4243{
4244  if (FF.isZero() && degree(GG) > 0) return GG/Lc(GG);
4245  else if (GG.isZero() && degree (FF) > 0) return FF/Lc(FF);
4246  else if (FF.isZero() && GG.isZero()) return FF.genOne();
4247  if (FF.inBaseDomain() || GG.inBaseDomain()) return FF.genOne();
4248  if (FF.isUnivariate() && fdivides(FF, GG)) return FF/Lc(FF);
4249  if (GG.isUnivariate() && fdivides(GG, FF)) return GG/Lc(GG);
4250  if (FF == GG) return FF/Lc(FF);
4251
4252  int maxNumVars= tmax (getNumVars (FF), getNumVars (GG));
4253  Variable a, oldA;
4254  int sizeF= size (FF);
4255  int sizeG= size (GG);
4256
4257  if (sizeF/maxNumVars > sizePerVars1 && sizeG/maxNumVars > sizePerVars1)
4258  {
4259    if (hasFirstAlgVar (FF, a) || hasFirstAlgVar (GG, a))
4260      return GCD_Fp_extension (FF, GG, a);
4261    else if (CFFactory::gettype() == GaloisFieldDomain)
4262      return GCD_GF (FF, GG);
4263    else
4264      return GCD_small_p (FF, GG);
4265  }
4266
4267  CanonicalForm F, G, f, g, d, Fb, Gb, Db, Fbt, Gbt, Dbt, B0, B, D0, lcF, lcG,
4268                lcD;
4269  CFArray DD( 1, 2 ), lcDD( 1, 2 );
4270  int degF, degG, delta, count;
4271  int maxeval;
4272  maxeval= tmin((getCharacteristic()/
4273                (int)(ilog2(getCharacteristic())*log2exp))*2, maxNumEval);
4274  count= 0; // number of eval. used
4275  REvaluation b, bt;
4276  int gcdfound = 0;
4277  Variable x = Variable(1);
4278
4279  F= FF;
4280  G= GG;
4281
4282  CFMap M,N;
4283  int smallestDegLev;
4284  TIMING_START (ez_p_compress)
4285  int best_level= compress4EZGCD (F, G, M, N, smallestDegLev);
4286
4287  if (best_level == 0) return G.genOne();
4288
4289  F= M (F);
4290  G= M (G);
4291  TIMING_END_AND_PRINT (ez_p_compress, "time for compression in EZ_P: ")
4292
4293  TIMING_START (ez_p_content)
4294  f = content( F, x ); g = content( G, x );
4295  d = gcd( f, g );
4296  F /= f; G /= g;
4297  TIMING_END_AND_PRINT (ez_p_content, "time to extract content in EZ_P: ")
4298
4299  if( F.isUnivariate() && G.isUnivariate() )
4300  {
4301    if( F.mvar() == G.mvar() )
4302      d *= gcd( F, G );
4303    else
4304      return N (d);
4305    return N (d);
4306  }
4307  if ( F.isUnivariate())
4308  {
4309    g= content (G,G.mvar());
4310    return N(d*gcd(F,g));
4311  }
4312  if ( G.isUnivariate())
4313  {
4314    f= content (F,F.mvar());
4315    return N(d*gcd(G,f));
4316  }
4317
4318  maxNumVars= tmax (getNumVars (F), getNumVars (G));
4319  sizeF= size (F);
4320  sizeG= size (G);
4321
4322  if (sizeF/maxNumVars > sizePerVars1 && sizeG/maxNumVars > sizePerVars1)
4323  {
4324    if (hasFirstAlgVar (F, a) || hasFirstAlgVar (G, a))
4325      return N (d*GCD_Fp_extension (F, G, a));
4326    else if (CFFactory::gettype() == GaloisFieldDomain)
4327      return N (d*GCD_GF (F, G));
4328    else
4329      return N (d*GCD_small_p (F, G));
4330  }
4331
4332  int dummy= 0;
4333  if( gcd_test_one( F, G, false, dummy ) )
4334  {
4335    return N (d);
4336  }
4337
4338  bool passToGF= false;
4339  bool extOfExt= false;
4340  int p= getCharacteristic();
4341  bool algExtension= (hasFirstAlgVar(F,a) || hasFirstAlgVar(G,a));
4342  int k= 1;
4343  CanonicalForm primElem, imPrimElem;
4344  CFList source, dest;
4345  if (p < 50 && CFFactory::gettype() != GaloisFieldDomain && !algExtension)
4346  {
4347    if (p == 2)
4348      setCharacteristic (2, 12, 'Z');
4349    else if (p == 3)
4350      setCharacteristic (3, 4, 'Z');
4351    else if (p == 5 || p == 7)
4352      setCharacteristic (p, 3, 'Z');
4353    else
4354      setCharacteristic (p, 2, 'Z');
4355    passToGF= true;
4356    F= F.mapinto();
4357    G= G.mapinto();
4358    maxeval= 2*ipower (p, getGFDegree());
4359  }
4360  else if (CFFactory::gettype() == GaloisFieldDomain &&
4361           ipower (p , getGFDegree()) < 50)
4362  {
4363    k= getGFDegree();
4364    if (ipower (p, 2*k) > 50)
4365      setCharacteristic (p, 2*k, gf_name);
4366    else
4367      setCharacteristic (p, 3*k, gf_name);
4368    F= GFMapUp (F, k);
4369    G= GFMapUp (G, k);
4370    maxeval= tmin (2*ipower (p, getGFDegree()), maxNumEval);
4371  }
4372  else if (p < 50 && algExtension && CFFactory::gettype() != GaloisFieldDomain)
4373  {
4374    int d= degree (getMipo (a));
4375    oldA= a;
4376    Variable v2;
4377    if (p == 2 && d < 6)
4378    {
4379      if (fac_NTL_char != p)
4380      {
4381        fac_NTL_char= p;
4382        zz_p::init (p);
4383      }
4384      bool primFail= false;
4385      Variable vBuf;
4386      primElem= primitiveElement (a, vBuf, primFail);
4387      ASSERT (!primFail, "failure in integer factorizer");
4388      if (d < 3)
4389      {
4390        zz_pX NTLIrredpoly;
4391        BuildIrred (NTLIrredpoly, d*3);
4392        CanonicalForm newMipo= convertNTLzzpX2CF (NTLIrredpoly, Variable (1));
4393        v2= rootOf (newMipo);
4394      }
4395      else
4396      {
4397        zz_pX NTLIrredpoly;
4398        BuildIrred (NTLIrredpoly, d*2);
4399        CanonicalForm newMipo= convertNTLzzpX2CF (NTLIrredpoly, Variable (1));
4400        v2= rootOf (newMipo);
4401      }
4402      imPrimElem= mapPrimElem (primElem, a, v2);
4403      extOfExt= true;
4404    }
4405    else if ((p == 3 && d < 4) || ((p == 5 || p == 7) && d < 3))
4406    {
4407      if (fac_NTL_char != p)
4408      {
4409        fac_NTL_char= p;
4410        zz_p::init (p);
4411      }
4412      bool primFail= false;
4413      Variable vBuf;
4414      primElem= primitiveElement (a, vBuf, primFail);
4415      ASSERT (!primFail, "failure in integer factorizer");
4416      zz_pX NTLIrredpoly;
4417      BuildIrred (NTLIrredpoly, d*2);
4418      CanonicalForm newMipo= convertNTLzzpX2CF (NTLIrredpoly, Variable (1));
4419      v2= rootOf (newMipo);
4420      imPrimElem= mapPrimElem (primElem, a, v2);
4421      extOfExt= true;
4422    }
4423    if (extOfExt)
4424    {
4425      maxeval= tmin (2*ipower (p, degree (getMipo (v2))), maxNumEval);
4426      F= mapUp (F, a, v2, primElem, imPrimElem, source, dest);
4427      G= mapUp (G, a, v2, primElem, imPrimElem, source, dest);
4428      a= v2;
4429    }
4430  }
4431
4432  lcF = LC( F, x ); lcG = LC( G, x );
4433  lcD = gcd( lcF, lcG );
4434
4435  delta = 0;
4436  degF = degree( F, x ); degG = degree( G, x );
4437
4438  if(hasFirstAlgVar(G,a))
4439    b = REvaluation( 2, tmax(F.level(), G.level()), AlgExtRandomF( a ) );
4440  else
4441  { // both not in extension given by algebraic variable
4442    if (CFFactory::gettype() != GaloisFieldDomain)
4443      b = REvaluation( 2, tmax(F.level(), G.level()), FFRandom() );
4444    else
4445      b = REvaluation( 2, tmax(F.level(), G.level()), GFRandom() );
4446  }
4447
4448  CanonicalForm cand, contcand;
4449  CanonicalForm result;
4450  int o, t;
4451  o= 0;
4452  t= 1;
4453  int goodPointCount= 0;
4454  while( !gcdfound )
4455  {
4456    TIMING_START (ez_p_eval);
4457    if( !findeval_P( F, G, Fb, Gb, Db, b, delta, degF, degG, maxeval, count, o,
4458         maxeval/maxNumVars, t ))
4459    { // too many eval. used --> try another method
4460      Off (SW_USE_EZGCD_P);
4461      result= gcd (F,G);
4462      On (SW_USE_EZGCD_P);
4463      if (passToGF)
4464      {
4465        CanonicalForm mipo= gf_mipo;
4466        setCharacteristic (p);
4467        Variable alpha= rootOf (mipo.mapinto());
4468        result= GF2FalphaRep (result, alpha);
4469      }
4470      if (k > 1)
4471      {
4472        result= GFMapDown (result, k);
4473        setCharacteristic (p, k, gf_name);
4474      }
4475      if (extOfExt)
4476        result= mapDown (result, primElem, imPrimElem, oldA, dest, source);
4477      return N (d*result);
4478    }
4479    TIMING_END_AND_PRINT (ez_p_eval, "time for eval point search in EZ_P1: ");
4480    delta = degree( Db );
4481    if( delta == 0 )
4482    {
4483      if (passToGF)
4484        setCharacteristic (p);
4485      if (k > 1)
4486        setCharacteristic (p, k, gf_name);
4487      return N (d);
4488    }
4489    while( true )
4490    {
4491      bt = b;
4492      TIMING_START (ez_p_eval);
4493      if( !findeval_P(F,G,Fbt,Gbt,Dbt, bt, delta, degF, degG, maxeval, count, o,
4494           maxeval/maxNumVars, t ))
4495      { // too many eval. used --> try another method
4496        Off (SW_USE_EZGCD_P);
4497        result= gcd (F,G);
4498        On (SW_USE_EZGCD_P);
4499        if (passToGF)
4500        {
4501          CanonicalForm mipo= gf_mipo;
4502          setCharacteristic (p);
4503          Variable alpha= rootOf (mipo.mapinto());
4504          result= GF2FalphaRep (result, alpha);
4505        }
4506        if (k > 1)
4507        {
4508          result= GFMapDown (result, k);
4509          setCharacteristic (p, k, gf_name);
4510        }
4511        if (extOfExt)
4512          result= mapDown (result, primElem, imPrimElem, oldA, dest, source);
4513        return N (d*result);
4514      }
4515      TIMING_END_AND_PRINT (ez_p_eval, "time for eval point search in EZ_P2: ");
4516      int dd = degree( Dbt );
4517      if( dd == 0 )
4518      {
4519        if (passToGF)
4520          setCharacteristic (p);
4521        if (k > 1)
4522          setCharacteristic (p, k, gf_name);
4523        return N (d);
4524      }
4525      if( dd == delta )
4526      {
4527        goodPointCount++;
4528        if (goodPointCount == 5)
4529          break;
4530      }
4531      if( dd < delta )
4532      {
4533        goodPointCount= 0;
4534        delta = dd;
4535        b = bt;
4536        Db = Dbt; Fb = Fbt; Gb = Gbt;
4537      }
4538      if (delta == degF)
4539      {
4540        if (degF <= degG && fdivides (F, G))
4541        {
4542          if (passToGF)
4543          {
4544            CanonicalForm mipo= gf_mipo;
4545            setCharacteristic (p);
4546            Variable alpha= rootOf (mipo.mapinto());
4547            F= GF2FalphaRep (F, alpha);
4548          }
4549          if (k > 1)
4550          {
4551            F= GFMapDown (F, k);
4552            setCharacteristic (p, k, gf_name);
4553          }
4554          if (extOfExt)
4555            F= mapDown (F, primElem, imPrimElem, oldA, dest, source);
4556          return N (d*F);
4557        }
4558        else
4559          delta--;
4560      }
4561      else if (delta == degG)
4562      {
4563        if (degG <= degF && fdivides (G, F))
4564        {
4565          if (passToGF)
4566          {
4567            CanonicalForm mipo= gf_mipo;
4568            setCharacteristic (p);
4569            Variable alpha= rootOf (mipo.mapinto());
4570            G= GF2FalphaRep (G, alpha);
4571          }
4572          if (k > 1)
4573          {
4574            G= GFMapDown (G, k);
4575            setCharacteristic (p, k, gf_name);
4576          }
4577          if (extOfExt)
4578            G= mapDown (G, primElem, imPrimElem, oldA, dest, source);
4579          return N (d*G);
4580        }
4581        else
4582          delta--;
4583      }
4584    }
4585    if( delta != degF && delta != degG )
4586    {
4587      bool B_is_F;
4588      CanonicalForm xxx1, xxx2;
4589      CanonicalForm buf;
4590      DD[1] = Fb / Db;
4591      buf= Gb/Db;
4592      xxx1 = gcd( DD[1], Db );
4593      xxx2 = gcd( buf, Db );
4594      if (((xxx1.inCoeffDomain() && xxx2.inCoeffDomain()) &&
4595          (size (F) <= size (G)))
4596          || (xxx1.inCoeffDomain() && !xxx2.inCoeffDomain()))
4597      {
4598        B = F;
4599        DD[2] = Db;
4600        lcDD[1] = lcF;
4601        lcDD[2] = lcD;
4602        B_is_F = true;
4603      }
4604      else if (((xxx1.inCoeffDomain() && xxx2.inCoeffDomain()) &&
4605               (size (G) < size (F)))
4606               || (!xxx1.inCoeffDomain() && xxx2.inCoeffDomain()))
4607      {
4608        DD[1] = buf;
4609        B = G;
4610        DD[2] = Db;
4611        lcDD[1] = lcG;
4612        lcDD[2] = lcD;
4613        B_is_F = false;
4614      }
4615      else // special case handling
4616      {
4617        Off (SW_USE_EZGCD_P);
4618        result= gcd (F,G);
4619        On (SW_USE_EZGCD_P);
4620        if (passToGF)
4621        {
4622          CanonicalForm mipo= gf_mipo;
4623          setCharacteristic (p);
4624          Variable alpha= rootOf (mipo.mapinto());
4625          result= GF2FalphaRep (result, alpha);
4626        }
4627        if (k > 1)
4628        {
4629          result= GFMapDown (result, k);
4630          setCharacteristic (p, k, gf_name);
4631        }
4632        if (extOfExt)
4633          result= mapDown (result, primElem, imPrimElem, oldA, dest, source);
4634        return N (d*result);
4635      }
4636      DD[2] = DD[2] * ( b( lcDD[2] ) / lc( DD[2] ) );
4637      DD[1] = DD[1] * ( b( lcDD[1] ) / lc( DD[1] ) );
4638
4639      if (size (B*lcDD[2])/maxNumVars > sizePerVars1)
4640      {
4641        if (algExtension)
4642        {
4643          result= GCD_Fp_extension (F, G, a);
4644          if (extOfExt)
4645            result= mapDown (result, primElem, imPrimElem, oldA, dest, source);
4646          return N (d*result);
4647        }
4648        if (CFFactory::gettype() == GaloisFieldDomain)
4649        {
4650          result= GCD_GF (F, G);
4651          if (passToGF)
4652          {
4653            CanonicalForm mipo= gf_mipo;
4654            setCharacteristic (p);
4655            Variable alpha= rootOf (mipo.mapinto());
4656            result= GF2FalphaRep (result, alpha);
4657          }
4658          if (k > 1)
4659          {
4660            result= GFMapDown (result, k);
4661            setCharacteristic (p, k, gf_name);
4662          }
4663          return N (d*result);
4664        }
4665        else
4666          return N (d*GCD_small_p (F,G));
4667      }
4668
4669      TIMING_START (ez_p_hensel_lift);
4670      gcdfound= Hensel_P (B*lcD, DD, b, lcDD);
4671      TIMING_END_AND_PRINT (ez_p_hensel_lift, "time for Hensel lift in EZ_P: ");
4672
4673      if (gcdfound == -1) //things became dense
4674      {
4675        if (algExtension)
4676        {
4677          result= GCD_Fp_extension (F, G, a);
4678          if (extOfExt)
4679            result= mapDown (result, primElem, imPrimElem, oldA, dest, source);
4680          return N (d*result);
4681        }
4682        if (CFFactory::gettype() == GaloisFieldDomain)
4683        {
4684          result= GCD_GF (F, G);
4685          if (passToGF)
4686          {
4687            CanonicalForm mipo= gf_mipo;
4688            setCharacteristic (p);
4689            Variable alpha= rootOf (mipo.mapinto());
4690            result= GF2FalphaRep (result, alpha);
4691          }
4692          if (k > 1)
4693          {
4694            result= GFMapDown (result, k);
4695            setCharacteristic (p, k, gf_name);
4696          }
4697          return N (d*result);
4698        }
4699        else
4700        {
4701          if (p >= cf_getBigPrime(0))
4702            return N (d*sparseGCDFp (F,G));
4703          else
4704            return N (d*GCD_small_p (F,G));
4705        }
4706      }
4707
4708      if (gcdfound == 1)
4709      {
4710        TIMING_START (termination_test);
4711        contcand= content (DD[2], Variable (1));
4712        cand = DD[2] / contcand;
4713        if (B_is_F)
4714          gcdfound = fdivides( cand, G ) && cand*(DD[1]/(lcD/contcand)) == F;
4715        else
4716          gcdfound = fdivides( cand, F ) && cand*(DD[1]/(lcD/contcand)) == G;
4717        TIMING_END_AND_PRINT (termination_test,
4718                              "time for termination test EZ_P: ");
4719
4720        if (passToGF && gcdfound)
4721        {
4722          CanonicalForm mipo= gf_mipo;
4723          setCharacteristic (p);
4724          Variable alpha= rootOf (mipo.mapinto());
4725          cand= GF2FalphaRep (cand, alpha);
4726        }
4727        if (k > 1 && gcdfound)
4728        {
4729          cand= GFMapDown (cand, k);
4730          setCharacteristic (p, k, gf_name);
4731        }
4732        if (extOfExt && gcdfound)
4733          cand= mapDown (cand, primElem, imPrimElem, oldA, dest, source);
4734      }
4735    }
4736    delta--;
4737    goodPointCount= 0;
4738  }
4739  return N (d*cand);
4740}
4741
4742
4743#endif
Note: See TracBrowser for help on using the repository browser.