Loading...
Searching...
No Matches
cfModGcd.cc
Go to the documentation of this file.
// -*- c++ -*-
//*****************************************************************************
/** @file cfModGcd.cc
*
* This file implements the GCD of two polynomials over \f$ F_{p} \f$ ,
* \f$ F_{p}(\alpha ) \f$, GF or Z based on Alg. 7.1. and 7.2. as described in
* "Algorithms for Computer Algebra" by Geddes, Czapor, Labahn via modular
* computations. And sparse modular variants as described in Zippel
* "Effective Polynomial Computation", deKleine, Monagan, Wittkopf
* "Algorithms for the non-monic case of the sparse modular GCD algorithm" and
* Javadi "A new solution to the normalization problem"
*
* @author Martin Lee
* @date 22.10.2009
*
* @par Copyright:
* (c) by The SINGULAR Team, see LICENSE file
*
**/
//*****************************************************************************
#include "config.h"
#include <math.h>
#include "cf_assert.h"
#include "debug.h"
#include "timing.h"
#include "canonicalform.h"
#include "cfGcdUtil.h"
#include "cf_map.h"
#include "cf_util.h"
#include "cf_irred.h"
#include "templates/ftmpl_functions.h"
#include "cf_random.h"
#include "cf_reval.h"
#include "facHensel.h"
#include "cf_iter.h"
#include "cfNewtonPolygon.h"
#include "cf_algorithm.h"
#include "cf_primes.h"
// inline helper functions:
#include "cf_map_ext.h"
#ifdef HAVE_NTL
#include "NTLconvert.h"
#endif
#ifdef HAVE_FLINT
#include "FLINTconvert.h"
#endif
#include "cfModGcd.h"
TIMING_DEFINE_PRINT(gcd_recursion)
TIMING_DEFINE_PRINT(newton_interpolation)
TIMING_DEFINE_PRINT(termination_test)
TIMING_DEFINE_PRINT(ez_p_compress)
TIMING_DEFINE_PRINT(ez_p_hensel_lift)
TIMING_DEFINE_PRINT(ez_p_eval)
TIMING_DEFINE_PRINT(ez_p_content)
bool
terminationTest (const CanonicalForm& F, const CanonicalForm& G,
const CanonicalForm& coF, const CanonicalForm& coG,
const CanonicalForm& cand)
{
CanonicalForm LCCand= abs (LC (cand));
if (LCCand*abs (LC (coF)) == abs (LC (F)))
{
if (LCCand*abs (LC (coG)) == abs (LC (G)))
{
if (abs (cand)*abs (coF) == abs (F))
{
if (abs (cand)*abs (coG) == abs (G))
return true;
}
return false;
}
return false;
}
return false;
}
#if defined(HAVE_NTL)|| defined(HAVE_FLINT)
/// compressing two polynomials F and G, M is used for compressing,
/// N to reverse the compression
int myCompress (const CanonicalForm& F, const CanonicalForm& G, CFMap & M,
CFMap & N, bool topLevel)
{
int n= tmax (F.level(), G.level());
int * degsf= NEW_ARRAY(int,n + 1);
int * degsg= NEW_ARRAY(int,n + 1);
for (int i = n; i >= 0; i--)
degsf[i]= degsg[i]= 0;
degsf= degrees (F, degsf);
degsg= degrees (G, degsg);
int both_non_zero= 0;
int f_zero= 0;
int g_zero= 0;
int both_zero= 0;
if (topLevel)
{
for (int i= 1; i <= n; i++)
{
if (degsf[i] != 0 && degsg[i] != 0)
{
both_non_zero++;
continue;
}
if (degsf[i] == 0 && degsg[i] != 0 && i <= G.level())
{
f_zero++;
continue;
}
if (degsg[i] == 0 && degsf[i] && i <= F.level())
{
g_zero++;
continue;
}
}
if (both_non_zero == 0)
{
DELETE_ARRAY(degsf);
DELETE_ARRAY(degsg);
return 0;
}
// map Variables which do not occur in both polynomials to higher levels
int k= 1;
int l= 1;
for (int i= 1; i <= n; i++)
{
if (degsf[i] != 0 && degsg[i] == 0 && i <= F.level())
{
if (k + both_non_zero != i)
{
M.newpair (Variable (i), Variable (k + both_non_zero));
N.newpair (Variable (k + both_non_zero), Variable (i));
}
k++;
}
if (degsf[i] == 0 && degsg[i] != 0 && i <= G.level())
{
if (l + g_zero + both_non_zero != i)
{
M.newpair (Variable (i), Variable (l + g_zero + both_non_zero));
N.newpair (Variable (l + g_zero + both_non_zero), Variable (i));
}
l++;
}
}
// sort Variables x_{i} in increasing order of
// min(deg_{x_{i}}(f),deg_{x_{i}}(g))
int m= tmax (F.level(), G.level());
int min_max_deg;
k= both_non_zero;
l= 0;
int i= 1;
while (k > 0)
{
if (degsf [i] != 0 && degsg [i] != 0)
min_max_deg= tmax (degsf[i], degsg[i]);
else
min_max_deg= 0;
while (min_max_deg == 0)
{
i++;
if (degsf [i] != 0 && degsg [i] != 0)
min_max_deg= tmax (degsf[i], degsg[i]);
else
min_max_deg= 0;
}
for (int j= i + 1; j <= m; j++)
{
if (degsf[j] != 0 && degsg [j] != 0 &&
tmax (degsf[j],degsg[j]) <= min_max_deg)
{
min_max_deg= tmax (degsf[j], degsg[j]);
l= j;
}
}
if (l != 0)
{
if (l != k)
{
M.newpair (Variable (l), Variable(k));
N.newpair (Variable (k), Variable(l));
degsf[l]= 0;
degsg[l]= 0;
l= 0;
}
else
{
degsf[l]= 0;
degsg[l]= 0;
l= 0;
}
}
else if (l == 0)
{
if (i != k)
{
M.newpair (Variable (i), Variable (k));
N.newpair (Variable (k), Variable (i));
degsf[i]= 0;
degsg[i]= 0;
}
else
{
degsf[i]= 0;
degsg[i]= 0;
}
i++;
}
k--;
}
}
else
{
//arrange Variables such that no gaps occur
for (int i= 1; i <= n; i++)
{
if (degsf[i] == 0 && degsg[i] == 0)
{
both_zero++;
continue;
}
else
{
if (both_zero != 0)
{
M.newpair (Variable (i), Variable (i - both_zero));
N.newpair (Variable (i - both_zero), Variable (i));
}
}
}
}
DELETE_ARRAY(degsf);
DELETE_ARRAY(degsg);
return 1;
}
static inline CanonicalForm
uni_content (const CanonicalForm & F);
CanonicalForm
uni_content (const CanonicalForm& F, const Variable& x)
{
if (F.inCoeffDomain())
return F.genOne();
if (F.level() == x.level() && F.isUnivariate())
return F;
if (F.level() != x.level() && F.isUnivariate())
return F.genOne();
if (x.level() != 1)
{
CanonicalForm f= swapvar (F, x, Variable (1));
CanonicalForm result= uni_content (f);
return swapvar (result, x, Variable (1));
}
else
return uni_content (F);
}
/// compute the content of F, where F is considered as an element of
/// \f$ R[x_{1}][x_{2},\ldots ,x_{n}] \f$
static inline CanonicalForm
uni_content (const CanonicalForm & F)
{
if (F.inBaseDomain())
return F.genOne();
if (F.level() == 1 && F.isUnivariate())
return F;
if (F.level() != 1 && F.isUnivariate())
return F.genOne();
if (degree (F,1) == 0) return F.genOne();
int l= F.level();
if (l == 2)
return content(F);
else
{
CanonicalForm pol, c = 0;
CFIterator i = F;
for (; i.hasTerms(); i++)
{
pol= i.coeff();
pol= uni_content (pol);
c= gcd (c, pol);
if (c.isOne())
return c;
}
return c;
}
}
CanonicalForm
extractContents (const CanonicalForm& F, const CanonicalForm& G,
CanonicalForm& contentF, CanonicalForm& contentG,
CanonicalForm& ppF, CanonicalForm& ppG, const int d)
{
CanonicalForm uniContentF, uniContentG, gcdcFcG;
contentF= 1;
contentG= 1;
ppF= F;
ppG= G;
CanonicalForm result= 1;
for (int i= 1; i <= d; i++)
{
uniContentF= uni_content (F, Variable (i));
uniContentG= uni_content (G, Variable (i));
gcdcFcG= gcd (uniContentF, uniContentG);
contentF *= uniContentF;
contentG *= uniContentG;
ppF /= uniContentF;
ppG /= uniContentG;
result *= gcdcFcG;
}
return result;
}
/// compute the leading coefficient of F, where F is considered as an element
/// of \f$ R[x_{1}][x_{2},\ldots ,x_{n}] \f$, order on
/// \f$ Mon (x_{2},\ldots ,x_{n}) \f$ is dp.
static inline
CanonicalForm uni_lcoeff (const CanonicalForm& F)
{
if (F.level() > 1)
{
Variable x= Variable (2);
int deg= totaldegree (F, x, F.mvar());
for (CFIterator i= F; i.hasTerms(); i++)
{
if (i.exp() + totaldegree (i.coeff(), x, i.coeff().mvar()) == deg)
return uni_lcoeff (i.coeff());
}
}
return F;
}
/// Newton interpolation - Incremental algorithm.
/// Given a list of values alpha_i and a list of polynomials u_i, 1 <= i <= n,
/// computes the interpolation polynomial assuming that
/// the polynomials in u are the results of evaluating the variabe x
/// of the unknown polynomial at the alpha values.
/// This incremental version receives only the values of alpha_n and u_n and
/// the previous interpolation polynomial for points 1 <= i <= n-1, and computes
/// the polynomial interpolating in all the points.
/// newtonPoly must be equal to (x - alpha_1) * ... * (x - alpha_{n-1})
static inline CanonicalForm
newtonInterp(const CanonicalForm & alpha, const CanonicalForm & u,
const CanonicalForm & newtonPoly, const CanonicalForm & oldInterPoly,
const Variable & x)
{
CanonicalForm interPoly;
interPoly= oldInterPoly + ((u - oldInterPoly(alpha, x))/newtonPoly(alpha, x))
*newtonPoly;
return interPoly;
}
/// compute a random element a of \f$ F_{p}(\alpha ) \f$ ,
/// s.t. F(a) \f$ \neq 0 \f$ , F is a univariate polynomial, returns
/// fail if there are no field elements left which have not been used before
static inline CanonicalForm
randomElement (const CanonicalForm & F, const Variable & alpha, CFList & list,
bool & fail)
{
fail= false;
Variable x= F.mvar();
AlgExtRandomF genAlgExt (alpha);
FFRandom genFF;
CanonicalForm random, mipo;
mipo= getMipo (alpha);
int p= getCharacteristic ();
int d= degree (mipo);
double bound= pow ((double) p, (double) d);
do
{
if (list.length() == bound)
{
fail= true;
break;
}
if (list.length() < p)
{
random= genFF.generate();
while (find (list, random))
random= genFF.generate();
}
else
{
random= genAlgExt.generate();
while (find (list, random))
random= genAlgExt.generate();
}
if (F (random, x) == 0)
{
list.append (random);
continue;
}
} while (find (list, random));
return random;
}
static inline
Variable chooseExtension (const Variable & alpha)
{
int i, m;
// extension of F_p needed
if (alpha.level() == 1)
{
i= 1;
m= 2;
} //extension of F_p(alpha)
if (alpha.level() != 1)
{
i= 4;
m= degree (getMipo (alpha));
}
#ifdef HAVE_FLINT
nmod_poly_t Irredpoly;
nmod_poly_init(Irredpoly,getCharacteristic());
nmod_poly_randtest_monic_irreducible(Irredpoly, FLINTrandom, i*m+1);
CanonicalForm newMipo=convertnmod_poly_t2FacCF(Irredpoly,Variable(1));
nmod_poly_clear(Irredpoly);
#else
if (fac_NTL_char != getCharacteristic())
{
fac_NTL_char= getCharacteristic();
zz_p::init (getCharacteristic());
}
zz_pX NTLIrredpoly;
BuildIrred (NTLIrredpoly, i*m);
CanonicalForm newMipo= convertNTLzzpX2CF (NTLIrredpoly, Variable (1));
#endif
return rootOf (newMipo);
}
#if defined(HAVE_NTL) || defined(HAVE_FLINT)
CanonicalForm
modGCDFq (const CanonicalForm& F, const CanonicalForm& G,
CanonicalForm& coF, CanonicalForm& coG,
Variable & alpha, CFList& l, bool& topLevel);
#endif
#if defined(HAVE_NTL) || defined(HAVE_FLINT)
CanonicalForm
modGCDFq (const CanonicalForm& F, const CanonicalForm& G,
Variable & alpha, CFList& l, bool& topLevel)
{
CanonicalForm dummy1, dummy2;
CanonicalForm result= modGCDFq (F, G, dummy1, dummy2, alpha, l,
topLevel);
return result;
}
#endif
/// GCD of F and G over \f$ F_{p}(\alpha ) \f$ ,
/// l and topLevel are only used internally, output is monic
/// based on Alg. 7.2. as described in "Algorithms for
/// Computer Algebra" by Geddes, Czapor, Labahn
#if defined(HAVE_NTL) || defined(HAVE_FLINT)
CanonicalForm
modGCDFq (const CanonicalForm& F, const CanonicalForm& G,
CanonicalForm& coF, CanonicalForm& coG,
Variable & alpha, CFList& l, bool& topLevel)
{
CanonicalForm A= F;
CanonicalForm B= G;
if (F.isZero() && degree(G) > 0)
{
coF= 0;
coG= Lc (G);
return G/Lc(G);
}
else if (G.isZero() && degree (F) > 0)
{
coF= Lc (F);
coG= 0;
return F/Lc(F);
}
else if (F.isZero() && G.isZero())
{
coF= coG= 0;
return F.genOne();
}
if (F.inBaseDomain() || G.inBaseDomain())
{
coF= F;
coG= G;
return F.genOne();
}
if (F.isUnivariate() && fdivides(F, G, coG))
{
coF= Lc (F);
return F/Lc(F);
}
if (G.isUnivariate() && fdivides(G, F, coF))
{
coG= Lc (G);
return G/Lc(G);
}
if (F == G)
{
coF= coG= Lc (F);
return F/Lc(F);
}
CFMap M,N;
int best_level= myCompress (A, B, M, N, topLevel);
if (best_level == 0)
{
coF= F;
coG= G;
return B.genOne();
}
A= M(A);
B= M(B);
Variable x= Variable(1);
//univariate case
if (A.isUnivariate() && B.isUnivariate())
{
CanonicalForm result= gcd (A, B);
coF= N (A/result);
coG= N (B/result);
return N (result);
}
CanonicalForm cA, cB; // content of A and B
CanonicalForm ppA, ppB; // primitive part of A and B
CanonicalForm gcdcAcB;
cA = uni_content (A);
cB = uni_content (B);
gcdcAcB= gcd (cA, cB);
ppA= A/cA;
ppB= B/cB;
CanonicalForm lcA, lcB; // leading coefficients of A and B
CanonicalForm gcdlcAlcB;
lcA= uni_lcoeff (ppA);
lcB= uni_lcoeff (ppB);
gcdlcAlcB= gcd (lcA, lcB);
int d= totaldegree (ppA, Variable(2), Variable (ppA.level()));
if (d == 0)
{
coF= N (ppA*(cA/gcdcAcB));
coG= N (ppB*(cB/gcdcAcB));
return N(gcdcAcB);
}
int d0= totaldegree (ppB, Variable(2), Variable (ppB.level()));
if (d0 < d)
d= d0;
if (d == 0)
{
coF= N (ppA*(cA/gcdcAcB));
coG= N (ppB*(cB/gcdcAcB));
return N(gcdcAcB);
}
CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH;
CanonicalForm newtonPoly, coF_random_element, coG_random_element, coF_m,
coG_m, ppCoF, ppCoG;
newtonPoly= 1;
m= gcdlcAlcB;
G_m= 0;
coF= 0;
coG= 0;
H= 0;
bool fail= false;
topLevel= false;
bool inextension= false;
Variable V_buf= alpha, V_buf4= alpha;
CanonicalForm prim_elem, im_prim_elem;
CanonicalForm prim_elem_alpha, im_prim_elem_alpha;
CFList source, dest;
int bound1= degree (ppA, 1);
int bound2= degree (ppB, 1);
do
{
random_element= randomElement (m*lcA*lcB, V_buf, l, fail);
if (fail)
{
source= CFList();
dest= CFList();
Variable V_buf3= V_buf;
V_buf= chooseExtension (V_buf);
bool prim_fail= false;
Variable V_buf2;
prim_elem= primitiveElement (V_buf4, V_buf2, prim_fail);
if (V_buf4 == alpha)
prim_elem_alpha= prim_elem;
if (V_buf3 != V_buf4)
{
m= mapDown (m, prim_elem, im_prim_elem, V_buf4, source, dest);
G_m= mapDown (G_m, prim_elem, im_prim_elem, V_buf4, source, dest);
coF_m= mapDown (coF_m, prim_elem, im_prim_elem, V_buf4, source, dest);
coG_m= mapDown (coG_m, prim_elem, im_prim_elem, V_buf4, source, dest);
newtonPoly= mapDown (newtonPoly, prim_elem, im_prim_elem, V_buf4,
source, dest);
ppA= mapDown (ppA, prim_elem, im_prim_elem, V_buf4, source, dest);
ppB= mapDown (ppB, prim_elem, im_prim_elem, V_buf4, source, dest);
gcdlcAlcB= mapDown (gcdlcAlcB, prim_elem, im_prim_elem, V_buf4,
source, dest);
lcA= mapDown (lcA, prim_elem, im_prim_elem, V_buf4, source, dest);
lcB= mapDown (lcB, prim_elem, im_prim_elem, V_buf4, source, dest);
for (CFListIterator i= l; i.hasItem(); i++)
i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, V_buf4,
source, dest);
}
ASSERT (!prim_fail, "failure in integer factorizer");
if (prim_fail)
; //ERROR
else
im_prim_elem= mapPrimElem (prim_elem, V_buf4, V_buf);
if (V_buf4 == alpha)
im_prim_elem_alpha= im_prim_elem;
else
im_prim_elem_alpha= mapUp (im_prim_elem_alpha, V_buf4, V_buf, prim_elem,
im_prim_elem, source, dest);
DEBOUTLN (cerr, "getMipo (V_buf4)= " << getMipo (V_buf4));
DEBOUTLN (cerr, "getMipo (V_buf2)= " << getMipo (V_buf2));
inextension= true;
for (CFListIterator i= l; i.hasItem(); i++)
i.getItem()= mapUp (i.getItem(), V_buf4, V_buf, prim_elem,
im_prim_elem, source, dest);
m= mapUp (m, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
G_m= mapUp (G_m, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
coF_m= mapUp (coF_m, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
coG_m= mapUp (coG_m, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
newtonPoly= mapUp (newtonPoly, V_buf4, V_buf, prim_elem, im_prim_elem,
source, dest);
ppA= mapUp (ppA, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
ppB= mapUp (ppB, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
gcdlcAlcB= mapUp (gcdlcAlcB, V_buf4, V_buf, prim_elem, im_prim_elem,
source, dest);
lcA= mapUp (lcA, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
lcB= mapUp (lcB, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
fail= false;
random_element= randomElement (m*lcA*lcB, V_buf, l, fail );
DEBOUTLN (cerr, "fail= " << fail);
CFList list;
TIMING_START (gcd_recursion);
G_random_element=
modGCDFq (ppA (random_element, x), ppB (random_element, x),
coF_random_element, coG_random_element, V_buf,
list, topLevel);
TIMING_END_AND_PRINT (gcd_recursion,
"time for recursive call: ");
DEBOUTLN (cerr, "G_random_element= " << G_random_element);
V_buf4= V_buf;
}
else
{
CFList list;
TIMING_START (gcd_recursion);
G_random_element=
modGCDFq (ppA(random_element, x), ppB(random_element, x),
coF_random_element, coG_random_element, V_buf,
list, topLevel);
TIMING_END_AND_PRINT (gcd_recursion,
"time for recursive call: ");
DEBOUTLN (cerr, "G_random_element= " << G_random_element);
}
if (!G_random_element.inCoeffDomain())
d0= totaldegree (G_random_element, Variable(2),
Variable (G_random_element.level()));
else
d0= 0;
if (d0 == 0)
{
if (inextension)
{
CFList u, v;
ppA= mapDown (ppA, prim_elem_alpha, im_prim_elem_alpha, alpha, u, v);
ppB= mapDown (ppB, prim_elem_alpha, im_prim_elem_alpha, alpha, u, v);
prune1 (alpha);
}
coF= N (ppA*(cA/gcdcAcB));
coG= N (ppB*(cB/gcdcAcB));
return N(gcdcAcB);
}
if (d0 > d)
{
if (!find (l, random_element))
l.append (random_element);
continue;
}
G_random_element=
(gcdlcAlcB(random_element, x)/uni_lcoeff (G_random_element))
* G_random_element;
coF_random_element= (lcA(random_element,x)/uni_lcoeff(coF_random_element))
*coF_random_element;
coG_random_element= (lcB(random_element,x)/uni_lcoeff(coG_random_element))
*coG_random_element;
if (!G_random_element.inCoeffDomain())
d0= totaldegree (G_random_element, Variable(2),
Variable (G_random_element.level()));
else
d0= 0;
if (d0 < d)
{
m= gcdlcAlcB;
newtonPoly= 1;
G_m= 0;
d= d0;
coF_m= 0;
coG_m= 0;
}
TIMING_START (newton_interpolation);
H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
coF= newtonInterp (random_element, coF_random_element, newtonPoly, coF_m,x);
coG= newtonInterp (random_element, coG_random_element, newtonPoly, coG_m,x);
TIMING_END_AND_PRINT (newton_interpolation,
"time for newton interpolation: ");
//termination test
if ((uni_lcoeff (H) == gcdlcAlcB) || (G_m == H))
{
TIMING_START (termination_test);
if (gcdlcAlcB.isOne())
cH= 1;
else
cH= uni_content (H);
ppH= H/cH;
ppH /= Lc (ppH);
CanonicalForm lcppH= gcdlcAlcB/cH;
CanonicalForm ccoF= lcppH/Lc (lcppH);
CanonicalForm ccoG= lcppH/Lc (lcppH);
ppCoF= coF/ccoF;
ppCoG= coG/ccoG;
if (inextension)
{
if (((degree (ppCoF,1)+degree (ppH,1) == bound1) &&
(degree (ppCoG,1)+degree (ppH,1) == bound2) &&
terminationTest (ppA, ppB, ppCoF, ppCoG, ppH)) ||
(fdivides (ppH, ppA, ppCoF) && fdivides (ppH, ppB, ppCoG)))
{
CFList u, v;
DEBOUTLN (cerr, "ppH before mapDown= " << ppH);
ppH= mapDown (ppH, prim_elem_alpha, im_prim_elem_alpha, alpha, u, v);
ppCoF= mapDown (ppCoF, prim_elem_alpha, im_prim_elem_alpha, alpha, u, v);
ppCoG= mapDown (ppCoG, prim_elem_alpha, im_prim_elem_alpha, alpha, u, v);
DEBOUTLN (cerr, "ppH after mapDown= " << ppH);
coF= N ((cA/gcdcAcB)*ppCoF);
coG= N ((cB/gcdcAcB)*ppCoG);
TIMING_END_AND_PRINT (termination_test,
"time for successful termination test Fq: ");
prune1 (alpha);
return N(gcdcAcB*ppH);
}
}
else if (((degree (ppCoF,1)+degree (ppH,1) == bound1) &&
(degree (ppCoG,1)+degree (ppH,1) == bound2) &&
terminationTest (ppA, ppB, ppCoF, ppCoG, ppH)) ||
(fdivides (ppH, ppA, ppCoF) && fdivides (ppH, ppB, ppCoG)))
{
coF= N ((cA/gcdcAcB)*ppCoF);
coG= N ((cB/gcdcAcB)*ppCoG);
TIMING_END_AND_PRINT (termination_test,
"time for successful termination test Fq: ");
return N(gcdcAcB*ppH);
}
TIMING_END_AND_PRINT (termination_test,
"time for unsuccessful termination test Fq: ");
}
G_m= H;
coF_m= coF;
coG_m= coG;
newtonPoly= newtonPoly*(x - random_element);
m= m*(x - random_element);
if (!find (l, random_element))
l.append (random_element);
} while (1);
}
#endif
/// compute a random element a of GF, s.t. F(a) \f$ \neq 0 \f$ , F is a
/// univariate polynomial, returns fail if there are no field elements left
/// which have not been used before
static inline
CanonicalForm
GFRandomElement (const CanonicalForm& F, CFList& list, bool& fail)
{
fail= false;
Variable x= F.mvar();
GFRandom genGF;
CanonicalForm random;
int p= getCharacteristic();
int d= getGFDegree();
int bound= ipower (p, d);
do
{
if (list.length() == bound)
{
fail= true;
break;
}
if (list.length() < 1)
random= 0;
else
{
random= genGF.generate();
while (find (list, random))
random= genGF.generate();
}
if (F (random, x) == 0)
{
list.append (random);
continue;
}
} while (find (list, random));
return random;
}
CanonicalForm
modGCDGF (const CanonicalForm& F, const CanonicalForm& G,
CanonicalForm& coF, CanonicalForm& coG,
CFList& l, bool& topLevel);
CanonicalForm
modGCDGF (const CanonicalForm& F, const CanonicalForm& G, CFList& l,
bool& topLevel)
{
CanonicalForm dummy1, dummy2;
CanonicalForm result= modGCDGF (F, G, dummy1, dummy2, l, topLevel);
return result;
}
/// GCD of F and G over GF, based on Alg. 7.2. as described in "Algorithms for
/// Computer Algebra" by Geddes, Czapor, Labahn
/// Usually this algorithm will be faster than modGCDFq since GF has
/// faster field arithmetics, however it might fail if the input is large since
/// the size of the base field is bounded by 2^16, output is monic
CanonicalForm
modGCDGF (const CanonicalForm& F, const CanonicalForm& G,
CanonicalForm& coF, CanonicalForm& coG,
CFList& l, bool& topLevel)
{
CanonicalForm A= F;
CanonicalForm B= G;
if (F.isZero() && degree(G) > 0)
{
coF= 0;
coG= Lc (G);
return G/Lc(G);
}
else if (G.isZero() && degree (F) > 0)
{
coF= Lc (F);
coG= 0;
return F/Lc(F);
}
else if (F.isZero() && G.isZero())
{
coF= coG= 0;
return F.genOne();
}
if (F.inBaseDomain() || G.inBaseDomain())
{
coF= F;
coG= G;
return F.genOne();
}
if (F.isUnivariate() && fdivides(F, G, coG))
{
coF= Lc (F);
return F/Lc(F);
}
if (G.isUnivariate() && fdivides(G, F, coF))
{
coG= Lc (G);
return G/Lc(G);
}
if (F == G)
{
coF= coG= Lc (F);
return F/Lc(F);
}
CFMap M,N;
int best_level= myCompress (A, B, M, N, topLevel);
if (best_level == 0)
{
coF= F;
coG= G;
return B.genOne();
}
A= M(A);
B= M(B);
Variable x= Variable(1);
//univariate case
if (A.isUnivariate() && B.isUnivariate())
{
CanonicalForm result= gcd (A, B);
coF= N (A/result);
coG= N (B/result);
return N (result);
}
CanonicalForm cA, cB; // content of A and B
CanonicalForm ppA, ppB; // primitive part of A and B
CanonicalForm gcdcAcB;
cA = uni_content (A);
cB = uni_content (B);
gcdcAcB= gcd (cA, cB);
ppA= A/cA;
ppB= B/cB;
CanonicalForm lcA, lcB; // leading coefficients of A and B
CanonicalForm gcdlcAlcB;
lcA= uni_lcoeff (ppA);
lcB= uni_lcoeff (ppB);
gcdlcAlcB= gcd (lcA, lcB);
int d= totaldegree (ppA, Variable(2), Variable (ppA.level()));
if (d == 0)
{
coF= N (ppA*(cA/gcdcAcB));
coG= N (ppB*(cB/gcdcAcB));
return N(gcdcAcB);
}
int d0= totaldegree (ppB, Variable(2), Variable (ppB.level()));
if (d0 < d)
d= d0;
if (d == 0)
{
coF= N (ppA*(cA/gcdcAcB));
coG= N (ppB*(cB/gcdcAcB));
return N(gcdcAcB);
}
CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH;
CanonicalForm newtonPoly, coF_random_element, coG_random_element, coF_m,
coG_m, ppCoF, ppCoG;
newtonPoly= 1;
m= gcdlcAlcB;
G_m= 0;
coF= 0;
coG= 0;
H= 0;
bool fail= false;
topLevel= false;
bool inextension= false;
int p=-1;
int k= getGFDegree();
int kk;
int expon;
char gf_name_buf= gf_name;
int bound1= degree (ppA, 1);
int bound2= degree (ppB, 1);
do
{
random_element= GFRandomElement (m*lcA*lcB, l, fail);
if (fail)
{
p= getCharacteristic();
expon= 2;
kk= getGFDegree();
if (ipower (p, kk*expon) < (1 << 16))
setCharacteristic (p, kk*(int)expon, 'b');
else
{
expon= (int) floor((log ((double)((1<<16) - 1)))/(log((double)p)*kk));
ASSERT (expon >= 2, "not enough points in modGCDGF");
setCharacteristic (p, (int)(kk*expon), 'b');
}
inextension= true;
fail= false;
for (CFListIterator i= l; i.hasItem(); i++)
i.getItem()= GFMapUp (i.getItem(), kk);
m= GFMapUp (m, kk);
G_m= GFMapUp (G_m, kk);
newtonPoly= GFMapUp (newtonPoly, kk);
coF_m= GFMapUp (coF_m, kk);
coG_m= GFMapUp (coG_m, kk);
ppA= GFMapUp (ppA, kk);
ppB= GFMapUp (ppB, kk);
gcdlcAlcB= GFMapUp (gcdlcAlcB, kk);
lcA= GFMapUp (lcA, kk);
lcB= GFMapUp (lcB, kk);
random_element= GFRandomElement (m*lcA*lcB, l, fail);
DEBOUTLN (cerr, "fail= " << fail);
CFList list;
TIMING_START (gcd_recursion);
G_random_element= modGCDGF (ppA(random_element, x), ppB(random_element, x),
coF_random_element, coG_random_element,
list, topLevel);
TIMING_END_AND_PRINT (gcd_recursion,
"time for recursive call: ");
DEBOUTLN (cerr, "G_random_element= " << G_random_element);
}
else
{
CFList list;
TIMING_START (gcd_recursion);
G_random_element= modGCDGF (ppA(random_element, x), ppB(random_element, x),
coF_random_element, coG_random_element,
list, topLevel);
TIMING_END_AND_PRINT (gcd_recursion,
"time for recursive call: ");
DEBOUTLN (cerr, "G_random_element= " << G_random_element);
}
if (!G_random_element.inCoeffDomain())
d0= totaldegree (G_random_element, Variable(2),
Variable (G_random_element.level()));
else
d0= 0;
if (d0 == 0)
{
if (inextension)
{
ppA= GFMapDown (ppA, k);
ppB= GFMapDown (ppB, k);
setCharacteristic (p, k, gf_name_buf);
}
coF= N (ppA*(cA/gcdcAcB));
coG= N (ppB*(cB/gcdcAcB));
return N(gcdcAcB);
}
if (d0 > d)
{
if (!find (l, random_element))
l.append (random_element);
continue;
}
G_random_element=
(gcdlcAlcB(random_element, x)/uni_lcoeff(G_random_element)) *
G_random_element;
coF_random_element= (lcA(random_element,x)/uni_lcoeff(coF_random_element))
*coF_random_element;
coG_random_element= (lcB(random_element,x)/uni_lcoeff(coG_random_element))
*coG_random_element;
if (!G_random_element.inCoeffDomain())
d0= totaldegree (G_random_element, Variable(2),
Variable (G_random_element.level()));
else
d0= 0;
if (d0 < d)
{
m= gcdlcAlcB;
newtonPoly= 1;
G_m= 0;
d= d0;
coF_m= 0;
coG_m= 0;
}
TIMING_START (newton_interpolation);
H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
coF= newtonInterp (random_element, coF_random_element, newtonPoly, coF_m,x);
coG= newtonInterp (random_element, coG_random_element, newtonPoly, coG_m,x);
TIMING_END_AND_PRINT (newton_interpolation,
"time for newton interpolation: ");
//termination test
if ((uni_lcoeff (H) == gcdlcAlcB) || (G_m == H))
{
TIMING_START (termination_test);
if (gcdlcAlcB.isOne())
cH= 1;
else
cH= uni_content (H);
ppH= H/cH;
ppH /= Lc (ppH);
CanonicalForm lcppH= gcdlcAlcB/cH;
CanonicalForm ccoF= lcppH/Lc (lcppH);
CanonicalForm ccoG= lcppH/Lc (lcppH);
ppCoF= coF/ccoF;
ppCoG= coG/ccoG;
if (inextension)
{
if (((degree (ppCoF,1)+degree (ppH,1) == bound1) &&
(degree (ppCoG,1)+degree (ppH,1) == bound2) &&
terminationTest (ppA, ppB, ppCoF, ppCoG, ppH)) ||
(fdivides (ppH, ppA, ppCoF) && fdivides (ppH, ppB, ppCoG)))
{
DEBOUTLN (cerr, "ppH before mapDown= " << ppH);
ppH= GFMapDown (ppH, k);
ppCoF= GFMapDown (ppCoF, k);
ppCoG= GFMapDown (ppCoG, k);
DEBOUTLN (cerr, "ppH after mapDown= " << ppH);
coF= N ((cA/gcdcAcB)*ppCoF);
coG= N ((cB/gcdcAcB)*ppCoG);
setCharacteristic (p, k, gf_name_buf);
TIMING_END_AND_PRINT (termination_test,
"time for successful termination GF: ");
return N(gcdcAcB*ppH);
}
}
else
{
if (((degree (ppCoF,1)+degree (ppH,1) == bound1) &&
(degree (ppCoG,1)+degree (ppH,1) == bound2) &&
terminationTest (ppA, ppB, ppCoF, ppCoG, ppH)) ||
(fdivides (ppH, ppA, ppCoF) && fdivides (ppH, ppB, ppCoG)))
{
coF= N ((cA/gcdcAcB)*ppCoF);
coG= N ((cB/gcdcAcB)*ppCoG);
TIMING_END_AND_PRINT (termination_test,
"time for successful termination GF: ");
return N(gcdcAcB*ppH);
}
}
TIMING_END_AND_PRINT (termination_test,
"time for unsuccessful termination GF: ");
}
G_m= H;
coF_m= coF;
coG_m= coG;
newtonPoly= newtonPoly*(x - random_element);
m= m*(x - random_element);
if (!find (l, random_element))
l.append (random_element);
} while (1);
}
static inline
CanonicalForm
FpRandomElement (const CanonicalForm& F, CFList& list, bool& fail)
{
fail= false;
Variable x= F.mvar();
FFRandom genFF;
CanonicalForm random;
int p= getCharacteristic();
int bound= p;
do
{
if (list.length() == bound)
{
fail= true;
break;
}
if (list.length() < 1)
random= 0;
else
{
random= genFF.generate();
while (find (list, random))
random= genFF.generate();
}
if (F (random, x) == 0)
{
list.append (random);
continue;
}
} while (find (list, random));
return random;
}
#if defined(HAVE_NTL) || defined(HAVE_FLINT)
CanonicalForm
modGCDFp (const CanonicalForm& F, const CanonicalForm& G,
CanonicalForm& coF, CanonicalForm& coG,
bool& topLevel, CFList& l);
#endif
#if defined(HAVE_NTL) || defined(HAVE_FLINT)
CanonicalForm
modGCDFp (const CanonicalForm& F, const CanonicalForm& G,
bool& topLevel, CFList& l)
{
CanonicalForm dummy1, dummy2;
CanonicalForm result= modGCDFp (F, G, dummy1, dummy2, topLevel, l);
return result;
}
#endif
#if defined(HAVE_NTL) || defined(HAVE_FLINT)
CanonicalForm
modGCDFp (const CanonicalForm& F, const CanonicalForm& G,
CanonicalForm& coF, CanonicalForm& coG,
bool& topLevel, CFList& l)
{
CanonicalForm A= F;
CanonicalForm B= G;
if (F.isZero() && degree(G) > 0)
{
coF= 0;
coG= Lc (G);
return G/Lc(G);
}
else if (G.isZero() && degree (F) > 0)
{
coF= Lc (F);
coG= 0;
return F/Lc(F);
}
else if (F.isZero() && G.isZero())
{
coF= coG= 0;
return F.genOne();
}
if (F.inBaseDomain() || G.inBaseDomain())
{
coF= F;
coG= G;
return F.genOne();
}
if (F.isUnivariate() && fdivides(F, G, coG))
{
coF= Lc (F);
return F/Lc(F);
}
if (G.isUnivariate() && fdivides(G, F, coF))
{
coG= Lc (G);
return G/Lc(G);
}
if (F == G)
{
coF= coG= Lc (F);
return F/Lc(F);
}
CFMap M,N;
int best_level= myCompress (A, B, M, N, topLevel);
if (best_level == 0)
{
coF= F;
coG= G;
return B.genOne();
}
A= M(A);
B= M(B);
Variable x= Variable (1);
//univariate case
if (A.isUnivariate() && B.isUnivariate())
{
CanonicalForm result= gcd (A, B);
coF= N (A/result);
coG= N (B/result);
return N (result);
}
CanonicalForm cA, cB; // content of A and B
CanonicalForm ppA, ppB; // primitive part of A and B
CanonicalForm gcdcAcB;
cA = uni_content (A);
cB = uni_content (B);
gcdcAcB= gcd (cA, cB);
ppA= A/cA;
ppB= B/cB;
CanonicalForm lcA, lcB; // leading coefficients of A and B
CanonicalForm gcdlcAlcB;
lcA= uni_lcoeff (ppA);
lcB= uni_lcoeff (ppB);
gcdlcAlcB= gcd (lcA, lcB);
int d= totaldegree (ppA, Variable (2), Variable (ppA.level()));
int d0;
if (d == 0)
{
coF= N (ppA*(cA/gcdcAcB));
coG= N (ppB*(cB/gcdcAcB));
return N(gcdcAcB);
}
d0= totaldegree (ppB, Variable (2), Variable (ppB.level()));
if (d0 < d)
d= d0;
if (d == 0)
{
coF= N (ppA*(cA/gcdcAcB));
coG= N (ppB*(cB/gcdcAcB));
return N(gcdcAcB);
}
CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH;
CanonicalForm newtonPoly, coF_random_element, coG_random_element,
coF_m, coG_m, ppCoF, ppCoG;
newtonPoly= 1;
m= gcdlcAlcB;
H= 0;
coF= 0;
coG= 0;
G_m= 0;
Variable alpha, V_buf, cleanUp;
bool fail= false;
bool inextension= false;
topLevel= false;
CFList source, dest;
int bound1= degree (ppA, 1);
int bound2= degree (ppB, 1);
do
{
if (inextension)
random_element= randomElement (m*lcA*lcB, V_buf, l, fail);
else
random_element= FpRandomElement (m*lcA*lcB, l, fail);
if (!fail && !inextension)
{
CFList list;
TIMING_START (gcd_recursion);
G_random_element=
modGCDFp (ppA (random_element,x), ppB (random_element,x),
coF_random_element, coG_random_element, topLevel,
list);
TIMING_END_AND_PRINT (gcd_recursion,
"time for recursive call: ");
DEBOUTLN (cerr, "G_random_element= " << G_random_element);
}
else if (!fail && inextension)
{
CFList list;
TIMING_START (gcd_recursion);
G_random_element=
modGCDFq (ppA (random_element, x), ppB (random_element, x),
coF_random_element, coG_random_element, V_buf,
list, topLevel);
TIMING_END_AND_PRINT (gcd_recursion,
"time for recursive call: ");
DEBOUTLN (cerr, "G_random_element= " << G_random_element);
}
else if (fail && !inextension)
{
source= CFList();
dest= CFList();
CFList list;
CanonicalForm mipo;
int deg= 2;
bool initialized= false;
do
{
mipo= randomIrredpoly (deg, x);
if (initialized)
setMipo (alpha, mipo);
else
alpha= rootOf (mipo);
inextension= true;
initialized= true;
fail= false;
random_element= randomElement (m*lcA*lcB, alpha, l, fail);
deg++;
} while (fail);
list= CFList();
V_buf= alpha;
cleanUp= alpha;
TIMING_START (gcd_recursion);
G_random_element=
modGCDFq (ppA (random_element, x), ppB (random_element, x),
coF_random_element, coG_random_element, alpha,
list, topLevel);
TIMING_END_AND_PRINT (gcd_recursion,
"time for recursive call: ");
DEBOUTLN (cerr, "G_random_element= " << G_random_element);
}
else if (fail && inextension)
{
source= CFList();
dest= CFList();
Variable V_buf3= V_buf;
V_buf= chooseExtension (V_buf);
bool prim_fail= false;
Variable V_buf2;
CanonicalForm prim_elem, im_prim_elem;
prim_elem= primitiveElement (alpha, V_buf2, prim_fail);
if (V_buf3 != alpha)
{
m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
G_m= mapDown (G_m, prim_elem, im_prim_elem, alpha, source, dest);
coF_m= mapDown (coF_m, prim_elem, im_prim_elem, alpha, source, dest);
coG_m= mapDown (coG_m, prim_elem, im_prim_elem, alpha, source, dest);
newtonPoly= mapDown (newtonPoly, prim_elem, im_prim_elem, alpha,
source, dest);
ppA= mapDown (ppA, prim_elem, im_prim_elem, alpha, source, dest);
ppB= mapDown (ppB, prim_elem, im_prim_elem, alpha, source, dest);
gcdlcAlcB= mapDown (gcdlcAlcB, prim_elem, im_prim_elem, alpha, source,
dest);
lcA= mapDown (lcA, prim_elem, im_prim_elem, alpha, source, dest);
lcB= mapDown (lcB, prim_elem, im_prim_elem, alpha, source, dest);
for (CFListIterator i= l; i.hasItem(); i++)
i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, alpha,
source, dest);
}
ASSERT (!prim_fail, "failure in integer factorizer");
if (prim_fail)
; //ERROR
else
im_prim_elem= mapPrimElem (prim_elem, alpha, V_buf);
DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (alpha));
DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf2));
for (CFListIterator i= l; i.hasItem(); i++)
i.getItem()= mapUp (i.getItem(), alpha, V_buf, prim_elem,
im_prim_elem, source, dest);
m= mapUp (m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
G_m= mapUp (G_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
coF_m= mapUp (coF_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
coG_m= mapUp (coG_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
newtonPoly= mapUp (newtonPoly, alpha, V_buf, prim_elem, im_prim_elem,
source, dest);
ppA= mapUp (ppA, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
ppB= mapUp (ppB, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
gcdlcAlcB= mapUp (gcdlcAlcB, alpha, V_buf, prim_elem, im_prim_elem,
source, dest);
lcA= mapUp (lcA, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
lcB= mapUp (lcB, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
fail= false;
random_element= randomElement (m*lcA*lcB, V_buf, l, fail );
DEBOUTLN (cerr, "fail= " << fail);
CFList list;
TIMING_START (gcd_recursion);
G_random_element=
modGCDFq (ppA (random_element, x), ppB (random_element, x),
coF_random_element, coG_random_element, V_buf,
list, topLevel);
TIMING_END_AND_PRINT (gcd_recursion,
"time for recursive call: ");
DEBOUTLN (cerr, "G_random_element= " << G_random_element);
alpha= V_buf;
}
if (!G_random_element.inCoeffDomain())
d0= totaldegree (G_random_element, Variable(2),
Variable (G_random_element.level()));
else
d0= 0;
if (d0 == 0)
{
if (inextension)
prune (cleanUp);
coF= N (ppA*(cA/gcdcAcB));
coG= N (ppB*(cB/gcdcAcB));
return N(gcdcAcB);
}
if (d0 > d)
{
if (!find (l, random_element))
l.append (random_element);
continue;
}
G_random_element= (gcdlcAlcB(random_element,x)/uni_lcoeff(G_random_element))
*G_random_element;
coF_random_element= (lcA(random_element,x)/uni_lcoeff(coF_random_element))
*coF_random_element;
coG_random_element= (lcB(random_element,x)/uni_lcoeff(coG_random_element))
*coG_random_element;
if (!G_random_element.inCoeffDomain())
d0= totaldegree (G_random_element, Variable(2),
Variable (G_random_element.level()));
else
d0= 0;
if (d0 < d)
{
m= gcdlcAlcB;
newtonPoly= 1;
G_m= 0;
d= d0;
coF_m= 0;
coG_m= 0;
}
TIMING_START (newton_interpolation);
H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
coF= newtonInterp (random_element, coF_random_element, newtonPoly, coF_m,x);
coG= newtonInterp (random_element, coG_random_element, newtonPoly, coG_m,x);
TIMING_END_AND_PRINT (newton_interpolation,
"time for newton_interpolation: ");
//termination test
if ((uni_lcoeff (H) == gcdlcAlcB) || (G_m == H))
{
TIMING_START (termination_test);
if (gcdlcAlcB.isOne())
cH= 1;
else
cH= uni_content (H);
ppH= H/cH;
ppH /= Lc (ppH);
CanonicalForm lcppH= gcdlcAlcB/cH;
CanonicalForm ccoF= lcppH/Lc (lcppH);
CanonicalForm ccoG= lcppH/Lc (lcppH);
ppCoF= coF/ccoF;
ppCoG= coG/ccoG;
DEBOUTLN (cerr, "ppH= " << ppH);
if (((degree (ppCoF,1)+degree (ppH,1) == bound1) &&
(degree (ppCoG,1)+degree (ppH,1) == bound2) &&
terminationTest (ppA, ppB, ppCoF, ppCoG, ppH)) ||
(fdivides (ppH, ppA, ppCoF) && fdivides (ppH, ppB, ppCoG)))
{
if (inextension)
prune (cleanUp);
coF= N ((cA/gcdcAcB)*ppCoF);
coG= N ((cB/gcdcAcB)*ppCoG);
TIMING_END_AND_PRINT (termination_test,
"time for successful termination Fp: ");
return N(gcdcAcB*ppH);
}
TIMING_END_AND_PRINT (termination_test,
"time for unsuccessful termination Fp: ");
}
G_m= H;
coF_m= coF;
coG_m= coG;
newtonPoly= newtonPoly*(x - random_element);
m= m*(x - random_element);
if (!find (l, random_element))
l.append (random_element);
} while (1);
}
#endif
CFArray
solveVandermonde (const CFArray& M, const CFArray& A)
{
int r= M.size();
ASSERT (A.size() == r, "vector does not have right size");
if (r == 1)
{
CFArray result= CFArray (1);
result [0]= A [0] / M [0];
return result;
}
// check solvability
bool notDistinct= false;
for (int i= 0; i < r - 1; i++)
{
for (int j= i + 1; j < r; j++)
{
if (M [i] == M [j])
{
notDistinct= true;
break;
}
}
}
if (notDistinct)
return CFArray();
CanonicalForm master= 1;
Variable x= Variable (1);
for (int i= 0; i < r; i++)
master *= x - M [i];
CFList Pj;
CanonicalForm tmp;
for (int i= 0; i < r; i++)
{
tmp= master/(x - M [i]);
tmp /= tmp (M [i], 1);
Pj.append (tmp);
}
CFArray result= CFArray (r);
CFListIterator j= Pj;
for (int i= 1; i <= r; i++, j++)
{
tmp= 0;
for (int l= 0; l < A.size(); l++)
tmp += A[l]*j.getItem()[l];
result[i - 1]= tmp;
}
return result;
}
CFArray
solveGeneralVandermonde (const CFArray& M, const CFArray& A)
{
int r= M.size();
ASSERT (A.size() == r, "vector does not have right size");
if (r == 1)
{
CFArray result= CFArray (1);
result [0]= A[0] / M [0];
return result;
}
// check solvability
bool notDistinct= false;
for (int i= 0; i < r - 1; i++)
{
for (int j= i + 1; j < r; j++)
{
if (M [i] == M [j])
{
notDistinct= true;
break;
}
}
}
if (notDistinct)
return CFArray();
CanonicalForm master= 1;
Variable x= Variable (1);
for (int i= 0; i < r; i++)
master *= x - M [i];
master *= x;
CFList Pj;
CanonicalForm tmp;
for (int i= 0; i < r; i++)
{
tmp= master/(x - M [i]);
tmp /= tmp (M [i], 1);
Pj.append (tmp);
}
CFArray result= CFArray (r);
CFListIterator j= Pj;
for (int i= 1; i <= r; i++, j++)
{
tmp= 0;
for (int l= 1; l <= A.size(); l++)
tmp += A[l - 1]*j.getItem()[l];
result[i - 1]= tmp;
}
return result;
}
/// M in row echolon form, rk rank of M
CFArray
readOffSolution (const CFMatrix& M, const long rk)
{
CFArray result= CFArray (rk);
CanonicalForm tmp1, tmp2, tmp3;
for (int i= rk; i >= 1; i--)
{
tmp3= 0;
tmp1= M (i, M.columns());
for (int j= M.columns() - 1; j >= 1; j--)
{
tmp2= M (i, j);
if (j == i)
break;
else
tmp3 += tmp2*result[j - 1];
}
result[i - 1]= (tmp1 - tmp3)/tmp2;
}
return result;
}
CFArray
readOffSolution (const CFMatrix& M, const CFArray& L, const CFArray& partialSol)
{
CFArray result= CFArray (M.rows());
CanonicalForm tmp1, tmp2, tmp3;
int k;
for (int i= M.rows(); i >= 1; i--)
{
tmp3= 0;
tmp1= L[i - 1];
k= 0;
for (int j= M.columns(); j >= 1; j--, k++)
{
tmp2= M (i, j);
if (j == i)
break;
else
{
if (k > partialSol.size() - 1)
tmp3 += tmp2*result[j - 1];
else
tmp3 += tmp2*partialSol[partialSol.size() - k - 1];
}
}
result [i - 1]= (tmp1 - tmp3)/tmp2;
}
return result;
}
long
gaussianElimFp (CFMatrix& M, CFArray& L)
{
ASSERT (L.size() <= M.rows(), "dimension exceeded");
CFMatrix *N;
N= new CFMatrix (M.rows(), M.columns() + 1);
for (int i= 1; i <= M.rows(); i++)
for (int j= 1; j <= M.columns(); j++)
(*N) (i, j)= M (i, j);
int j= 1;
for (int i= 0; i < L.size(); i++, j++)
(*N) (j, M.columns() + 1)= L[i];
#ifdef HAVE_FLINT
nmod_mat_t FLINTN;
convertFacCFMatrix2nmod_mat_t (FLINTN, *N);
long rk= nmod_mat_rref (FLINTN);
delete N;
N= convertNmod_mat_t2FacCFMatrix (FLINTN);
nmod_mat_clear (FLINTN);
#else
int p= getCharacteristic ();
if (fac_NTL_char != p)
{
fac_NTL_char= p;
zz_p::init (p);
}
mat_zz_p *NTLN= convertFacCFMatrix2NTLmat_zz_p(*N);
delete N;
long rk= gauss (*NTLN);
N= convertNTLmat_zz_p2FacCFMatrix (*NTLN);
delete NTLN;
#endif
L= CFArray (M.rows());
for (int i= 0; i < M.rows(); i++)
L[i]= (*N) (i + 1, M.columns() + 1);
M= (*N) (1, M.rows(), 1, M.columns());
delete N;
return rk;
}
long
gaussianElimFq (CFMatrix& M, CFArray& L, const Variable& alpha)
{
ASSERT (L.size() <= M.rows(), "dimension exceeded");
CFMatrix *N;
N= new CFMatrix (M.rows(), M.columns() + 1);
for (int i= 1; i <= M.rows(); i++)
for (int j= 1; j <= M.columns(); j++)
(*N) (i, j)= M (i, j);
int j= 1;
for (int i= 0; i < L.size(); i++, j++)
(*N) (j, M.columns() + 1)= L[i];
#ifdef HAVE_FLINT
// convert mipo
nmod_poly_t mipo1;
convertFacCF2nmod_poly_t(mipo1,getMipo(alpha));
fq_nmod_ctx_t ctx;
fq_nmod_ctx_init_modulus(ctx,mipo1,"t");
nmod_poly_clear(mipo1);
// convert matrix
fq_nmod_mat_t FLINTN;
convertFacCFMatrix2Fq_nmod_mat_t (FLINTN, ctx, *N);
// rank
#if __FLINT_RELEASE >= 30100
long rk= fq_nmod_mat_rref (FLINTN,FLINTN,ctx);
#else
long rk= fq_nmod_mat_rref (FLINTN,ctx);
#endif
// clean up
fq_nmod_mat_clear (FLINTN,ctx);
fq_nmod_ctx_clear(ctx);
#elif defined(HAVE_NTL)
int p= getCharacteristic ();
if (fac_NTL_char != p)
{
fac_NTL_char= p;
zz_p::init (p);
}
zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
zz_pE::init (NTLMipo);
mat_zz_pE *NTLN= convertFacCFMatrix2NTLmat_zz_pE(*N);
long rk= gauss (*NTLN);
N= convertNTLmat_zz_pE2FacCFMatrix (*NTLN, alpha);
delete NTLN;
#else
factoryError("NTL/FLINT missing: gaussianElimFq");
#endif
M= (*N) (1, M.rows(), 1, M.columns());
L= CFArray (M.rows());
for (int i= 0; i < M.rows(); i++)
L[i]= (*N) (i + 1, M.columns() + 1);
delete N;
return rk;
}
CFArray
solveSystemFp (const CFMatrix& M, const CFArray& L)
{
ASSERT (L.size() <= M.rows(), "dimension exceeded");
CFMatrix *N;
N= new CFMatrix (M.rows(), M.columns() + 1);
for (int i= 1; i <= M.rows(); i++)
for (int j= 1; j <= M.columns(); j++)
(*N) (i, j)= M (i, j);
int j= 1;
for (int i= 0; i < L.size(); i++, j++)
(*N) (j, M.columns() + 1)= L[i];
#ifdef HAVE_FLINT
nmod_mat_t FLINTN;
convertFacCFMatrix2nmod_mat_t (FLINTN, *N);
long rk= nmod_mat_rref (FLINTN);
#else
int p= getCharacteristic ();
if (fac_NTL_char != p)
{
fac_NTL_char= p;
zz_p::init (p);
}
mat_zz_p *NTLN= convertFacCFMatrix2NTLmat_zz_p(*N);
long rk= gauss (*NTLN);
#endif
delete N;
if (rk != M.columns())
{
#ifdef HAVE_FLINT
nmod_mat_clear (FLINTN);
#else
delete NTLN;
#endif
return CFArray();
}
#ifdef HAVE_FLINT
N= convertNmod_mat_t2FacCFMatrix (FLINTN);
nmod_mat_clear (FLINTN);
#else
N= convertNTLmat_zz_p2FacCFMatrix (*NTLN);
delete NTLN;
#endif
CFArray A= readOffSolution (*N, rk);
delete N;
return A;
}
CFArray
solveSystemFq (const CFMatrix& M, const CFArray& L, const Variable& alpha)
{
ASSERT (L.size() <= M.rows(), "dimension exceeded");
CFMatrix *N;
N= new CFMatrix (M.rows(), M.columns() + 1);
for (int i= 1; i <= M.rows(); i++)
for (int j= 1; j <= M.columns(); j++)
(*N) (i, j)= M (i, j);
int j= 1;
for (int i= 0; i < L.size(); i++, j++)
(*N) (j, M.columns() + 1)= L[i];
#ifdef HAVE_FLINT
// convert mipo
nmod_poly_t mipo1;
convertFacCF2nmod_poly_t(mipo1,getMipo(alpha));
fq_nmod_ctx_t ctx;
fq_nmod_ctx_init_modulus(ctx,mipo1,"t");
nmod_poly_clear(mipo1);
// convert matrix
fq_nmod_mat_t FLINTN;
convertFacCFMatrix2Fq_nmod_mat_t (FLINTN, ctx, *N);
// rank
#if __FLINT_RELEASE >= 30100
long rk= fq_nmod_mat_rref (FLINTN,FLINTN,ctx);
#else
long rk= fq_nmod_mat_rref (FLINTN,ctx);
#endif
#elif defined(HAVE_NTL)
int p= getCharacteristic ();
if (fac_NTL_char != p)
{
fac_NTL_char= p;
zz_p::init (p);
}
zz_pX NTLMipo= convertFacCF2NTLzzpX (getMipo (alpha));
zz_pE::init (NTLMipo);
mat_zz_pE *NTLN= convertFacCFMatrix2NTLmat_zz_pE(*N);
long rk= gauss (*NTLN);
#else
factoryError("NTL/FLINT missing: solveSystemFq");
#endif
delete N;
if (rk != M.columns())
{
#if defined(HAVE_NTL) && !defined(HAVE_FLINT)
delete NTLN;
#endif
return CFArray();
}
#ifdef HAVE_FLINT
// convert and clean up
N=convertFq_nmod_mat_t2FacCFMatrix(FLINTN,ctx,alpha);
fq_nmod_mat_clear (FLINTN,ctx);
fq_nmod_ctx_clear(ctx);
#elif defined(HAVE_NTL)
N= convertNTLmat_zz_pE2FacCFMatrix (*NTLN, alpha);
delete NTLN;
#endif
CFArray A= readOffSolution (*N, rk);
delete N;
return A;
}
#endif
CFArray
getMonoms (const CanonicalForm& F)
{
if (F.inCoeffDomain())
{
CFArray result= CFArray (1);
result [0]= 1;
return result;
}
if (F.isUnivariate())
{
CFArray result= CFArray (size(F));
int j= 0;
for (CFIterator i= F; i.hasTerms(); i++, j++)
result[j]= power (F.mvar(), i.exp());
return result;
}
int numMon= size (F);
CFArray result= CFArray (numMon);
int j= 0;
CFArray recResult;
Variable x= F.mvar();
CanonicalForm powX;
for (CFIterator i= F; i.hasTerms(); i++)
{
powX= power (x, i.exp());
recResult= getMonoms (i.coeff());
for (int k= 0; k < recResult.size(); k++)
result[j+k]= powX*recResult[k];
j += recResult.size();
}
return result;
}
#if defined(HAVE_NTL) || defined(HAVE_FLINT)
CFArray
evaluateMonom (const CanonicalForm& F, const CFList& evalPoints)
{
if (F.inCoeffDomain())
{
CFArray result= CFArray (1);
result [0]= F;
return result;
}
if (F.isUnivariate())
{
ASSERT (evalPoints.length() == 1,
"expected an eval point with only one component");
CFArray result= CFArray (size(F));
int j= 0;
CanonicalForm evalPoint= evalPoints.getLast();
for (CFIterator i= F; i.hasTerms(); i++, j++)
result[j]= power (evalPoint, i.exp());
return result;
}
int numMon= size (F);
CFArray result= CFArray (numMon);
int j= 0;
CanonicalForm evalPoint= evalPoints.getLast();
CFList buf= evalPoints;
buf.removeLast();
CFArray recResult;
CanonicalForm powEvalPoint;
for (CFIterator i= F; i.hasTerms(); i++)
{
powEvalPoint= power (evalPoint, i.exp());
recResult= evaluateMonom (i.coeff(), buf);
for (int k= 0; k < recResult.size(); k++)
result[j+k]= powEvalPoint*recResult[k];
j += recResult.size();
}
return result;
}
CFArray
evaluate (const CFArray& A, const CFList& evalPoints)
{
CFArray result= A.size();
CanonicalForm tmp;
int k;
for (int i= 0; i < A.size(); i++)
{
tmp= A[i];
k= 1;
for (CFListIterator j= evalPoints; j.hasItem(); j++, k++)
tmp= tmp (j.getItem(), k);
result[i]= tmp;
}
return result;
}
CFList
evaluationPoints (const CanonicalForm& F, const CanonicalForm& G,
CanonicalForm& Feval, CanonicalForm& Geval,
const CanonicalForm& LCF, const bool& GF,
const Variable& alpha, bool& fail, CFList& list
)
{
int k= tmax (F.level(), G.level()) - 1;
Variable x= Variable (1);
CFList result;
FFRandom genFF;
GFRandom genGF;
int p= getCharacteristic ();
double bound;
if (alpha != Variable (1))
{
bound= pow ((double) p, (double) degree (getMipo(alpha)));
bound= pow (bound, (double) k);
}
else if (GF)
{
bound= pow ((double) p, (double) getGFDegree());
bound= pow ((double) bound, (double) k);
}
else
bound= pow ((double) p, (double) k);
CanonicalForm random;
int j;
bool zeroOneOccured= false;
bool allEqual= false;
CanonicalForm buf;
do
{
random= 0;
// possible overflow if list.length() does not fit into a int
if (list.length() >= bound)
{
fail= true;
break;
}
for (int i= 0; i < k; i++)
{
if (GF)
{
result.append (genGF.generate());
random += result.getLast()*power (x, i);
}
else if (alpha.level() != 1)
{
AlgExtRandomF genAlgExt (alpha);
result.append (genAlgExt.generate());
random += result.getLast()*power (x, i);
}
else
{
result.append (genFF.generate());
random += result.getLast()*power (x, i);
}
if (result.getLast().isOne() || result.getLast().isZero())
zeroOneOccured= true;
}
if (find (list, random))
{
zeroOneOccured= false;
allEqual= false;
result= CFList();
continue;
}
if (zeroOneOccured)
{
list.append (random);
zeroOneOccured= false;
allEqual= false;
result= CFList();
continue;
}
// no zero at this point
if (k > 1)
{
allEqual= true;
CFIterator iter= random;
buf= iter.coeff();
iter++;
for (; iter.hasTerms(); iter++)
if (buf != iter.coeff())
allEqual= false;
}
if (allEqual)
{
list.append (random);
allEqual= false;
zeroOneOccured= false;
result= CFList();
continue;
}
Feval= F;
Geval= G;
CanonicalForm LCeval= LCF;
j= 1;
for (CFListIterator i= result; i.hasItem(); i++, j++)
{
Feval= Feval (i.getItem(), j);
Geval= Geval (i.getItem(), j);
LCeval= LCeval (i.getItem(), j);
}
if (LCeval.isZero())
{
if (!find (list, random))
list.append (random);
zeroOneOccured= false;
allEqual= false;
result= CFList();
continue;
}
if (list.length() >= bound)
{
fail= true;
break;
}
} while (find (list, random));
return result;
}
/// multiply two lists componentwise
void mult (CFList& L1, const CFList& L2)
{
ASSERT (L1.length() == L2.length(), "lists of the same size expected");
CFListIterator j= L2;
for (CFListIterator i= L1; i.hasItem(); i++, j++)
i.getItem() *= j.getItem();
}
void eval (const CanonicalForm& A, const CanonicalForm& B, CanonicalForm& Aeval,
CanonicalForm& Beval, const CFList& L)
{
Aeval= A;
Beval= B;
int j= 1;
for (CFListIterator i= L; i.hasItem(); i++, j++)
{
Aeval= Aeval (i.getItem(), j);
Beval= Beval (i.getItem(), j);
}
}
CanonicalForm
monicSparseInterpol (const CanonicalForm& F, const CanonicalForm& G,
const CanonicalForm& skeleton, const Variable& alpha,
bool& fail, CFArray*& coeffMonoms, CFArray& Monoms
)
{
CanonicalForm A= F;
CanonicalForm B= G;
if (F.isZero() && degree(G) > 0) return G/Lc(G);
else if (G.isZero() && degree (F) > 0) return F/Lc(F);
else if (F.isZero() && G.isZero()) return F.genOne();
if (F.inBaseDomain() || G.inBaseDomain()) return F.genOne();
if (F.isUnivariate() && fdivides(F, G)) return F/Lc(F);
if (G.isUnivariate() && fdivides(G, F)) return G/Lc(G);
if (F == G) return F/Lc(F);
ASSERT (degree (A, 1) == 0, "expected degree (F, 1) == 0");
ASSERT (degree (B, 1) == 0, "expected degree (G, 1) == 0");
CFMap M,N;
int best_level= myCompress (A, B, M, N, false);
if (best_level == 0)
return B.genOne();
A= M(A);
B= M(B);
Variable x= Variable (1);
//univariate case
if (A.isUnivariate() && B.isUnivariate())
return N (gcd (A, B));
CanonicalForm skel= M(skeleton);
CanonicalForm cA, cB; // content of A and B
CanonicalForm ppA, ppB; // primitive part of A and B
CanonicalForm gcdcAcB;
cA = uni_content (A);
cB = uni_content (B);
gcdcAcB= gcd (cA, cB);
ppA= A/cA;
ppB= B/cB;
CanonicalForm lcA, lcB; // leading coefficients of A and B
CanonicalForm gcdlcAlcB;
lcA= uni_lcoeff (ppA);
lcB= uni_lcoeff (ppB);
if (fdivides (lcA, lcB))
{
if (fdivides (A, B))
return F/Lc(F);
}
if (fdivides (lcB, lcA))
{
if (fdivides (B, A))
return G/Lc(G);
}
gcdlcAlcB= gcd (lcA, lcB);
int skelSize= size (skel, skel.mvar());
int j= 0;
int biggestSize= 0;
for (CFIterator i= skel; i.hasTerms(); i++, j++)
biggestSize= tmax (biggestSize, size (i.coeff()));
CanonicalForm g, Aeval, Beval;
CFList evalPoints;
bool evalFail= false;
CFList list;
bool GF= false;
CanonicalForm LCA= LC (A);
CanonicalForm tmp;
CFArray gcds= CFArray (biggestSize);
CFList * pEvalPoints= new CFList [biggestSize];
Variable V_buf= alpha, V_buf4= alpha;
CFList source, dest;
CanonicalForm prim_elem, im_prim_elem;
CanonicalForm prim_elem_alpha, im_prim_elem_alpha;
for (int i= 0; i < biggestSize; i++)
{
if (i == 0)
evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf, evalFail,
list);
else
{
mult (evalPoints, pEvalPoints [0]);
eval (A, B, Aeval, Beval, evalPoints);
}
if (evalFail)
{
if (V_buf.level() != 1)
{
do
{
Variable V_buf3= V_buf;
V_buf= chooseExtension (V_buf);
source= CFList();
dest= CFList();
bool prim_fail= false;
Variable V_buf2;
prim_elem= primitiveElement (V_buf4, V_buf2, prim_fail);
if (V_buf4 == alpha && alpha.level() != 1)
prim_elem_alpha= prim_elem;
ASSERT (!prim_fail, "failure in integer factorizer");
if (prim_fail)
; //ERROR
else
im_prim_elem= mapPrimElem (prim_elem, V_buf4, V_buf);
DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf));
DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf2));
if (V_buf4 == alpha && alpha.level() != 1)
im_prim_elem_alpha= im_prim_elem;
else if (alpha.level() != 1)
im_prim_elem_alpha= mapUp (im_prim_elem_alpha, V_buf4, V_buf,
prim_elem, im_prim_elem, source, dest);
for (CFListIterator j= list; j.hasItem(); j++)
j.getItem()= mapUp (j.getItem(), V_buf4, V_buf, prim_elem,
im_prim_elem, source, dest);
for (int k= 0; k < i; k++)
{
for (CFListIterator j= pEvalPoints[k]; j.hasItem(); j++)
j.getItem()= mapUp (j.getItem(), V_buf4, V_buf, prim_elem,
im_prim_elem, source, dest);
gcds[k]= mapUp (gcds[k], V_buf4, V_buf, prim_elem, im_prim_elem,
source, dest);
}
if (alpha.level() != 1)
{
A= mapUp (A, V_buf4, V_buf, prim_elem, im_prim_elem, source,dest);
B= mapUp (B, V_buf4, V_buf, prim_elem, im_prim_elem, source,dest);
}
V_buf4= V_buf;
evalFail= false;
evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf,
evalFail, list);
} while (evalFail);
}
else
{
CanonicalForm mipo;
int deg= 2;
bool initialized= false;
do
{
mipo= randomIrredpoly (deg, x);
if (initialized)
setMipo (V_buf, mipo);
else
V_buf= rootOf (mipo);
evalFail= false;
initialized= true;
evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf,
evalFail, list);
deg++;
} while (evalFail);
V_buf4= V_buf;
}
}
g= gcd (Aeval, Beval);
g /= Lc (g);
if (degree (g) != degree (skel) || (skelSize != size (g)))
{
delete[] pEvalPoints;
fail= true;
if (alpha.level() != 1 && V_buf != alpha)
prune1 (alpha);
return 0;
}
CFIterator l= skel;
for (CFIterator k= g; k.hasTerms(); k++, l++)
{
if (k.exp() != l.exp())
{
delete[] pEvalPoints;
fail= true;
if (alpha.level() != 1 && V_buf != alpha)
prune1 (alpha);
return 0;
}
}
pEvalPoints[i]= evalPoints;
gcds[i]= g;
tmp= 0;
int j= 0;
for (CFListIterator k= evalPoints; k.hasItem(); k++, j++)
tmp += k.getItem()*power (x, j);
list.append (tmp);
}
if (Monoms.size() == 0)
Monoms= getMonoms (skel);
coeffMonoms= new CFArray [skelSize];
j= 0;
for (CFIterator i= skel; i.hasTerms(); i++, j++)
coeffMonoms[j]= getMonoms (i.coeff());
CFArray* pL= new CFArray [skelSize];
CFArray* pM= new CFArray [skelSize];
for (int i= 0; i < biggestSize; i++)
{
CFIterator l= gcds [i];
evalPoints= pEvalPoints [i];
for (int k= 0; k < skelSize; k++, l++)
{
if (i == 0)
pL[k]= CFArray (biggestSize);
pL[k] [i]= l.coeff();
if (i == 0)
pM[k]= evaluate (coeffMonoms [k], evalPoints);
}
}
CFArray solution;
CanonicalForm result= 0;
int ind= 0;
CFArray bufArray;
CFMatrix Mat;
for (int k= 0; k < skelSize; k++)
{
if (biggestSize != coeffMonoms[k].size())
{
bufArray= CFArray (coeffMonoms[k].size());
for (int i= 0; i < coeffMonoms[k].size(); i++)
bufArray [i]= pL[k] [i];
solution= solveGeneralVandermonde (pM [k], bufArray);
}
else
solution= solveGeneralVandermonde (pM [k], pL[k]);
if (solution.size() == 0)
{
delete[] pEvalPoints;
delete[] pM;
delete[] pL;
delete[] coeffMonoms;
fail= true;
if (alpha.level() != 1 && V_buf != alpha)
prune1 (alpha);
return 0;
}
for (int l= 0; l < solution.size(); l++)
result += solution[l]*Monoms [ind + l];
ind += solution.size();
}
delete[] pEvalPoints;
delete[] pM;
delete[] pL;
delete[] coeffMonoms;
if (alpha.level() != 1 && V_buf != alpha)
{
CFList u, v;
result= mapDown (result, prim_elem_alpha, im_prim_elem_alpha, alpha, u, v);
prune1 (alpha);
}
result= N(result);
if (fdivides (result, F) && fdivides (result, G))
return result;
else
{
fail= true;
return 0;
}
}
CanonicalForm
nonMonicSparseInterpol (const CanonicalForm& F, const CanonicalForm& G,
const CanonicalForm& skeleton, const Variable& alpha,
bool& fail, CFArray*& coeffMonoms, CFArray& Monoms
)
{
CanonicalForm A= F;
CanonicalForm B= G;
if (F.isZero() && degree(G) > 0) return G/Lc(G);
else if (G.isZero() && degree (F) > 0) return F/Lc(F);
else if (F.isZero() && G.isZero()) return F.genOne();
if (F.inBaseDomain() || G.inBaseDomain()) return F.genOne();
if (F.isUnivariate() && fdivides(F, G)) return F/Lc(F);
if (G.isUnivariate() && fdivides(G, F)) return G/Lc(G);
if (F == G) return F/Lc(F);
ASSERT (degree (A, 1) == 0, "expected degree (F, 1) == 0");
ASSERT (degree (B, 1) == 0, "expected degree (G, 1) == 0");
CFMap M,N;
int best_level= myCompress (A, B, M, N, false);
if (best_level == 0)
return B.genOne();
A= M(A);
B= M(B);
Variable x= Variable (1);
//univariate case
if (A.isUnivariate() && B.isUnivariate())
return N (gcd (A, B));
CanonicalForm skel= M(skeleton);
CanonicalForm cA, cB; // content of A and B
CanonicalForm ppA, ppB; // primitive part of A and B
CanonicalForm gcdcAcB;
cA = uni_content (A);
cB = uni_content (B);
gcdcAcB= gcd (cA, cB);
ppA= A/cA;
ppB= B/cB;
CanonicalForm lcA, lcB; // leading coefficients of A and B
CanonicalForm gcdlcAlcB;
lcA= uni_lcoeff (ppA);
lcB= uni_lcoeff (ppB);
if (fdivides (lcA, lcB))
{
if (fdivides (A, B))
return F/Lc(F);
}
if (fdivides (lcB, lcA))
{
if (fdivides (B, A))
return G/Lc(G);
}
gcdlcAlcB= gcd (lcA, lcB);
int skelSize= size (skel, skel.mvar());
int j= 0;
int biggestSize= 0;
int bufSize;
int numberUni= 0;
for (CFIterator i= skel; i.hasTerms(); i++, j++)
{
bufSize= size (i.coeff());
biggestSize= tmax (biggestSize, bufSize);
numberUni += bufSize;
}
numberUni--;
numberUni= (int) ceil ( (double) numberUni / (skelSize - 1));
biggestSize= tmax (biggestSize , numberUni);
numberUni= biggestSize + size (LC(skel)) - 1;
int biggestSize2= tmax (numberUni, biggestSize);
CanonicalForm g, Aeval, Beval;
CFList evalPoints;
CFArray coeffEval;
bool evalFail= false;
CFList list;
bool GF= false;
CanonicalForm LCA= LC (A);
CanonicalForm tmp;
CFArray gcds= CFArray (biggestSize);
CFList * pEvalPoints= new CFList [biggestSize];
Variable V_buf= alpha, V_buf4= alpha;
CFList source, dest;
CanonicalForm prim_elem, im_prim_elem;
CanonicalForm prim_elem_alpha, im_prim_elem_alpha;
for (int i= 0; i < biggestSize; i++)
{
if (i == 0)
{
if (getCharacteristic() > 3)
evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf,
evalFail, list);
else
evalFail= true;
if (evalFail)
{
if (V_buf.level() != 1)
{
do
{
Variable V_buf3= V_buf;
V_buf= chooseExtension (V_buf);
source= CFList();
dest= CFList();
bool prim_fail= false;
Variable V_buf2;
prim_elem= primitiveElement (V_buf4, V_buf2, prim_fail);
if (V_buf4 == alpha && alpha.level() != 1)
prim_elem_alpha= prim_elem;
ASSERT (!prim_fail, "failure in integer factorizer");
if (prim_fail)
; //ERROR
else
im_prim_elem= mapPrimElem (prim_elem, V_buf4, V_buf);
DEBOUTLN (cerr, "getMipo (V_buf)= " << getMipo (V_buf));
DEBOUTLN (cerr, "getMipo (V_buf2)= " << getMipo (V_buf2));
if (V_buf4 == alpha && alpha.level() != 1)
im_prim_elem_alpha= im_prim_elem;
else if (alpha.level() != 1)
im_prim_elem_alpha= mapUp (im_prim_elem_alpha, V_buf4, V_buf,
prim_elem, im_prim_elem, source, dest);
for (CFListIterator i= list; i.hasItem(); i++)
i.getItem()= mapUp (i.getItem(), V_buf4, V_buf, prim_elem,
im_prim_elem, source, dest);
if (alpha.level() != 1)
{
A= mapUp (A, V_buf4, V_buf, prim_elem, im_prim_elem, source,dest);
B= mapUp (B, V_buf4, V_buf, prim_elem, im_prim_elem, source,dest);
}
evalFail= false;
V_buf4= V_buf;
evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf,
evalFail, list);
} while (evalFail);
}
else
{
CanonicalForm mipo;
int deg= 2;
bool initialized= false;
do
{
mipo= randomIrredpoly (deg, x);
if (initialized)
setMipo (V_buf, mipo);
else
V_buf= rootOf (mipo);
evalFail= false;
initialized= true;
evalPoints= evaluationPoints (A, B, Aeval, Beval, LCA, GF, V_buf,
evalFail, list);
deg++;
} while (evalFail);
V_buf4= V_buf;
}
}
}
else
{
mult (evalPoints, pEvalPoints[0]);
eval (A, B, Aeval, Beval, evalPoints);
}
g= gcd (Aeval, Beval);
g /= Lc (g);
if (degree (g) != degree (skel) || (skelSize != size (g)))
{
delete[] pEvalPoints;
fail= true;
if (alpha.level() != 1 && V_buf != alpha)
prune1 (alpha);
return 0;
}
CFIterator l= skel;
for (CFIterator k= g; k.hasTerms(); k++, l++)
{
if (k.exp() != l.exp())
{
delete[] pEvalPoints;
fail= true;
if (alpha.level() != 1 && V_buf != alpha)
prune1 (alpha);
return 0;
}
}
pEvalPoints[i]= evalPoints;
gcds[i]= g;
tmp= 0;
int j= 0;
for (CFListIterator k= evalPoints; k.hasItem(); k++, j++)
tmp += k.getItem()*power (x, j);
list.append (tmp);
}
if (Monoms.size() == 0)
Monoms= getMonoms (skel);
coeffMonoms= new CFArray [skelSize];
j= 0;
for (CFIterator i= skel; i.hasTerms(); i++, j++)
coeffMonoms[j]= getMonoms (i.coeff());
int minimalColumnsIndex;
if (skelSize > 1)
minimalColumnsIndex= 1;
else
minimalColumnsIndex= 0;
int minimalColumns=-1;
CFArray* pM= new CFArray [skelSize];
CFMatrix Mat;
// find the Matrix with minimal number of columns
for (int i= 0; i < skelSize; i++)
{
pM[i]= CFArray (coeffMonoms[i].size());
if (i == 1)
minimalColumns= coeffMonoms[i].size();
if (i > 1)
{
if (minimalColumns > coeffMonoms[i].size())
{
minimalColumns= coeffMonoms[i].size();
minimalColumnsIndex= i;
}
}
}
CFMatrix* pMat= new CFMatrix [2];
pMat[0]= CFMatrix (biggestSize, coeffMonoms[0].size());
pMat[1]= CFMatrix (biggestSize, coeffMonoms[minimalColumnsIndex].size());
CFArray* pL= new CFArray [skelSize];
for (int i= 0; i < biggestSize; i++)
{
CFIterator l= gcds [i];
evalPoints= pEvalPoints [i];
for (int k= 0; k < skelSize; k++, l++)
{
if (i == 0)
pL[k]= CFArray (biggestSize);
pL[k] [i]= l.coeff();
if (i == 0 && k != 0 && k != minimalColumnsIndex)
pM[k]= evaluate (coeffMonoms[k], evalPoints);
else if (i == 0 && (k == 0 || k == minimalColumnsIndex))
coeffEval= evaluate (coeffMonoms[k], evalPoints);
if (i > 0 && (k == 0 || k == minimalColumnsIndex))
coeffEval= evaluate (coeffMonoms[k], evalPoints);
if (k == 0)
{
if (pMat[k].rows() >= i + 1)
{
for (int ind= 1; ind < coeffEval.size() + 1; ind++)
pMat[k] (i + 1, ind)= coeffEval[ind - 1];
}
}
if (k == minimalColumnsIndex)
{
if (pMat[1].rows() >= i + 1)
{
for (int ind= 1; ind < coeffEval.size() + 1; ind++)
pMat[1] (i + 1, ind)= coeffEval[ind - 1];
}
}
}
}
CFArray solution;
CanonicalForm result= 0;
int ind= 1;
int matRows, matColumns;
matRows= pMat[1].rows();
matColumns= pMat[0].columns() - 1;
matColumns += pMat[1].columns();
Mat= CFMatrix (matRows, matColumns);
for (int i= 1; i <= matRows; i++)
for (int j= 1; j <= pMat[1].columns(); j++)
Mat (i, j)= pMat[1] (i, j);
for (int j= pMat[1].columns() + 1; j <= matColumns;
j++, ind++)
{
for (int i= 1; i <= matRows; i++)
Mat (i, j)= (-pMat [0] (i, ind + 1))*pL[minimalColumnsIndex] [i - 1];
}
CFArray firstColumn= CFArray (pMat[0].rows());
for (int i= 0; i < pMat[0].rows(); i++)
firstColumn [i]= pMat[0] (i + 1, 1);
CFArray bufArray= pL[minimalColumnsIndex];
for (int i= 0; i < biggestSize; i++)
pL[minimalColumnsIndex] [i] *= firstColumn[i];
CFMatrix bufMat= pMat[1];
pMat[1]= Mat;
if (V_buf.level() != 1)
solution= solveSystemFq (pMat[1],
pL[minimalColumnsIndex], V_buf);
else
solution= solveSystemFp (pMat[1],
pL[minimalColumnsIndex]);
if (solution.size() == 0)
{ //Javadi's method failed, try deKleine, Monagan, Wittkopf instead
CFMatrix bufMat0= pMat[0];
delete [] pMat;
pMat= new CFMatrix [skelSize];
pL[minimalColumnsIndex]= bufArray;
CFList* bufpEvalPoints= NULL;
CFArray bufGcds;
if (biggestSize != biggestSize2)
{
bufpEvalPoints= pEvalPoints;
pEvalPoints= new CFList [biggestSize2];
bufGcds= gcds;
gcds= CFArray (biggestSize2);
for (int i= 0; i < biggestSize; i++)
{
pEvalPoints[i]= bufpEvalPoints [i];
gcds[i]= bufGcds[i];
}
for (int i= 0; i < biggestSize2 - biggestSize; i++)
{
mult (evalPoints, pEvalPoints[0]);
eval (A, B, Aeval, Beval, evalPoints);
g= gcd (Aeval, Beval);
g /= Lc (g);
if (degree (g) != degree (skel) || (skelSize != size (g)))
{
delete[] pEvalPoints;
delete[] pMat;
delete[] pL;
delete[] coeffMonoms;
delete[] pM;
if (bufpEvalPoints != NULL)
delete [] bufpEvalPoints;
fail= true;
if (alpha.level() != 1 && V_buf != alpha)
prune1 (alpha);
return 0;
}
CFIterator l= skel;
for (CFIterator k= g; k.hasTerms(); k++, l++)
{
if (k.exp() != l.exp())
{
delete[] pEvalPoints;
delete[] pMat;
delete[] pL;
delete[] coeffMonoms;
delete[] pM;
if (bufpEvalPoints != NULL)
delete [] bufpEvalPoints;
fail= true;
if (alpha.level() != 1 && V_buf != alpha)
prune1 (alpha);
return 0;
}
}
pEvalPoints[i + biggestSize]= evalPoints;
gcds[i + biggestSize]= g;
}
}
for (int i= 0; i < biggestSize; i++)
{
CFIterator l= gcds [i];
evalPoints= pEvalPoints [i];
for (int k= 1; k < skelSize; k++, l++)
{
if (i == 0)
pMat[k]= CFMatrix (biggestSize2,coeffMonoms[k].size()+biggestSize2-1);
if (k == minimalColumnsIndex)
continue;
coeffEval= evaluate (coeffMonoms[k], evalPoints);
if (pMat[k].rows() >= i + 1)
{
for (int ind= 1; ind < coeffEval.size() + 1; ind++)
pMat[k] (i + 1, ind)= coeffEval[ind - 1];
}
}
}
Mat= bufMat0;
pMat[0]= CFMatrix (biggestSize2, coeffMonoms[0].size() + biggestSize2 - 1);
for (int j= 1; j <= Mat.rows(); j++)
for (int k= 1; k <= Mat.columns(); k++)
pMat [0] (j,k)= Mat (j,k);
Mat= bufMat;
for (int j= 1; j <= Mat.rows(); j++)
for (int k= 1; k <= Mat.columns(); k++)
pMat [minimalColumnsIndex] (j,k)= Mat (j,k);
// write old matrix entries into new matrices
for (int i= 0; i < skelSize; i++)
{
bufArray= pL[i];
pL[i]= CFArray (biggestSize2);
for (int j= 0; j < bufArray.size(); j++)
pL[i] [j]= bufArray [j];
}
//write old vector entries into new and add new entries to old matrices
for (int i= 0; i < biggestSize2 - biggestSize; i++)
{
CFIterator l= gcds [i + biggestSize];
evalPoints= pEvalPoints [i + biggestSize];
for (int k= 0; k < skelSize; k++, l++)
{
pL[k] [i + biggestSize]= l.coeff();
coeffEval= evaluate (coeffMonoms[k], evalPoints);
if (pMat[k].rows() >= i + biggestSize + 1)
{
for (int ind= 1; ind < coeffEval.size() + 1; ind++)
pMat[k] (i + biggestSize + 1, ind)= coeffEval[ind - 1];
}
}
}
// begin new
for (int i= 0; i < skelSize; i++)
{
if (pL[i].size() > 1)
{
for (int j= 2; j <= pMat[i].rows(); j++)
pMat[i] (j, coeffMonoms[i].size() + j - 1)=
-pL[i] [j - 1];
}
}
matColumns= biggestSize2 - 1;
matRows= 0;
for (int i= 0; i < skelSize; i++)
{
if (V_buf.level() == 1)
(void) gaussianElimFp (pMat[i], pL[i]);
else
(void) gaussianElimFq (pMat[i], pL[i], V_buf);
if (pMat[i] (coeffMonoms[i].size(), coeffMonoms[i].size()) == 0)
{
delete[] pEvalPoints;
delete[] pMat;
delete[] pL;
delete[] coeffMonoms;
delete[] pM;
if (bufpEvalPoints != NULL)
delete [] bufpEvalPoints;
fail= true;
if (alpha.level() != 1 && V_buf != alpha)
prune1 (alpha);
return 0;
}
matRows += pMat[i].rows() - coeffMonoms[i].size();
}
CFMatrix bufMat;
Mat= CFMatrix (matRows, matColumns);
ind= 0;
bufArray= CFArray (matRows);
CFArray bufArray2;
for (int i= 0; i < skelSize; i++)
{
if (coeffMonoms[i].size() + 1 >= pMat[i].rows() || coeffMonoms[i].size() + 1 >= pMat[i].columns())
{
delete[] pEvalPoints;
delete[] pMat;
delete[] pL;
delete[] coeffMonoms;
delete[] pM;
if (bufpEvalPoints != NULL)
delete [] bufpEvalPoints;
fail= true;
if (alpha.level() != 1 && V_buf != alpha)
prune1 (alpha);
return 0;
}
bufMat= pMat[i] (coeffMonoms[i].size() + 1, pMat[i].rows(),
coeffMonoms[i].size() + 1, pMat[i].columns());
for (int j= 1; j <= bufMat.rows(); j++)
for (int k= 1; k <= bufMat.columns(); k++)
Mat (j + ind, k)= bufMat(j, k);
bufArray2= coeffMonoms[i].size();
for (int j= 1; j <= pMat[i].rows(); j++)
{
if (j > coeffMonoms[i].size())
bufArray [j-coeffMonoms[i].size() + ind - 1]= pL[i] [j - 1];
else
bufArray2 [j - 1]= pL[i] [j - 1];
}
pL[i]= bufArray2;
ind += bufMat.rows();
pMat[i]= pMat[i] (1, coeffMonoms[i].size(), 1, pMat[i].columns());
}
if (V_buf.level() != 1)
solution= solveSystemFq (Mat, bufArray, V_buf);
else
solution= solveSystemFp (Mat, bufArray);
if (solution.size() == 0)
{
delete[] pEvalPoints;
delete[] pMat;
delete[] pL;
delete[] coeffMonoms;
delete[] pM;
if (bufpEvalPoints != NULL)
delete [] bufpEvalPoints;
fail= true;
if (alpha.level() != 1 && V_buf != alpha)
prune1 (alpha);
return 0;
}
ind= 0;
result= 0;
CFArray bufSolution;
for (int i= 0; i < skelSize; i++)
{
bufSolution= readOffSolution (pMat[i], pL[i], solution);
for (int i= 0; i < bufSolution.size(); i++, ind++)
result += Monoms [ind]*bufSolution[i];
}
if (alpha.level() != 1 && V_buf != alpha)
{
CFList u, v;
result= mapDown (result,prim_elem_alpha, im_prim_elem_alpha, alpha, u, v);
prune1 (alpha);
}
result= N(result);
delete[] pEvalPoints;
delete[] pMat;
delete[] pL;
delete[] coeffMonoms;
delete[] pM;
if (bufpEvalPoints != NULL)
delete [] bufpEvalPoints;
if (fdivides (result, F) && fdivides (result, G))
return result;
else
{
fail= true;
return 0;
}
} // end of deKleine, Monagan & Wittkopf
result += Monoms[0];
int ind2= 0, ind3= 2;
ind= 0;
for (int l= 0; l < minimalColumnsIndex; l++)
ind += coeffMonoms[l].size();
for (int l= solution.size() - pMat[0].columns() + 1; l < solution.size();
l++, ind2++, ind3++)
{
result += solution[l]*Monoms [1 + ind2];
for (int i= 0; i < pMat[0].rows(); i++)
firstColumn[i] += solution[l]*pMat[0] (i + 1, ind3);
}
for (int l= 0; l < solution.size() + 1 - pMat[0].columns(); l++)
result += solution[l]*Monoms [ind + l];
ind= coeffMonoms[0].size();
for (int k= 1; k < skelSize; k++)
{
if (k == minimalColumnsIndex)
{
ind += coeffMonoms[k].size();
continue;
}
if (k != minimalColumnsIndex)
{
for (int i= 0; i < biggestSize; i++)
pL[k] [i] *= firstColumn [i];
}
if (biggestSize != coeffMonoms[k].size() && k != minimalColumnsIndex)
{
bufArray= CFArray (coeffMonoms[k].size());
for (int i= 0; i < bufArray.size(); i++)
bufArray [i]= pL[k] [i];
solution= solveGeneralVandermonde (pM[k], bufArray);
}
else
solution= solveGeneralVandermonde (pM[k], pL[k]);
if (solution.size() == 0)
{
delete[] pEvalPoints;
delete[] pMat;
delete[] pL;
delete[] coeffMonoms;
delete[] pM;
fail= true;
if (alpha.level() != 1 && V_buf != alpha)
prune1 (alpha);
return 0;
}
if (k != minimalColumnsIndex)
{
for (int l= 0; l < solution.size(); l++)
result += solution[l]*Monoms [ind + l];
ind += solution.size();
}
}
delete[] pEvalPoints;
delete[] pMat;
delete[] pL;
delete[] pM;
delete[] coeffMonoms;
if (alpha.level() != 1 && V_buf != alpha)
{
CFList u, v;
result= mapDown (result, prim_elem, im_prim_elem, alpha, u, v);
prune1 (alpha);
}
result= N(result);
if (fdivides (result, F) && fdivides (result, G))
return result;
else
{
fail= true;
return 0;
}
}
CanonicalForm sparseGCDFq (const CanonicalForm& F, const CanonicalForm& G,
const Variable & alpha, CFList& l, bool& topLevel)
{
CanonicalForm A= F;
CanonicalForm B= G;
if (F.isZero() && degree(G) > 0) return G/Lc(G);
else if (G.isZero() && degree (F) > 0) return F/Lc(F);
else if (F.isZero() && G.isZero()) return F.genOne();
if (F.inBaseDomain() || G.inBaseDomain()) return F.genOne();
if (F.isUnivariate() && fdivides(F, G)) return F/Lc(F);
if (G.isUnivariate() && fdivides(G, F)) return G/Lc(G);
if (F == G) return F/Lc(F);
CFMap M,N;
int best_level= myCompress (A, B, M, N, topLevel);
if (best_level == 0) return B.genOne();
A= M(A);
B= M(B);
Variable x= Variable (1);
//univariate case
if (A.isUnivariate() && B.isUnivariate())
return N (gcd (A, B));
CanonicalForm cA, cB; // content of A and B
CanonicalForm ppA, ppB; // primitive part of A and B
CanonicalForm gcdcAcB;
cA = uni_content (A);
cB = uni_content (B);
gcdcAcB= gcd (cA, cB);
ppA= A/cA;
ppB= B/cB;
CanonicalForm lcA, lcB; // leading coefficients of A and B
CanonicalForm gcdlcAlcB;
lcA= uni_lcoeff (ppA);
lcB= uni_lcoeff (ppB);
if (fdivides (lcA, lcB))
{
if (fdivides (A, B))
return F/Lc(F);
}
if (fdivides (lcB, lcA))
{
if (fdivides (B, A))
return G/Lc(G);
}
gcdlcAlcB= gcd (lcA, lcB);
int d= totaldegree (ppA, Variable (2), Variable (ppA.level()));
int d0;
if (d == 0)
return N(gcdcAcB);
d0= totaldegree (ppB, Variable (2), Variable (ppB.level()));
if (d0 < d)
d= d0;
if (d == 0)
return N(gcdcAcB);
CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH, skeleton;
CanonicalForm newtonPoly= 1;
m= gcdlcAlcB;
G_m= 0;
H= 0;
bool fail= false;
topLevel= false;
bool inextension= false;
Variable V_buf= alpha, V_buf4= alpha;
CanonicalForm prim_elem, im_prim_elem;
CanonicalForm prim_elem_alpha, im_prim_elem_alpha;
CFList source, dest;
do // first do
{
random_element= randomElement (m, V_buf, l, fail);
if (random_element == 0 && !fail)
{
if (!find (l, random_element))
l.append (random_element);
continue;
}
if (fail)
{
source= CFList();
dest= CFList();
Variable V_buf3= V_buf;
V_buf= chooseExtension (V_buf);
bool prim_fail= false;
Variable V_buf2;
prim_elem= primitiveElement (V_buf4, V_buf2, prim_fail);
if (V_buf4 == alpha)
prim_elem_alpha= prim_elem;
if (V_buf3 != V_buf4)
{
m= mapDown (m, prim_elem, im_prim_elem, V_buf4, source, dest);
G_m= mapDown (m, prim_elem, im_prim_elem, V_buf4, source, dest);
newtonPoly= mapDown (newtonPoly, prim_elem, im_prim_elem, V_buf4,
source, dest);
ppA= mapDown (ppA, prim_elem, im_prim_elem, V_buf4, source, dest);
ppB= mapDown (ppB, prim_elem, im_prim_elem, V_buf4, source, dest);
gcdlcAlcB= mapDown (gcdlcAlcB, prim_elem, im_prim_elem, V_buf4, source,
dest);
for (CFListIterator i= l; i.hasItem(); i++)
i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, V_buf4,
source, dest);
}
ASSERT (!prim_fail, "failure in integer factorizer");
if (prim_fail)
; //ERROR
else
im_prim_elem= mapPrimElem (prim_elem, V_buf4, V_buf);
if (V_buf4 == alpha)
im_prim_elem_alpha= im_prim_elem;
else
im_prim_elem_alpha= mapUp (im_prim_elem_alpha, V_buf4, V_buf, prim_elem,
im_prim_elem, source, dest);
DEBOUTLN (cerr, "getMipo (V_buf4)= " << getMipo (V_buf4));
DEBOUTLN (cerr, "getMipo (V_buf2)= " << getMipo (V_buf2));
inextension= true;
for (CFListIterator i= l; i.hasItem(); i++)
i.getItem()= mapUp (i.getItem(), V_buf4, V_buf, prim_elem,
im_prim_elem, source, dest);
m= mapUp (m, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
G_m= mapUp (G_m, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
newtonPoly= mapUp (newtonPoly, V_buf4, V_buf, prim_elem, im_prim_elem,
source, dest);
ppA= mapUp (ppA, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
ppB= mapUp (ppB, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
gcdlcAlcB= mapUp (gcdlcAlcB, V_buf4, V_buf, prim_elem, im_prim_elem,
source, dest);
fail= false;
random_element= randomElement (m, V_buf, l, fail );
DEBOUTLN (cerr, "fail= " << fail);
CFList list;
TIMING_START (gcd_recursion);
G_random_element=
sparseGCDFq (ppA (random_element, x), ppB (random_element, x), V_buf,
list, topLevel);
TIMING_END_AND_PRINT (gcd_recursion,
"time for recursive call: ");
DEBOUTLN (cerr, "G_random_element= " << G_random_element);
V_buf4= V_buf;
}
else
{
CFList list;
TIMING_START (gcd_recursion);
G_random_element=
sparseGCDFq (ppA(random_element, x), ppB(random_element, x), V_buf,
list, topLevel);
TIMING_END_AND_PRINT (gcd_recursion,
"time for recursive call: ");
DEBOUTLN (cerr, "G_random_element= " << G_random_element);
}
if (!G_random_element.inCoeffDomain())
d0= totaldegree (G_random_element, Variable(2),
Variable (G_random_element.level()));
else
d0= 0;
if (d0 == 0)
{
if (inextension)
prune1 (alpha);
return N(gcdcAcB);
}
if (d0 > d)
{
if (!find (l, random_element))
l.append (random_element);
continue;
}
G_random_element=
(gcdlcAlcB(random_element, x)/uni_lcoeff (G_random_element))
* G_random_element;
skeleton= G_random_element;
if (!G_random_element.inCoeffDomain())
d0= totaldegree (G_random_element, Variable(2),
Variable (G_random_element.level()));
else
d0= 0;
if (d0 < d)
{
m= gcdlcAlcB;
newtonPoly= 1;
G_m= 0;
d= d0;
}
H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
if (uni_lcoeff (H) == gcdlcAlcB)
{
cH= uni_content (H);
ppH= H/cH;
if (inextension)
{
CFList u, v;
//maybe it's better to test if ppH is an element of F(\alpha) before
//mapping down
if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
{
DEBOUTLN (cerr, "ppH before mapDown= " << ppH);
ppH= mapDown (ppH, prim_elem_alpha, im_prim_elem_alpha, alpha, u, v);
ppH /= Lc(ppH);
DEBOUTLN (cerr, "ppH after mapDown= " << ppH);
prune1 (alpha);
return N(gcdcAcB*ppH);
}
}
else if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
return N(gcdcAcB*ppH);
}
G_m= H;
newtonPoly= newtonPoly*(x - random_element);
m= m*(x - random_element);
if (!find (l, random_element))
l.append (random_element);
if (getCharacteristic () > 3 && size (skeleton) < 100)
{
CFArray Monoms;
CFArray *coeffMonoms;
do //second do
{
random_element= randomElement (m, V_buf, l, fail);
if (random_element == 0 && !fail)
{
if (!find (l, random_element))
l.append (random_element);
continue;
}
if (fail)
{
source= CFList();
dest= CFList();
Variable V_buf3= V_buf;
V_buf= chooseExtension (V_buf);
bool prim_fail= false;
Variable V_buf2;
prim_elem= primitiveElement (V_buf4, V_buf2, prim_fail);
if (V_buf4 == alpha)
prim_elem_alpha= prim_elem;
if (V_buf3 != V_buf4)
{
m= mapDown (m, prim_elem, im_prim_elem, V_buf4, source, dest);
G_m= mapDown (m, prim_elem, im_prim_elem, V_buf4, source, dest);
newtonPoly= mapDown (newtonPoly, prim_elem, im_prim_elem, V_buf4,
source, dest);
ppA= mapDown (ppA, prim_elem, im_prim_elem, V_buf4, source, dest);
ppB= mapDown (ppB, prim_elem, im_prim_elem, V_buf4, source, dest);
gcdlcAlcB= mapDown (gcdlcAlcB, prim_elem, im_prim_elem, V_buf4,
source, dest);
for (CFListIterator i= l; i.hasItem(); i++)
i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, V_buf4,
source, dest);
}
ASSERT (!prim_fail, "failure in integer factorizer");
if (prim_fail)
; //ERROR
else
im_prim_elem= mapPrimElem (prim_elem, V_buf4, V_buf);
if (V_buf4 == alpha)
im_prim_elem_alpha= im_prim_elem;
else
im_prim_elem_alpha= mapUp (im_prim_elem_alpha, V_buf4, V_buf,
prim_elem, im_prim_elem, source, dest);
DEBOUTLN (cerr, "getMipo (V_buf4)= " << getMipo (V_buf4));
DEBOUTLN (cerr, "getMipo (V_buf2)= " << getMipo (V_buf2));
inextension= true;
for (CFListIterator i= l; i.hasItem(); i++)
i.getItem()= mapUp (i.getItem(), V_buf4, V_buf, prim_elem,
im_prim_elem, source, dest);
m= mapUp (m, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
G_m= mapUp (G_m, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
newtonPoly= mapUp (newtonPoly, V_buf4, V_buf, prim_elem, im_prim_elem,
source, dest);
ppA= mapUp (ppA, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
ppB= mapUp (ppB, V_buf4, V_buf, prim_elem, im_prim_elem, source, dest);
gcdlcAlcB= mapUp (gcdlcAlcB, V_buf4, V_buf, prim_elem, im_prim_elem,
source, dest);
fail= false;
random_element= randomElement (m, V_buf, l, fail);
DEBOUTLN (cerr, "fail= " << fail);
CFList list;
TIMING_START (gcd_recursion);
V_buf4= V_buf;
//sparseInterpolation
bool sparseFail= false;
if (LC (skeleton).inCoeffDomain())
G_random_element=
monicSparseInterpol (ppA (random_element, x), ppB(random_element,x),
skeleton,V_buf, sparseFail, coeffMonoms,Monoms);
else
G_random_element=
nonMonicSparseInterpol (ppA(random_element,x),ppB(random_element,x),
skeleton, V_buf, sparseFail, coeffMonoms,
Monoms);
if (sparseFail)
break;
TIMING_END_AND_PRINT (gcd_recursion,
"time for recursive call: ");
DEBOUTLN (cerr, "G_random_element= " << G_random_element);
}
else
{
CFList list;
TIMING_START (gcd_recursion);
bool sparseFail= false;
if (LC (skeleton).inCoeffDomain())
G_random_element=
monicSparseInterpol (ppA (random_element, x),ppB(random_element, x),
skeleton,V_buf, sparseFail,coeffMonoms, Monoms);
else
G_random_element=
nonMonicSparseInterpol (ppA(random_element,x),ppB(random_element,x),
skeleton, V_buf, sparseFail, coeffMonoms,
Monoms);
if (sparseFail)
break;
TIMING_END_AND_PRINT (gcd_recursion,
"time for recursive call: ");
DEBOUTLN (cerr, "G_random_element= " << G_random_element);
}
if (!G_random_element.inCoeffDomain())
d0= totaldegree (G_random_element, Variable(2),
Variable (G_random_element.level()));
else
d0= 0;
if (d0 == 0)
{
if (inextension)
prune1 (alpha);
return N(gcdcAcB);
}
if (d0 > d)
{
if (!find (l, random_element))
l.append (random_element);
continue;
}
G_random_element=
(gcdlcAlcB(random_element, x)/uni_lcoeff (G_random_element))
* G_random_element;
if (!G_random_element.inCoeffDomain())
d0= totaldegree (G_random_element, Variable(2),
Variable (G_random_element.level()));
else
d0= 0;
if (d0 < d)
{
m= gcdlcAlcB;
newtonPoly= 1;
G_m= 0;
d= d0;
}
TIMING_START (newton_interpolation);
H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
TIMING_END_AND_PRINT (newton_interpolation,
"time for newton interpolation: ");
//termination test
if (uni_lcoeff (H) == gcdlcAlcB)
{
cH= uni_content (H);
ppH= H/cH;
if (inextension)
{
CFList u, v;
//maybe it's better to test if ppH is an element of F(\alpha) before
//mapping down
if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
{
DEBOUTLN (cerr, "ppH before mapDown= " << ppH);
ppH= mapDown (ppH, prim_elem_alpha, im_prim_elem_alpha, alpha, u, v);
ppH /= Lc(ppH);
DEBOUTLN (cerr, "ppH after mapDown= " << ppH);
prune1 (alpha);
return N(gcdcAcB*ppH);
}
}
else if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
{
return N(gcdcAcB*ppH);
}
}
G_m= H;
newtonPoly= newtonPoly*(x - random_element);
m= m*(x - random_element);
if (!find (l, random_element))
l.append (random_element);
} while (1);
}
} while (1);
}
CanonicalForm sparseGCDFp (const CanonicalForm& F, const CanonicalForm& G,
bool& topLevel, CFList& l)
{
CanonicalForm A= F;
CanonicalForm B= G;
if (F.isZero() && degree(G) > 0) return G/Lc(G);
else if (G.isZero() && degree (F) > 0) return F/Lc(F);
else if (F.isZero() && G.isZero()) return F.genOne();
if (F.inBaseDomain() || G.inBaseDomain()) return F.genOne();
if (F.isUnivariate() && fdivides(F, G)) return F/Lc(F);
if (G.isUnivariate() && fdivides(G, F)) return G/Lc(G);
if (F == G) return F/Lc(F);
CFMap M,N;
int best_level= myCompress (A, B, M, N, topLevel);
if (best_level == 0) return B.genOne();
A= M(A);
B= M(B);
Variable x= Variable (1);
//univariate case
if (A.isUnivariate() && B.isUnivariate())
return N (gcd (A, B));
CanonicalForm cA, cB; // content of A and B
CanonicalForm ppA, ppB; // primitive part of A and B
CanonicalForm gcdcAcB;
cA = uni_content (A);
cB = uni_content (B);
gcdcAcB= gcd (cA, cB);
ppA= A/cA;
ppB= B/cB;
CanonicalForm lcA, lcB; // leading coefficients of A and B
CanonicalForm gcdlcAlcB;
lcA= uni_lcoeff (ppA);
lcB= uni_lcoeff (ppB);
if (fdivides (lcA, lcB))
{
if (fdivides (A, B))
return F/Lc(F);
}
if (fdivides (lcB, lcA))
{
if (fdivides (B, A))
return G/Lc(G);
}
gcdlcAlcB= gcd (lcA, lcB);
int d= totaldegree (ppA, Variable (2), Variable (ppA.level()));
int d0;
if (d == 0)
return N(gcdcAcB);
d0= totaldegree (ppB, Variable (2), Variable (ppB.level()));
if (d0 < d)
d= d0;
if (d == 0)
return N(gcdcAcB);
CanonicalForm m, random_element, G_m, G_random_element, H, cH, ppH, skeleton;
CanonicalForm newtonPoly= 1;
m= gcdlcAlcB;
G_m= 0;
H= 0;
bool fail= false;
topLevel= false;
bool inextension= false;
Variable V_buf, alpha, cleanUp;
CanonicalForm prim_elem, im_prim_elem;
CFList source, dest;
do //first do
{
if (inextension)
random_element= randomElement (m, V_buf, l, fail);
else
random_element= FpRandomElement (m, l, fail);
if (random_element == 0 && !fail)
{
if (!find (l, random_element))
l.append (random_element);
continue;
}
if (!fail && !inextension)
{
CFList list;
TIMING_START (gcd_recursion);
G_random_element=
sparseGCDFp (ppA (random_element,x), ppB (random_element,x), topLevel,
list);
TIMING_END_AND_PRINT (gcd_recursion,
"time for recursive call: ");
DEBOUTLN (cerr, "G_random_element= " << G_random_element);
}
else if (!fail && inextension)
{
CFList list;
TIMING_START (gcd_recursion);
G_random_element=
sparseGCDFq (ppA (random_element, x), ppB (random_element, x), alpha,
list, topLevel);
TIMING_END_AND_PRINT (gcd_recursion,
"time for recursive call: ");
DEBOUTLN (cerr, "G_random_element= " << G_random_element);
}
else if (fail && !inextension)
{
source= CFList();
dest= CFList();
CFList list;
CanonicalForm mipo;
int deg= 2;
bool initialized= false;
do
{
mipo= randomIrredpoly (deg, x);
if (initialized)
setMipo (alpha, mipo);
else
alpha= rootOf (mipo);
inextension= true;
fail= false;
initialized= true;
random_element= randomElement (m, alpha, l, fail);
deg++;
} while (fail);
cleanUp= alpha;
V_buf= alpha;
list= CFList();
TIMING_START (gcd_recursion);
G_random_element=
sparseGCDFq (ppA (random_element, x), ppB (random_element, x), alpha,
list, topLevel);
TIMING_END_AND_PRINT (gcd_recursion,
"time for recursive call: ");
DEBOUTLN (cerr, "G_random_element= " << G_random_element);
}
else if (fail && inextension)
{
source= CFList();
dest= CFList();
Variable V_buf3= V_buf;
V_buf= chooseExtension (V_buf);
bool prim_fail= false;
Variable V_buf2;
CanonicalForm prim_elem, im_prim_elem;
prim_elem= primitiveElement (alpha, V_buf2, prim_fail);
if (V_buf3 != alpha)
{
m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
G_m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
newtonPoly= mapDown (newtonPoly, prim_elem, im_prim_elem, alpha, source,
dest);
ppA= mapDown (ppA, prim_elem, im_prim_elem, alpha, source, dest);
ppB= mapDown (ppB, prim_elem, im_prim_elem, alpha, source, dest);
gcdlcAlcB= mapDown (gcdlcAlcB, prim_elem, im_prim_elem, alpha, source,
dest);
for (CFListIterator i= l; i.hasItem(); i++)
i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, alpha,
source, dest);
}
ASSERT (!prim_fail, "failure in integer factorizer");
if (prim_fail)
; //ERROR
else
im_prim_elem= mapPrimElem (prim_elem, alpha, V_buf);
DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (alpha));
DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf2));
for (CFListIterator i= l; i.hasItem(); i++)
i.getItem()= mapUp (i.getItem(), alpha, V_buf, prim_elem,
im_prim_elem, source, dest);
m= mapUp (m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
G_m= mapUp (G_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
newtonPoly= mapUp (newtonPoly, alpha, V_buf, prim_elem, im_prim_elem,
source, dest);
ppA= mapUp (ppA, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
ppB= mapUp (ppB, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
gcdlcAlcB= mapUp (gcdlcAlcB, alpha, V_buf, prim_elem, im_prim_elem,
source, dest);
fail= false;
random_element= randomElement (m, V_buf, l, fail );
DEBOUTLN (cerr, "fail= " << fail);
CFList list;
TIMING_START (gcd_recursion);
G_random_element=
sparseGCDFq (ppA (random_element, x), ppB (random_element, x), V_buf,
list, topLevel);
TIMING_END_AND_PRINT (gcd_recursion,
"time for recursive call: ");
DEBOUTLN (cerr, "G_random_element= " << G_random_element);
alpha= V_buf;
}
if (!G_random_element.inCoeffDomain())
d0= totaldegree (G_random_element, Variable(2),
Variable (G_random_element.level()));
else
d0= 0;
if (d0 == 0)
{
if (inextension)
prune (cleanUp);
return N(gcdcAcB);
}
if (d0 > d)
{
if (!find (l, random_element))
l.append (random_element);
continue;
}
G_random_element=
(gcdlcAlcB(random_element, x)/uni_lcoeff (G_random_element))
* G_random_element;
skeleton= G_random_element;
if (!G_random_element.inCoeffDomain())
d0= totaldegree (G_random_element, Variable(2),
Variable (G_random_element.level()));
else
d0= 0;
if (d0 < d)
{
m= gcdlcAlcB;
newtonPoly= 1;
G_m= 0;
d= d0;
}
H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
if (uni_lcoeff (H) == gcdlcAlcB)
{
cH= uni_content (H);
ppH= H/cH;
ppH /= Lc (ppH);
DEBOUTLN (cerr, "ppH= " << ppH);
if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
{
if (inextension)
prune (cleanUp);
return N(gcdcAcB*ppH);
}
}
G_m= H;
newtonPoly= newtonPoly*(x - random_element);
m= m*(x - random_element);
if (!find (l, random_element))
l.append (random_element);
if ((getCharacteristic() > 3 && size (skeleton) < 200))
{
CFArray Monoms;
CFArray* coeffMonoms;
do //second do
{
if (inextension)
random_element= randomElement (m, alpha, l, fail);
else
random_element= FpRandomElement (m, l, fail);
if (random_element == 0 && !fail)
{
if (!find (l, random_element))
l.append (random_element);
continue;
}
bool sparseFail= false;
if (!fail && !inextension)
{
CFList list;
TIMING_START (gcd_recursion);
if (LC (skeleton).inCoeffDomain())
G_random_element=
monicSparseInterpol(ppA(random_element, x), ppB (random_element, x),
skeleton, x, sparseFail, coeffMonoms,
Monoms);
else
G_random_element=
nonMonicSparseInterpol(ppA(random_element,x), ppB(random_element,x),
skeleton, x, sparseFail,
coeffMonoms, Monoms);
TIMING_END_AND_PRINT (gcd_recursion,
"time for recursive call: ");
DEBOUTLN (cerr, "G_random_element= " << G_random_element);
}
else if (!fail && inextension)
{
CFList list;
TIMING_START (gcd_recursion);
if (LC (skeleton).inCoeffDomain())
G_random_element=
monicSparseInterpol(ppA (random_element,x), ppB (random_element, x),
skeleton, alpha, sparseFail, coeffMonoms,
Monoms);
else
G_random_element=
nonMonicSparseInterpol(ppA(random_element,x), ppB(random_element,x),
skeleton, alpha, sparseFail, coeffMonoms,
Monoms);
TIMING_END_AND_PRINT (gcd_recursion,
"time for recursive call: ");
DEBOUTLN (cerr, "G_random_element= " << G_random_element);
}
else if (fail && !inextension)
{
source= CFList();
dest= CFList();
CFList list;
CanonicalForm mipo;
int deg= 2;
bool initialized= false;
do
{
mipo= randomIrredpoly (deg, x);
if (initialized)
setMipo (alpha, mipo);
else
alpha= rootOf (mipo);
inextension= true;
fail= false;
initialized= true;
random_element= randomElement (m, alpha, l, fail);
deg++;
} while (fail);
cleanUp= alpha;
V_buf= alpha;
list= CFList();
TIMING_START (gcd_recursion);
if (LC (skeleton).inCoeffDomain())
G_random_element=
monicSparseInterpol (ppA (random_element,x), ppB (random_element,x),
skeleton, alpha, sparseFail, coeffMonoms,
Monoms);
else
G_random_element=
nonMonicSparseInterpol(ppA(random_element,x), ppB(random_element,x),
skeleton, alpha, sparseFail, coeffMonoms,
Monoms);
TIMING_END_AND_PRINT (gcd_recursion,
"time for recursive call: ");
DEBOUTLN (cerr, "G_random_element= " << G_random_element);
}
else if (fail && inextension)
{
source= CFList();
dest= CFList();
Variable V_buf3= V_buf;
V_buf= chooseExtension (V_buf);
bool prim_fail= false;
Variable V_buf2;
CanonicalForm prim_elem, im_prim_elem;
prim_elem= primitiveElement (alpha, V_buf2, prim_fail);
if (V_buf3 != alpha)
{
m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
G_m= mapDown (m, prim_elem, im_prim_elem, alpha, source, dest);
newtonPoly= mapDown (newtonPoly, prim_elem, im_prim_elem, alpha,
source, dest);
ppA= mapDown (ppA, prim_elem, im_prim_elem, alpha, source, dest);
ppB= mapDown (ppB, prim_elem, im_prim_elem, alpha, source, dest);
gcdlcAlcB= mapDown (gcdlcAlcB, prim_elem, im_prim_elem, alpha,
source, dest);
for (CFListIterator i= l; i.hasItem(); i++)
i.getItem()= mapDown (i.getItem(), prim_elem, im_prim_elem, alpha,
source, dest);
}
ASSERT (!prim_fail, "failure in integer factorizer");
if (prim_fail)
; //ERROR
else
im_prim_elem= mapPrimElem (prim_elem, alpha, V_buf);
DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (alpha));
DEBOUTLN (cerr, "getMipo (alpha)= " << getMipo (V_buf2));
for (CFListIterator i= l; i.hasItem(); i++)
i.getItem()= mapUp (i.getItem(), alpha, V_buf, prim_elem,
im_prim_elem, source, dest);
m= mapUp (m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
G_m= mapUp (G_m, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
newtonPoly= mapUp (newtonPoly, alpha, V_buf, prim_elem, im_prim_elem,
source, dest);
ppA= mapUp (ppA, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
ppB= mapUp (ppB, alpha, V_buf, prim_elem, im_prim_elem, source, dest);
gcdlcAlcB= mapUp (gcdlcAlcB, alpha, V_buf, prim_elem, im_prim_elem,
source, dest);
fail= false;
random_element= randomElement (m, V_buf, l, fail );
DEBOUTLN (cerr, "fail= " << fail);
CFList list;
TIMING_START (gcd_recursion);
if (LC (skeleton).inCoeffDomain())
G_random_element=
monicSparseInterpol (ppA (random_element, x), ppB (random_element, x),
skeleton, V_buf, sparseFail, coeffMonoms,
Monoms);
else
G_random_element=
nonMonicSparseInterpol (ppA(random_element,x), ppB(random_element,x),
skeleton, V_buf, sparseFail,
coeffMonoms, Monoms);
TIMING_END_AND_PRINT (gcd_recursion,
"time for recursive call: ");
DEBOUTLN (cerr, "G_random_element= " << G_random_element);
alpha= V_buf;
}
if (sparseFail)
break;
if (!G_random_element.inCoeffDomain())
d0= totaldegree (G_random_element, Variable(2),
Variable (G_random_element.level()));
else
d0= 0;
if (d0 == 0)
{
if (inextension)
prune (cleanUp);
return N(gcdcAcB);
}
if (d0 > d)
{
if (!find (l, random_element))
l.append (random_element);
continue;
}
G_random_element=
(gcdlcAlcB(random_element, x)/uni_lcoeff (G_random_element))
* G_random_element;
if (!G_random_element.inCoeffDomain())
d0= totaldegree (G_random_element, Variable(2),
Variable (G_random_element.level()));
else
d0= 0;
if (d0 < d)
{
m= gcdlcAlcB;
newtonPoly= 1;
G_m= 0;
d= d0;
}
TIMING_START (newton_interpolation);
H= newtonInterp (random_element, G_random_element, newtonPoly, G_m, x);
TIMING_END_AND_PRINT (newton_interpolation,
"time for newton interpolation: ");
//termination test
if (uni_lcoeff (H) == gcdlcAlcB)
{
cH= uni_content (H);
ppH= H/cH;
ppH /= Lc (ppH);
DEBOUTLN (cerr, "ppH= " << ppH);
if (fdivides (ppH, ppA) && fdivides (ppH, ppB))
{
if (inextension)
prune (cleanUp);
return N(gcdcAcB*ppH);
}
}
G_m= H;
newtonPoly= newtonPoly*(x - random_element);
m= m*(x - random_element);
if (!find (l, random_element))
l.append (random_element);
} while (1); //end of second do
}
else
{
if (inextension)
prune (cleanUp);
return N(gcdcAcB*modGCDFp (ppA, ppB));
}
} while (1); //end of first do
}
TIMING_DEFINE_PRINT(modZ_termination)
TIMING_DEFINE_PRINT(modZ_recursion)
/// modular gcd algorithm, see Keith, Czapor, Geddes "Algorithms for Computer
/// Algebra", Algorithm 7.1
CanonicalForm modGCDZ ( const CanonicalForm & FF, const CanonicalForm & GG )
{
CanonicalForm f, g, cl, q(0), Dp, newD, D, newq, newqh;
int p, i, dp_deg, d_deg=-1;
CanonicalForm cd ( bCommonDen( FF ));
f=cd*FF;
Variable x= Variable (1);
CanonicalForm cf, cg;
cf= icontent (f);
f /= cf;
//cd = bCommonDen( f ); f *=cd;
//f /=vcontent(f,Variable(1));
cd = bCommonDen( GG );
g=cd*GG;
cg= icontent (g);
g /= cg;
//cd = bCommonDen( g ); g *=cd;
//g /=vcontent(g,Variable(1));
CanonicalForm Dn, test= 0;
CanonicalForm lcf, lcg;
lcf= f.lc();
lcg= g.lc();
cl = gcd (f.lc(),g.lc());
CanonicalForm gcdcfcg= gcd (cf, cg);
CanonicalForm fp, gp;
CanonicalForm b= 1;
int minCommonDeg= 0;
for (i= tmax (f.level(), g.level()); i > 0; i--)
{
if (degree (f, i) <= 0 || degree (g, i) <= 0)
continue;
else
{
minCommonDeg= tmin (degree (g, i), degree (f, i));
break;
}
}
if (i == 0)
return gcdcfcg;
for (; i > 0; i--)
{
if (degree (f, i) <= 0 || degree (g, i) <= 0)
continue;
else
minCommonDeg= tmin (minCommonDeg, tmin (degree (g, i), degree (f, i)));
}
b= 2*tmin (maxNorm (f), maxNorm (g))*abs (cl)*
power (CanonicalForm (2), minCommonDeg);
bool equal= false;
i = cf_getNumBigPrimes() - 1;
CanonicalForm cof, cog, cofp, cogp, newCof, newCog, cofn, cogn, cDn;
int maxNumVars= tmax (getNumVars (f), getNumVars (g));
//Off (SW_RATIONAL);
while ( true )
{
p = cf_getBigPrime( i );
i--;
while ( i >= 0 && mod( cl*(lc(f)/cl)*(lc(g)/cl), p ) == 0 )
{
p = cf_getBigPrime( i );
i--;
}
//printf("try p=%d\n",p);
setCharacteristic( p );
fp= mapinto (f);
gp= mapinto (g);
TIMING_START (modZ_recursion)
#if defined(HAVE_NTL) || defined(HAVE_FLINT)
if (size (fp)/maxNumVars > 500 && size (gp)/maxNumVars > 500)
Dp = modGCDFp (fp, gp, cofp, cogp);
else
{
Dp= gcd_poly (fp, gp);
cofp= fp/Dp;
cogp= gp/Dp;
}
#else
Dp= gcd_poly (fp, gp);
cofp= fp/Dp;
cogp= gp/Dp;
#endif
TIMING_END_AND_PRINT (modZ_recursion,
"time for gcd mod p in modular gcd: ");
Dp /=Dp.lc();
Dp *= mapinto (cl);
cofp /= lc (cofp);
cofp *= mapinto (lcf);
cogp /= lc (cogp);
cogp *= mapinto (lcg);
setCharacteristic( 0 );
dp_deg=totaldegree(Dp);
if ( dp_deg == 0 )
{
//printf(" -> 1\n");
return CanonicalForm(gcdcfcg);
}
if ( q.isZero() )
{
D = mapinto( Dp );
cof= mapinto (cofp);
cog= mapinto (cogp);
d_deg=dp_deg;
q = p;
Dn= balance_p (D, p);
cofn= balance_p (cof, p);
cogn= balance_p (cog, p);
}
else
{
if ( dp_deg == d_deg )
{
chineseRemainder( D, q, mapinto( Dp ), p, newD, newq );
chineseRemainder( cof, q, mapinto (cofp), p, newCof, newq);
chineseRemainder( cog, q, mapinto (cogp), p, newCog, newq);
cof= newCof;
cog= newCog;
newqh= newq/2;
Dn= balance_p (newD, newq, newqh);
cofn= balance_p (newCof, newq, newqh);
cogn= balance_p (newCog, newq, newqh);
if (test != Dn) //balance_p (newD, newq))
test= balance_p (newD, newq);
else
equal= true;
q = newq;
D = newD;
}
else if ( dp_deg < d_deg )
{
//printf(" were all bad, try more\n");
// all previous p's are bad primes
q = p;
D = mapinto( Dp );
cof= mapinto (cof);
cog= mapinto (cog);
d_deg=dp_deg;
test= 0;
equal= false;
Dn= balance_p (D, p);
cofn= balance_p (cof, p);
cogn= balance_p (cog, p);
}
else
{
//printf(" was bad, try more\n");
}
//else dp_deg > d_deg: bad prime
}
if ( i >= 0 )
{
cDn= icontent (Dn);
Dn /= cDn;
cofn /= cl/cDn;
//cofn /= icontent (cofn);
cogn /= cl/cDn;
//cogn /= icontent (cogn);
TIMING_START (modZ_termination);
if ((terminationTest (f,g, cofn, cogn, Dn)) ||
((equal || q > b) && fdivides (Dn, f) && fdivides (Dn, g)))
{
TIMING_END_AND_PRINT (modZ_termination,
"time for successful termination in modular gcd: ");
//printf(" -> success\n");
return Dn*gcdcfcg;
}
TIMING_END_AND_PRINT (modZ_termination,
"time for unsuccessful termination in modular gcd: ");
equal= false;
//else: try more primes
}
else
{ // try other method
//printf("try other gcd\n");
Off(SW_USE_CHINREM_GCD);
D=gcd_poly( f, g );
On(SW_USE_CHINREM_GCD);
return D*gcdcfcg;
}
}
}
#endif
convertNmod_mat_t2FacCFMatrix
CFMatrix * convertNmod_mat_t2FacCFMatrix(const nmod_mat_t m)
conversion of a FLINT matrix over Z/p to a factory matrix
Definition FLINTconvert.cc:630
convertFacCFMatrix2nmod_mat_t
void convertFacCFMatrix2nmod_mat_t(nmod_mat_t M, const CFMatrix &m)
conversion of a factory matrix over Z/p to a nmod_mat_t
Definition FLINTconvert.cc:612
convertFacCFMatrix2Fq_nmod_mat_t
void convertFacCFMatrix2Fq_nmod_mat_t(fq_nmod_mat_t M, const fq_nmod_ctx_t fq_con, const CFMatrix &m)
conversion of a factory matrix over F_q to a fq_nmod_mat_t
Definition FLINTconvert.cc:646
convertnmod_poly_t2FacCF
CanonicalForm convertnmod_poly_t2FacCF(const nmod_poly_t poly, const Variable &x)
conversion of a FLINT poly over Z/p to CanonicalForm
Definition FLINTconvert.cc:208
convertFq_nmod_mat_t2FacCFMatrix
CFMatrix * convertFq_nmod_mat_t2FacCFMatrix(const fq_nmod_mat_t m, const fq_nmod_ctx_t &fq_con, const Variable &alpha)
conversion of a FLINT matrix over F_q to a factory matrix
Definition FLINTconvert.cc:661
FLINTconvert.h
This file defines functions for conversion to FLINT (www.flintlib.org) and back.
pow
Rational pow(const Rational &a, int e)
Definition GMPrat.cc:411
abs
Rational abs(const Rational &a)
Definition GMPrat.cc:436
convertNTLmat_zz_p2FacCFMatrix
CFMatrix * convertNTLmat_zz_p2FacCFMatrix(const mat_zz_p &m)
Definition NTLconvert.cc:1183
convertNTLmat_zz_pE2FacCFMatrix
CFMatrix * convertNTLmat_zz_pE2FacCFMatrix(const mat_zz_pE &m, const Variable &alpha)
Definition NTLconvert.cc:1212
convertNTLzzpX2CF
CanonicalForm convertNTLzzpX2CF(const zz_pX &poly, const Variable &x)
Definition NTLconvert.cc:255
convertFacCFMatrix2NTLmat_zz_pE
mat_zz_pE * convertFacCFMatrix2NTLmat_zz_pE(const CFMatrix &m)
Definition NTLconvert.cc:1196
convertFacCF2NTLzzpX
zz_pX convertFacCF2NTLzzpX(const CanonicalForm &f)
Definition NTLconvert.cc:105
convertFacCFMatrix2NTLmat_zz_p
mat_zz_p * convertFacCFMatrix2NTLmat_zz_p(const CFMatrix &m)
Definition NTLconvert.cc:1167
fac_NTL_char
VAR long fac_NTL_char
Definition NTLconvert.cc:46
NTLconvert.h
Conversion to and from NTL.
On
void On(int sw)
switches
Definition canonicalform.cc:1957
Off
void Off(int sw)
switches
Definition canonicalform.cc:1964
power
CanonicalForm power(const CanonicalForm &f, int n)
exponentiation
Definition canonicalform.cc:1896
canonicalform.h
Header for factory's main class CanonicalForm.
content
CanonicalForm FACTORY_PUBLIC content(const CanonicalForm &)
CanonicalForm content ( const CanonicalForm & f ).
Definition cf_gcd.cc:603
mapinto
CanonicalForm mapinto(const CanonicalForm &f)
Definition canonicalform.h:355
size
int size(const CanonicalForm &f, const Variable &v)
int size ( const CanonicalForm & f, const Variable & v )
Definition cf_ops.cc:600
getNumVars
int getNumVars(const CanonicalForm &f)
int getNumVars ( const CanonicalForm & f )
Definition cf_ops.cc:314
lc
CanonicalForm lc(const CanonicalForm &f)
Definition canonicalform.h:304
icontent
CanonicalForm FACTORY_PUBLIC icontent(const CanonicalForm &f)
CanonicalForm icontent ( const CanonicalForm & f ).
Definition cf_gcd.cc:74
degrees
int * degrees(const CanonicalForm &f, int *degs=0)
int * degrees ( const CanonicalForm & f, int * degs )
Definition cf_ops.cc:493
mod
CF_NO_INLINE FACTORY_PUBLIC CanonicalForm mod(const CanonicalForm &, const CanonicalForm &)
degree
int degree(const CanonicalForm &f)
Definition canonicalform.h:316
getGFDegree
int getGFDegree()
Definition cf_char.cc:75
CFArray
Array< CanonicalForm > CFArray
Definition canonicalform.h:397
setCharacteristic
void FACTORY_PUBLIC setCharacteristic(int c)
Definition cf_char.cc:28
CFMatrix
Matrix< CanonicalForm > CFMatrix
Definition canonicalform.h:398
swapvar
CanonicalForm FACTORY_PUBLIC swapvar(const CanonicalForm &, const Variable &, const Variable &)
swapvar() - swap variables x1 and x2 in f.
Definition cf_ops.cc:168
CFListIterator
ListIterator< CanonicalForm > CFListIterator
Definition canonicalform.h:396
totaldegree
int totaldegree(const CanonicalForm &f)
int totaldegree ( const CanonicalForm & f )
Definition cf_ops.cc:523
Lc
CanonicalForm Lc(const CanonicalForm &f)
Definition canonicalform.h:307
gcd_poly
CanonicalForm FACTORY_PUBLIC gcd_poly(const CanonicalForm &f, const CanonicalForm &g)
CanonicalForm gcd_poly ( const CanonicalForm & f, const CanonicalForm & g ).
Definition cf_gcd.cc:492
CFList
List< CanonicalForm > CFList
Definition canonicalform.h:395
LC
CanonicalForm LC(const CanonicalForm &f)
Definition canonicalform.h:310
getCharacteristic
int FACTORY_PUBLIC getCharacteristic()
Definition cf_char.cc:70
N
const CanonicalForm CFMap CFMap & N
Definition cfEzgcd.cc:56
degsf
int * degsf
Definition cfEzgcd.cc:59
f_zero
int f_zero
Definition cfEzgcd.cc:69
l
int l
Definition cfEzgcd.cc:100
m
int m
Definition cfEzgcd.cc:128
both_non_zero
const CanonicalForm CFMap CFMap int & both_non_zero
Definition cfEzgcd.cc:57
i
int i
Definition cfEzgcd.cc:132
g_zero
int g_zero
Definition cfEzgcd.cc:70
k
int k
Definition cfEzgcd.cc:99
degsg
int * degsg
Definition cfEzgcd.cc:60
topLevel
const CanonicalForm CFMap CFMap bool topLevel
Definition cfGcdAlgExt.cc:58
both_zero
int both_zero
Definition cfGcdAlgExt.cc:72
cfGcdUtil.h
coprimality check and change of representation mod n
sparseGCDFq
CanonicalForm sparseGCDFq(const CanonicalForm &F, const CanonicalForm &G, const Variable &alpha, CFList &l, bool &topLevel)
Definition cfModGcd.cc:3138
readOffSolution
CFArray readOffSolution(const CFMatrix &M, const long rk)
M in row echolon form, rk rank of M.
Definition cfModGcd.cc:1689
modGCDFq
CanonicalForm modGCDFq(const CanonicalForm &F, const CanonicalForm &G, CanonicalForm &coF, CanonicalForm &coG, Variable &alpha, CFList &l, bool &topLevel)
GCD of F and G over , l and topLevel are only used internally, output is monic based on Alg....
Definition cfModGcd.cc:479
GFRandomElement
static CanonicalForm GFRandomElement(const CanonicalForm &F, CFList &list, bool &fail)
compute a random element a of GF, s.t. F(a) , F is a univariate polynomial, returns fail if there ar...
Definition cfModGcd.cc:820
evaluateMonom
CFArray evaluateMonom(const CanonicalForm &F, const CFList &evalPoints)
Definition cfModGcd.cc:2000
coG
const CanonicalForm const CanonicalForm const CanonicalForm & coG
Definition cfModGcd.cc:68
getMonoms
CFArray getMonoms(const CanonicalForm &F)
extract monomials of F, parts in algebraic variable are considered coefficients
Definition cfModGcd.cc:1965
mult
void mult(CFList &L1, const CFList &L2)
multiply two lists componentwise
Definition cfModGcd.cc:2184
cand
const CanonicalForm const CanonicalForm const CanonicalForm const CanonicalForm & cand
Definition cfModGcd.cc:70
cofp
CanonicalForm cofp
Definition cfModGcd.cc:4137
gaussianElimFq
long gaussianElimFq(CFMatrix &M, CFArray &L, const Variable &alpha)
Definition cfModGcd.cc:1785
myCompress
int myCompress(const CanonicalForm &F, const CanonicalForm &G, CFMap &M, CFMap &N, bool topLevel)
compressing two polynomials F and G, M is used for compressing, N to reverse the compression
Definition cfModGcd.cc:92
x
Variable x
Definition cfModGcd.cc:4090
lcg
CanonicalForm lcg
Definition cfModGcd.cc:4105
chooseExtension
static Variable chooseExtension(const Variable &alpha)
Definition cfModGcd.cc:421
dp_deg
int dp_deg
Definition cfModGcd.cc:4086
newCog
CanonicalForm newCog
Definition cfModGcd.cc:4137
evaluate
CFArray evaluate(const CFArray &A, const CFList &evalPoints)
Definition cfModGcd.cc:2039
monicSparseInterpol
CanonicalForm monicSparseInterpol(const CanonicalForm &F, const CanonicalForm &G, const CanonicalForm &skeleton, const Variable &alpha, bool &fail, CFArray *&coeffMonoms, CFArray &Monoms)
Definition cfModGcd.cc:2207
gaussianElimFp
long gaussianElimFp(CFMatrix &M, CFArray &L)
Definition cfModGcd.cc:1740
cogn
CanonicalForm cogn
Definition cfModGcd.cc:4137
p
int p
Definition cfModGcd.cc:4086
uni_content
static CanonicalForm uni_content(const CanonicalForm &F)
compute the content of F, where F is considered as an element of
Definition cfModGcd.cc:282
solveGeneralVandermonde
CFArray solveGeneralVandermonde(const CFArray &M, const CFArray &A)
Definition cfModGcd.cc:1633
cl
cl
Definition cfModGcd.cc:4108
extractContents
CanonicalForm extractContents(const CanonicalForm &F, const CanonicalForm &G, CanonicalForm &contentF, CanonicalForm &contentG, CanonicalForm &ppF, CanonicalForm &ppG, const int d)
Definition cfModGcd.cc:312
lcf
CanonicalForm lcf
Definition cfModGcd.cc:4105
uni_lcoeff
static CanonicalForm uni_lcoeff(const CanonicalForm &F)
compute the leading coefficient of F, where F is considered as an element of , order on is dp.
Definition cfModGcd.cc:340
g
g
Definition cfModGcd.cc:4098
maxNumVars
int maxNumVars
Definition cfModGcd.cc:4138
fp
CanonicalForm fp
Definition cfModGcd.cc:4110
solveVandermonde
CFArray solveVandermonde(const CFArray &M, const CFArray &A)
Definition cfModGcd.cc:1580
cogp
CanonicalForm cogp
Definition cfModGcd.cc:4137
coF
const CanonicalForm const CanonicalForm & coF
Definition cfModGcd.cc:68
cofn
CanonicalForm cofn
Definition cfModGcd.cc:4137
cof
CanonicalForm cof
Definition cfModGcd.cc:4137
GG
const CanonicalForm & GG
Definition cfModGcd.cc:4084
cd
CanonicalForm cd(bCommonDen(FF))
Definition cfModGcd.cc:4097
cog
CanonicalForm cog
Definition cfModGcd.cc:4137
nonMonicSparseInterpol
CanonicalForm nonMonicSparseInterpol(const CanonicalForm &F, const CanonicalForm &G, const CanonicalForm &skeleton, const Variable &alpha, bool &fail, CFArray *&coeffMonoms, CFArray &Monoms)
Definition cfModGcd.cc:2491
minCommonDeg
int minCommonDeg
Definition cfModGcd.cc:4112
modGCDFp
CanonicalForm modGCDFp(const CanonicalForm &F, const CanonicalForm &G, CanonicalForm &coF, CanonicalForm &coG, bool &topLevel, CFList &l)
Definition cfModGcd.cc:1224
solveSystemFp
CFArray solveSystemFp(const CFMatrix &M, const CFArray &L)
Definition cfModGcd.cc:1844
sparseGCDFp
CanonicalForm sparseGCDFp(const CanonicalForm &F, const CanonicalForm &G, bool &topLevel, CFList &l)
Definition cfModGcd.cc:3570
randomElement
static CanonicalForm randomElement(const CanonicalForm &F, const Variable &alpha, CFList &list, bool &fail)
compute a random element a of , s.t. F(a) , F is a univariate polynomial, returns fail if there are...
Definition cfModGcd.cc:380
gcdcfcg
CanonicalForm gcdcfcg
Definition cfModGcd.cc:4109
cf
CanonicalForm cf
Definition cfModGcd.cc:4091
Dn
CanonicalForm Dn
Definition cfModGcd.cc:4104
newCof
CanonicalForm newCof
Definition cfModGcd.cc:4137
FpRandomElement
static CanonicalForm FpRandomElement(const CanonicalForm &F, CFList &list, bool &fail)
Definition cfModGcd.cc:1172
gp
CanonicalForm gp
Definition cfModGcd.cc:4110
modGCDGF
CanonicalForm modGCDGF(const CanonicalForm &F, const CanonicalForm &G, CanonicalForm &coF, CanonicalForm &coG, CFList &l, bool &topLevel)
GCD of F and G over GF, based on Alg. 7.2. as described in "Algorithms forComputer Algebra" by Geddes...
Definition cfModGcd.cc:873
newtonInterp
static CanonicalForm newtonInterp(const CanonicalForm &alpha, const CanonicalForm &u, const CanonicalForm &newtonPoly, const CanonicalForm &oldInterPoly, const Variable &x)
Newton interpolation - Incremental algorithm. Given a list of values alpha_i and a list of polynomial...
Definition cfModGcd.cc:365
equal
bool equal
Definition cfModGcd.cc:4134
solveSystemFq
CFArray solveSystemFq(const CFMatrix &M, const CFArray &L, const Variable &alpha)
Definition cfModGcd.cc:1896
test
CanonicalForm test
Definition cfModGcd.cc:4104
b
CanonicalForm b
Definition cfModGcd.cc:4111
cg
CanonicalForm cg
Definition cfModGcd.cc:4091
cDn
CanonicalForm cDn
Definition cfModGcd.cc:4137
d_deg
int d_deg
Definition cfModGcd.cc:4086
evaluationPoints
CFList evaluationPoints(const CanonicalForm &F, const CanonicalForm &G, CanonicalForm &Feval, CanonicalForm &Geval, const CanonicalForm &LCF, const bool &GF, const Variable &alpha, bool &fail, CFList &list)
Definition cfModGcd.cc:2056
cfModGcd.h
modular and sparse modular GCD algorithms over finite fields and Z.
terminationTest
bool terminationTest(const CanonicalForm &F, const CanonicalForm &G, const CanonicalForm &coF, const CanonicalForm &coG, const CanonicalForm &cand)
modGCDZ
CanonicalForm modGCDZ(const CanonicalForm &FF, const CanonicalForm &GG)
modular GCD over Z
evalPoint
static void evalPoint(const CanonicalForm &F, const CanonicalForm &G, CanonicalForm &FEval, CanonicalForm &GEval, CFGenerator &evalPoint)
Definition cfModResultant.cc:311
cfNewtonPolygon.h
This file provides functions to compute the Newton polygon of a bivariate polynomial.
bCommonDen
CanonicalForm bCommonDen(const CanonicalForm &f)
void psqr ( const CanonicalForm & f, const CanonicalForm & g, CanonicalForm & q, CanonicalForm & r,...
Definition cf_algorithm.cc:293
maxNorm
CanonicalForm maxNorm(const CanonicalForm &f)
bool fdivides ( const CanonicalForm & f, const CanonicalForm & g )
Definition cf_algorithm.cc:536
cf_algorithm.h
declarations of higher level algorithms.
fdivides
bool fdivides(const CanonicalForm &f, const CanonicalForm &g)
chineseRemainder
void FACTORY_PUBLIC chineseRemainder(const CanonicalForm &x1, const CanonicalForm &q1, const CanonicalForm &x2, const CanonicalForm &q2, CanonicalForm &xnew, CanonicalForm &qnew)
void chineseRemainder ( const CanonicalForm & x1, const CanonicalForm & q1, const CanonicalForm & x2,...
Definition cf_chinese.cc:57
cf_assert.h
assertions for Factory
ASSERT
#define ASSERT(expression, message)
Definition cf_assert.h:99
SW_USE_CHINREM_GCD
static const int SW_USE_CHINREM_GCD
set to 1 to use modular gcd over Z
Definition cf_defs.h:41
DELETE_ARRAY
#define DELETE_ARRAY(P)
Definition cf_defs.h:65
NEW_ARRAY
#define NEW_ARRAY(T, N)
Definition cf_defs.h:64
balance_p
static CanonicalForm balance_p(const CanonicalForm &f, const CanonicalForm &q, const CanonicalForm &qh)
Definition cf_gcd.cc:149
randomIrredpoly
CanonicalForm randomIrredpoly(int i, const Variable &x)
computes a random monic irreducible univariate polynomial in x over Fp of degree i via NTL/FLINT
Definition cf_irred.cc:26
cf_irred.h
generate random irreducible univariate polynomials
cf_iter.h
Iterators for CanonicalForm's.
bound
static CanonicalForm bound(const CFMatrix &M)
Definition cf_linsys.cc:460
cf_map.h
map polynomials
mapPrimElem
CanonicalForm mapPrimElem(const CanonicalForm &primElem, const Variable &alpha, const Variable &beta)
compute the image of a primitive element of in . We assume .
Definition cf_map_ext.cc:451
GFMapDown
CanonicalForm GFMapDown(const CanonicalForm &F, int k)
maps a polynomial over to a polynomial over , d needs to be a multiple of k
Definition cf_map_ext.cc:277
primitiveElement
CanonicalForm primitiveElement(const Variable &alpha, Variable &beta, bool &fail)
determine a primitive element of , is a primitive element of a field which is isomorphic to
Definition cf_map_ext.cc:343
mapDown
static CanonicalForm mapDown(const CanonicalForm &F, const Variable &alpha, const CanonicalForm &G, CFList &source, CFList &dest)
the CanonicalForm G is the output of map_up, returns F considered as an element over ,...
Definition cf_map_ext.cc:124
mapUp
static CanonicalForm mapUp(const Variable &alpha, const Variable &beta)
and is a primitive element, returns the image of
Definition cf_map_ext.cc:71
GFMapUp
CanonicalForm GFMapUp(const CanonicalForm &F, int k)
maps a polynomial over to a polynomial over , d needs to be a multiple of k
Definition cf_map_ext.cc:241
cf_map_ext.h
This file implements functions to map between extensions of finite fields.
cf_getBigPrime
int cf_getBigPrime(int i)
Definition cf_primes.cc:39
cf_getNumBigPrimes
int cf_getNumBigPrimes()
Definition cf_primes.cc:45
cf_primes.h
access to prime tables
FLINTrandom
GLOBAL_VAR flint_rand_t FLINTrandom
Definition cf_random.cc:23
cf_random.h
generate random integers, random elements of finite fields
cf_reval.h
generate random evaluation points
factoryError
VAR void(* factoryError)(const char *s)
Definition cf_util.cc:80
ipower
int ipower(int b, int m)
int ipower ( int b, int m )
Definition cf_util.cc:27
cf_util.h
f
FILE * f
Definition checklibs.c:9
AlgExtRandomF
generate random elements in F_p(alpha)
Definition cf_random.h:70
AlgExtRandomF::generate
CanonicalForm generate() const
Definition cf_random.cc:154
Array::size
int size() const
Definition ftmpl_array.cc:92
CFIterator
class to iterate through CanonicalForm's
Definition cf_iter.h:44
CFMap
class CFMap
Definition cf_map.h:85
CanonicalForm
factory's main class
Definition canonicalform.h:86
CanonicalForm::genOne
CanonicalForm genOne() const
Definition canonicalform.cc:1881
CanonicalForm::isZero
CF_NO_INLINE bool isZero() const
CanonicalForm::mvar
Variable mvar() const
mvar() returns the main variable of CO or Variable() if CO is in a base domain.
Definition canonicalform.cc:593
CanonicalForm::isOne
CF_NO_INLINE bool isOne() const
CanonicalForm::inCoeffDomain
bool inCoeffDomain() const
Definition canonicalform.cc:122
CanonicalForm::level
int level() const
level() returns the level of CO.
Definition canonicalform.cc:576
CanonicalForm::inBaseDomain
bool inBaseDomain() const
Definition canonicalform.cc:104
CanonicalForm::lc
CanonicalForm lc() const
CanonicalForm CanonicalForm::lc (), Lc (), LC (), LC ( v ) const.
Definition canonicalform.cc:337
CanonicalForm::isUnivariate
bool isUnivariate() const
Definition canonicalform.cc:155
FFRandom
generate random elements in F_p
Definition cf_random.h:44
FFRandom::generate
CanonicalForm generate() const
Definition cf_random.cc:65
GFRandom
generate random elements in GF
Definition cf_random.h:32
GFRandom::generate
CanonicalForm generate() const
Definition cf_random.cc:75
List::length
int length() const
Definition ftmpl_list.cc:273
List::append
void append(const T &)
Definition ftmpl_list.cc:256
Matrix::rows
int rows() const
Definition ftmpl_matrix.h:43
Matrix::columns
int columns() const
Definition ftmpl_matrix.h:44
Variable
factory's class for variables
Definition factory.h:127
Variable::level
int level() const
Definition factory.h:143
debug.h
functions to print debug output
DEBOUTLN
#define DEBOUTLN(stream, objects)
Definition debug.h:49
iter
CFFListIterator iter
Definition facAbsBiFact.cc:54
alpha
Variable alpha
Definition facAbsBiFact.cc:52
result
return result
Definition facAbsBiFact.cc:76
Feval
CanonicalForm Feval
Definition facAbsFact.cc:60
H
CanonicalForm H
Definition facAbsFact.cc:60
mipo
CanonicalForm mipo
Definition facAlgExt.cc:57
B
b *CanonicalForm B
Definition facBivar.cc:52
v
const Variable & v
< [in] a sqrfree bivariate poly
Definition facBivar.h:39
LCF
CanonicalForm LCF
Definition facFactorize.cc:52
eval
CFList & eval
Definition facFactorize.cc:47
Aeval
CFList *& Aeval
<[in] poly
Definition facFactorize.h:31
tmp1
CFList tmp1
Definition facFqBivar.cc:75
tmp2
CFList tmp2
Definition facFqBivar.cc:75
evalPoints
CFList evalPoints(const CanonicalForm &F, CFList &eval, const Variable &alpha, CFList &list, const bool &GF, bool &fail)
evaluation point search for multivariate factorization, looks for a (F.level() - 1)-tuple such that t...
Definition facFqFactorize.cc:751
j
int j
Definition facHensel.cc:110
fq_nmod_ctx_clear
fq_nmod_ctx_clear(fq_con)
nmod_poly_init
nmod_poly_init(FLINTmipo, getCharacteristic())
fq_nmod_ctx_init_modulus
fq_nmod_ctx_init_modulus(fq_con, FLINTmipo, "Z")
convertFacCF2nmod_poly_t
convertFacCF2nmod_poly_t(FLINTmipo, M)
nmod_poly_clear
nmod_poly_clear(FLINTmipo)
facHensel.h
This file defines functions for Hensel lifting.
rootOf
Variable FACTORY_PUBLIC rootOf(const CanonicalForm &, char name='@')
returns a symbolic root of polynomial with name name Use it to define algebraic variables
Definition variable.cc:162
prune1
void prune1(const Variable &alpha)
Definition variable.cc:291
getMipo
CanonicalForm getMipo(const Variable &alpha, const Variable &x)
Definition variable.cc:207
prune
void FACTORY_PUBLIC prune(Variable &alpha)
Definition variable.cc:261
setMipo
void setMipo(const Variable &alpha, const CanonicalForm &mipo)
Definition variable.cc:219
const
#define const
Definition fegetopt.c:39
ftmpl_functions.h
some useful template functions.
tmax
template CanonicalForm tmax(const CanonicalForm &, const CanonicalForm &)
tmin
template CanonicalForm tmin(const CanonicalForm &, const CanonicalForm &)
find
template bool find(const List< CanonicalForm > &, const CanonicalForm &)
D
#define D(A)
Definition gentable.cc:128
gf_name
VAR char gf_name
Definition gfops.cc:52
G
STATIC_VAR TreeM * G
Definition janet.cc:31
log
gmp_float log(const gmp_float &a)
Definition mpr_complex.cc:335
NULL
#define NULL
Definition omList.c:12
buf
int status int void * buf
Definition si_signals.h:69
A
#define A
Definition sirandom.c:24
M
#define M
Definition sirandom.c:25
timing.h
TIMING_DEFINE_PRINT
#define TIMING_DEFINE_PRINT(t)
Definition timing.h:95
TIMING_START
#define TIMING_START(t)
Definition timing.h:92
TIMING_END_AND_PRINT
#define TIMING_END_AND_PRINT(t, msg)
Definition timing.h:94
gcd
int gcd(int a, int b)
Definition walkSupport.cc:836