source: git/factory/cf_gcd_smallp.cc @ 473102

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