source: git/factory/cf_gcd_smallp.cc @ d08ed8

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