source: git/factory/cf_gcd_smallp.cc @ 26b3b7

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