launched benchs

parent cf4502f4
......@@ -45,28 +45,26 @@ BuildMonomialVector(size_t deg_curve, size_t deg_max) {
}
mat_ZZ_p
BuildInterpolationMatrixInDegree(const EffectiveDivisor& D, size_t deg) {
vector<ZZ_pX> powers_of_g;
vector<ZZ_pX> powers_of_x;
powers_of_g.push_back(ZZ_pX(1));
powers_of_x.push_back(ZZ_pX(1));
ZZ_pXModulus modf(D.get_f());
BuildInterpolationMatrixInDegree(const EffectiveDivisor& D, const EffectiveDivisor& E, size_t deg) {
vector<ZZ_pX> powers_of_gD;
vector<ZZ_pX> powers_of_xD;
powers_of_gD.push_back(ZZ_pX(1));
powers_of_xD.push_back(ZZ_pX(1));
ZZ_pXModulus modfD(D.get_f());
ZZ_pX x_pol(1);
x_pol <<= 1;
for (size_t i = 0; i <= min(deg, D.curve().degree()-1); ++i) {
ZZ_pX tmp;
MulMod(tmp, powers_of_g.back(), D.get_g(), modf);
powers_of_g.push_back(tmp);
MulMod(tmp, powers_of_gD.back(), D.get_g(), modfD);
powers_of_gD.push_back(tmp);
}
for (size_t i = 0; i <= deg; ++i) {
ZZ_pX tmp;
MulMod(tmp, powers_of_x.back(), x_pol % modf, modf);
powers_of_x.push_back(tmp);
MulMod(tmp, powers_of_xD.back(), x_pol % modfD, modfD);
powers_of_xD.push_back(tmp);
}
auto vec_mons = BuildMonomialVector(D.curve().degree(), deg);
......@@ -76,23 +74,50 @@ BuildInterpolationMatrixInDegree(const EffectiveDivisor& D, size_t deg) {
for (const auto& p : vec_mons) {
tmp_row.SetLength(0);
MulMod(tmp_pol, powers_of_x[p.first], powers_of_g[p.second], modf);
MulMod(tmp_pol, powers_of_xD[p.first], powers_of_gD[p.second], modfD);
for (size_t i = 0; i < D.degree(); ++i)
tmp_row.append(coeff(tmp_pol, i));
vec_vec_res.append(tmp_row);
}
if (E.degree() > 0) { // singular case
vector<ZZ_pX> powers_of_gE;
vector<ZZ_pX> powers_of_xE;
powers_of_gE.push_back(ZZ_pX(1));
powers_of_xE.push_back(ZZ_pX(1));
ZZ_pXModulus modfE(E.get_f());
for (size_t i = 0; i <= min(deg, D.curve().degree()-1); ++i) {
ZZ_pX tmp;
MulMod(tmp, powers_of_gE.back(), E.get_g(), modfE);
powers_of_gE.push_back(tmp);
}
for (size_t i = 0; i <= deg; ++i) {
ZZ_pX tmp;
MulMod(tmp, powers_of_xE.back(), x_pol % modfE, modfE);
powers_of_xE.push_back(tmp);
}
for (size_t i = 0; i < vec_mons.size(); ++i) {
tmp_row.SetLength(0);
MulMod(tmp_pol, powers_of_xE[vec_mons[i].first], powers_of_gE[vec_mons[i].second], modfE);
for (size_t ll = 0; ll < E.degree(); ++ll)
vec_vec_res[i].append(coeff(tmp_pol, ll));
}
}
mat_ZZ_p res;
res.SetDims(vec_mons.size(), D.degree());
res.SetDims(vec_mons.size(), D.degree()+E.degree());
MakeMatrix(res, vec_vec_res);
return res;
}
BivPol
Interpolate(const EffectiveDivisor& D) {
size_t deg_res = GetInterpolationDegree(D.curve().degree(), D.degree());
mat_ZZ_p M = BuildInterpolationMatrixInDegree(D, deg_res);
Interpolate(const EffectiveDivisor& D, const EffectiveDivisor& E) {
size_t deg_res = GetInterpolationDegree(D.curve().degree(), D.degree() + E.degree());
mat_ZZ_p M = BuildInterpolationMatrixInDegree(D, E, deg_res);
mat_ZZ_p kerM;
kernel(kerM, M);
mat_ZZ_p randmat;
......@@ -151,7 +176,12 @@ PrincipalDivisor(const Curve& C, const BivPol& h, const EffectiveDivisor& E) {
ZZ_pX subres0 = interpolate(eval_pts, eval_subres0);
ZZ_pX subres1 = interpolate(eval_pts, eval_subres1);
return PositiveDifference(EffectiveDivisor(C, res, -MulMod(InvMod(subres1, res), subres0, res)), E);
div(res, res, E.get_f());
div(res, res, E.get_f());
rem(subres0, subres0, res);
rem(subres1, subres1, res);
return EffectiveDivisor(C, res, -MulMod(InvMod(subres1, res), subres0, res));
}
EffectiveDivisor
......@@ -162,26 +192,26 @@ PositiveDifference(const EffectiveDivisor& D1, const EffectiveDivisor& D2) {
return EffectiveDivisor(D1.curve(), f/LeadCoeff(f), D1.get_g() % f);
}
// This simple sum only works if D1 and D2 do not have any point with the same
// x-coordinate in their support
EffectiveDivisor
SimpleSum(const EffectiveDivisor& D1, const EffectiveDivisor& D2) {
assert(D1.curve() == D2.curve());
if (deg(D1.get_f()) == 0)
return D2;
if (deg(D2.get_f()) == 0)
return D1;
NTL::ZZ_pX gcd, a1, a2;
XGCD(gcd, a1, a2, D1.get_f(), D2.get_f());
assert(deg(gcd) == 0);
ZZ_pX new_f = D1.get_f()*D2.get_f();
ZZ_pX new_g = (D1.get_g()*a2*D2.get_f() +
D2.get_g()*a1*D1.get_f()) % new_f;
assert(new_g % D1.get_f() == D1.get_g());
assert(new_g % D2.get_f() == D2.get_g());
return EffectiveDivisor(D1.curve(), new_f, new_g);
}
// // This simple sum only works if D1 and D2 do not have any point with the same
// // x-coordinate in their support
// EffectiveDivisor
// SimpleSum(const EffectiveDivisor& D1, const EffectiveDivisor& D2) {
// assert(D1.curve() == D2.curve());
// if (deg(D1.get_f()) == 0)
// return D2;
// if (deg(D2.get_f()) == 0)
// return D1;
// NTL::ZZ_pX gcd, a1, a2;
// XGCD(gcd, a1, a2, D1.get_f(), D2.get_f());
// assert(deg(gcd) == 0);
//
// ZZ_pX new_f = D1.get_f()*D2.get_f();
// ZZ_pX new_g = (D1.get_g()*a2*D2.get_f() +
// D2.get_g()*a1*D1.get_f()) % new_f;
// assert(new_g % D1.get_f() == D1.get_g());
// assert(new_g % D2.get_f() == D2.get_g());
// return EffectiveDivisor(D1.curve(), new_f, new_g);
// }
// This adding function works only if all points in the support are
// non-singular
......@@ -238,9 +268,9 @@ MultiplyByInt(const Divisor& D, unsigned int k) {
}
vector<BivPol>
EffectiveBasisRRinDegree(const EffectiveDivisor& D, size_t deg) {
EffectiveBasisRRinDegree(const EffectiveDivisor& D, const EffectiveDivisor& E, size_t deg) {
vector<BivPol> res;
mat_ZZ_p M = BuildInterpolationMatrixInDegree(D, deg);
mat_ZZ_p M = BuildInterpolationMatrixInDegree(D, E, deg);
mat_ZZ_p kerM;
kernel(kerM, M);
for (long k = 0; k < kerM.NumRows(); ++k) {
......@@ -266,10 +296,9 @@ RiemannRochBasis(const Divisor& D, const EffectiveDivisor& E) {
assert(D.curve().get_pdefpol()->mod_eval(D.get_pos().get_g(), D.get_pos().get_f()) == 0);
assert(D.curve().get_pdefpol()->mod_eval(D.get_neg().get_g(), D.get_neg().get_f()) == 0);
assert(D.curve().get_pdefpol()->mod_eval(E.get_g(), E.get_f()) == 0);
BivPol h = Interpolate(SimpleSum(D.get_pos(), E));
BivPol h = Interpolate(D.get_pos(), E);
EffectiveDivisor Dp = PrincipalDivisor(D.curve(), h, E);
EffectiveDivisor Dp2 = PositiveDifference(Dp, D.get_pos());
EffectiveDivisor Dp3 = PositiveDifference(Dp2, E);
EffectiveDivisor newDneg = SimpleSum(Sum(Dp3, D.get_neg()), E);
return RRspace(h, EffectiveBasisRRinDegree(newDneg, h.degree()), newDneg);
EffectiveDivisor newDneg = Sum(Dp2, D.get_neg());
return RRspace(h, EffectiveBasisRRinDegree(newDneg, E, h.degree()), newDneg);
}
......@@ -32,10 +32,10 @@ Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
// Given an effective divisor D, compute a form h on the curve such that (h) >= D
BivPol
Interpolate(const EffectiveDivisor& D);
Interpolate(const EffectiveDivisor& D, const EffectiveDivisor& E);
NTL::mat_ZZ_p
BuildInterpolationMatrixInDegree(const EffectiveDivisor& D, std::size_t deg);
BuildInterpolationMatrixInDegree(const EffectiveDivisor& D, const EffectiveDivisor& E, std::size_t deg);
// Given a form h on the curve, compute the principal effective divisor (h)
EffectiveDivisor
......@@ -55,14 +55,14 @@ MultiplyByInt(const EffectiveDivisor& D, unsigned int k);
Divisor
Sum(const Divisor& D1, const Divisor& D2);
Divisor
SimpleSum(const Divisor& D1, const Divisor& D2);
// Divisor
// SimpleSum(const Divisor& D1, const Divisor& D2);
Divisor
MultiplyByInt(const Divisor& D, unsigned int k);
std::vector<BivPol>
EffectiveBasisRRinDegree(const EffectiveDivisor& D, std::size_t deg);
EffectiveBasisRRinDegree(const EffectiveDivisor& D, const EffectiveDivisor& E, std::size_t deg);
struct RRspace {
BivPol denom;
......
Attach("magma_utils.m");
K := GF(NextPrime(2^32));
print "Disclaimer: the software rrspace implements a probabilistic algorithm. Consequently, it is possible that some of the tests fail. This is not necessarily a bug. Failure should happen very rarely (in particular the probability of failing is proportionnaly inverse to the size of the finite field), and lauching again the software should solve the problem.";
K := GF(NextPrime(2^31));
R<Y,X,Z> := PolynomialRing(K, 3);
curve_degree := 10;
print "# 16 bits finite field";
print "# smooth curve of degree 10, input divisor has no multiplicity.";
for jj in [4..100] do
Dpos_degrees := {* 10^^jj *};
repeat
......@@ -18,7 +24,7 @@ for jj in [4..100] do
print jj*10, t1, t2;
end for;
print "multiples of a point";
print "# smooth curve of degree 10, input divisor is a multiple of a rational point.";
for jj in [4..100] do
repeat
......@@ -32,7 +38,7 @@ for jj in [4..100] do
print jj*10, t1, t2;
end for;
print "nodal curve";
print "# nodal curve of degree 10, input divisor has no multiplicity.";
Q := - Y^2*Z^8 + X^2*Z^8 + Y^4*Z^6 -X^3*Z^7+X^10-5*Y^10+3*X^3*Y^7;
curve := Curve(ProjectiveSpace(K, 2), Q);
......@@ -44,3 +50,21 @@ for jj in [4..100] do
print jj*10, t1, t2;
end for;
print "# smooth curve of degree 10 over a 32 bits finite field.";
K := GF(NextPrime(2^31));
R<Y,X,Z> := PolynomialRing(K, 3);
curve_degree := 10;
for jj in [4..100] do
Dpos_degrees := {* 10^^jj *};
repeat
biv := &+[Random(K)*m : m in MonomialsOfDegree(R, curve_degree)];
until IsIrreducible(biv) and Dimension(Ideal([Derivative(biv, i) : i in [1..3]])) eq 0;
PP := ProjectiveSpace(K, 2);
RC := CoordinateRing(PP);
curve := Curve(PP, biv);
D_pos := &+[Divisor(RandomPlace(curve, d)) : d in Dpos_degrees];
t1, t2 := Bench_rrspace(curve, D_pos, false);
print jj*10, t1, t2;
end for;
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment