source: git/factory/cf_gcd_smallp.cc @ c1b9927

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