source: git/factory/cf_gcd_smallp.cc @ e368746

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