initial commit

parent 53e708ad
Pierre-Jean Spaenlehauer (pierre-jean.spaenlehauer@inria.fr)
This software implements the Brill-Noether algorithm for computing Riemann-Roch
spaces and is based on the work done by Aude Le Gluher in her Master's thesis.
The files ZZ_pXResultant.cc and ZZ_pXResultant.h have been written by Pierrick
Gaudry and are based on the NTL code for computing resultants. We thank
Pierrick Gaudry for allowing us to use these functions and to include them in
this software.
make rrspace class_grp_arith
./rrspace < example_rrspace_input
./class_grp_arith < example_class_grp_div_input
This diff is collapsed.
CPPFLAGS=--pedantic -Wall --std=c++11 -O3
LDLIBS=-lntl -lgmpxx -lgmp -lm
all: rrspace class_grp_arith tests
clean:
rm *.o rrspace class_grp_arith tests
algos.o: bivpol.h divisor.h curve.h algos.h algos.cc
bivpol.o: bivpol.h bivpol.cc
divisor.o: divisor.h curve.h bivpol.h divisor.cc
divisor_class_group.o: divisor.h curve.h bivpol.h divisor.h divisor_class_group.h divisor_class_group.cc
ZZ_pXResultant.o: ZZ_pXResultant.h ZZ_pXResultant.cc
rrspace: rrspace.cc algos.o bivpol.o ZZ_pXResultant.o divisor.o
class_grp_arith: class_grp_arith.cc algos.o bivpol.o ZZ_pXResultant.o divisor.o divisor_class_group.o
tests: tests.cc algos.o bivpol.o ZZ_pXResultant.o divisor.o divisor_class_group.o
$(CXX) $(CPPFLAGS) $^ $(LDLIBS) -lcppunit -o $@
.PHONY: all clean cleanall
This software is distributed under LGPL-2.1+ (i.e., LGPL version 2.1 or later).
Please report bugs to <pierre-jean.spaenlehauer@inria.fr>.
I - Installation
----------------
This software can be compiled using GNU make with g++.
Running the command
make
in the directory containing the source files builds three executable files:
rrspace, class_grp_arith and tests.
In order to build only one of these programs, use the command
make <PROGRAM>
where <PROGRAM> is rrspace, class_grp_arith or tests
Compiling this software requires the NTL library (https://www.shoup.net/ntl/).
In order to compile the program "tests", the library cppunit is also required
(https://freedesktop.org/wiki/Software/cppunit/).
After compiling, the three following commands should succeed without any error:
./tests
./rrspace < example_rrspace_input
./class_grp_arith < example_class_grp_div_input
II - Usage
----------
Examples of inputs for rr_space and class_grp_arith are provided in the files
example_rrspace_input and example_class_grp_div_input. The formats for the
mathematical objects are detailed in the next section.
Disclaimer: The algorithms implemented in this software are probabilistic, and
they have a probability of failure. Also, they rely on some genericity
assumptions (such as generic coordinates). Finally, some problems may arise in
small characteristic, for instance if the characteristic divides the degree of
the curve.
In case of failure, the most probable situation is that just starting again the
software with the same input will work. In the unlikely situation that it still
fails, changing the system of coordinates by doing a random linear change of
projective coordinates and by computing the curve equation and the divisors in
this new coordinate system may work. Finally, if none of this works, this means
that we are in the very unlikely case where all the possible denominators which
could be returned by the interpolation function vanish at singular points of
the curve (this case cannot occur if the curve is smooth). In this case, please
use another model of the curve.
1) rrspace: the program rrspace computes the Riemann-Roch space associated to a
given divisor D on a curve q. It reads its input on the standard input, and it
returns a basis of the Riemann-Roch space on the standard output.
The format for its input (see next section for more details) is
<PRIME>
<CURVE>
<DIVISOR>
where <DIVISOR> represents a divisor whose support does not contain any
singular point of the curve.
It outputs the following data:
- the dimension D of the Riemann-Roch space;
- a denominator h, which is a bivariate polynomial;
- a list of D bivariate polynomials g_1, ..., g_D
The set of rational functions {g_1/h, ..., g_D/h} is a basis of the
Riemann-Roch space.
2) class_grp_arith: this program takes as input two elements of the divisor
class group of a curve C of genus g and it returns a representation of the sum. Here, we
are given a curve, together with a degree 1 divisor on C. A class of divisors
(modulo principal divisors) is represented by an effective divisor D1 of degree
g: it represents the class of divisors equivalent to the degree-0 divisor D1 -
g*O.
class_grp_arith reads its input on the standard input, and it returns an
effective divisor of degree at most g on the standard output.
The format for its input (see next section for more details) is
<PRIME>
<CURVE>
<GENUS>
<DIVISOR>
<EFFECTIVEDIVISOR>
<EFFECTIVEDIVISOR>
The divisor provided on the 4th line must have degree 1. The two effective
divisors D1 and D2 provided on the last two lines must have degree <GENUS>.
The supports of all the divisors must not contain any singular point of the
curve.
The program returns a representation of a degree-g effective divisor D3 such that
(D1-g*O) + (D2-g*O) = (D3-g*O)
in the divisor class group of C.
3) tests: this software does not take any input. It just runs two tests: one
for the computation of Riemann-Roch spaces, and one for the arithmetic in the
divisor class group of a smooth quartic curve. It should output the string
"OK".
III - Data structures
---------------------
This section details the data structures for the input of the programs
described in Section II.
<PRIME>: a prime number p
<CURVE>: a bivariate polynomial q = sum_{i=0}^d q_i(x)*y^i, where for all i,
q_i is a polynomial of degree at most d-i. The program requires that q_d(x) = 1.
The polynomial q is encoded as follows
[ <DEGREE> [ <q_0> <q_1> ... <q_d> ] ]
where <DEGREE> is the degree of the curve and <q_i> is the list of coefficients
of the polynomial q_i. For instance if q_i(x) = sum_{j=0}^{r} q_{i,j} x^j,
then <q_i> should be
[ q_{i,0} q_{i,1} ... q_{i,r} ]
<GENUS>: the genus of the curve C described by the polynomial q. We decline all
responsability on the behavior of the programs if g is not the genus of the
curve.
<EFFECTIVEDIVISOR>: an effective divisor on C represented by two univariate
polynomials
< <POLY1> <POLY2> >
Each polynomial is described by the list of its coefficients
[ f_0 f_1 ... f_l ]
These two polynomials (referred to as f(X), g(X)) should satisfy the following
equality: q(X, g(X)) mod f(X) = 0.
They represent the effective divisor which is the sum of points (\alpha, \beta)
with multiplicity \gamma where \alpha is a root of multiplicity \gamma of f(X)
and \beta = g(\alpha).
<DIVISOR>: a divisor D on the curve C represented by two effective divisors D_+
and D_- such that D = D_+ - D_-.
{ <EFFECTIVEDIVISOR> <EFFECTIVEDIVISOR> }
# rrspace
The rrspace software computes bases of Riemann-Roch spaces for curves defined over Z/pZ.
It also provides functions to compute the group law on the Jacobain of such curves.
\ No newline at end of file
/* Common header file for the rrspace software
This file is part of the rrspace project.
This library is free software; you can redistribute it and/or modify it under
the terms of the GNU Lesser General Public License as published by the Free
Software Foundation; either version 2.1 of the License, or (at your option) any
later version.
This library is distributed in the hope that it will be useful, but WITHOUT ANY
WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
details.
You should have received a copy of the GNU Lesser General Public License along
with this library; if not, write to the Free Software Foundation, Inc., 51
Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
*/
/* This file has been written by Pierrick Gaudry <pierrick.gaudry@loria.fr>
* */
#include <NTL/ZZ_pX.h>
#include <NTL/new.h>
#include <assert.h>
using namespace std;
// include code from NTL's official ZZ_pX.c
NTL_OPEN_NNS
void mul(ZZ_pX& U, ZZ_pX& V, const ZZ_pXMatrix& M);
void mul(ZZ_pXMatrix& A, ZZ_pXMatrix& B, ZZ_pXMatrix& C);
void PlainResultant(ZZ_p& rres, const ZZ_pX& a, const ZZ_pX& b);
void ResIterHalfGCD(ZZ_pXMatrix& M_out, ZZ_pX& U, ZZ_pX& V, long d_red,
vec_ZZ_p& cvec, vec_long& dvec);
void ResHalfGCD(ZZ_pXMatrix& M_out, const ZZ_pX& U, const ZZ_pX& V,
long d_red, vec_ZZ_p& cvec, vec_long& dvec);
void ResHalfGCD(ZZ_pX& U, ZZ_pX& V, vec_ZZ_p& cvec, vec_long& dvec);
NTL_CLOSE_NNS
NTL_START_IMPL
// Euclidian algorithm between U and V
// Fill in cvec and dvec with degrees and LeadCoeffs of all the
// polynomials, starting with V and finishing with the last non-zero
// poly.
// U and V are modified: at the end, V is 0 and U is the previous poly
// (usually a constant). UU is filled with the ante-previous poly
// (usually of degree 1: the last subresultant up to a mult. factor)
void PlainResultant(ZZ_pX& UU, ZZ_pX& U, ZZ_pX& V, vec_ZZ_p& cvec,
vec_long& dvec) {
ZZ_pX Q;
ZZVec tmp(deg(U)+1, ZZ_pInfo->ExtendedModulusSize);
while (deg(V) >= 0) {
append(cvec, LeadCoeff(V));
append(dvec, deg(V));
UU = U;
PlainDivRem(Q, U, U, V, tmp);
swap(U, V);
}
}
// returns the resultant and the last non-trivial sub-resultant
// rres is the resultant
// subres1 is the last subresultant (of degree <= 1)
// subres2 is filled if the subres1 is degenerated (of degree 0). Then
// subres2 is the last non-trivial subresultant.
// All of this might fail, in the sense that subres1 and subres2 could
// be 0 instead of having meaningful values. In particular, if rres = 0,
// no further information is obtained. However, if a subres is <> 0, then
// it is (should be?) exact.
//
// Warning: subres1 and subres2 are not allowed to be aliases of the
// input parameters.
void PlainResultantWithSubRes(ZZ_p& rres, ZZ_pX& subres1, ZZ_pX& subres2,
const ZZ_pX& a, const ZZ_pX& b)
{
ZZ_p res, saved_res;
clear(subres1);
clear(subres2);
if (IsZero(a) || IsZero(b))
clear(res);
else if (deg(a) == 0 && deg(b) == 0)
set(res);
else {
long d0, d1, d2;
long sg = 1;
ZZ_p lc;
set(res);
long n = max(deg(a),deg(b)) + 1;
ZZ_pX u(INIT_SIZE, n), v(INIT_SIZE, n), saved_u(INIT_SIZE, n);
ZZVec tmp(n, ZZ_pInfo->ExtendedModulusSize);
u = a;
v = b;
for (;;) {
d0 = deg(u);
d1 = deg(v);
lc = LeadCoeff(v);
saved_u = u;
PlainRem(u, u, v, tmp);
swap(u, v);
d2 = deg(v);
if (d2 >= 0) {
power(lc, lc, d0-d2);
saved_res = res;
mul(res, res, lc);
if (d0 & d1 & 1) negate(res, res);
if (!((d0-d1) & 1)) {
sg = -sg;
cerr << "sg = " << sg << endl;
}
// cerr << saved_res*u << endl;
}
else {
if (d1 == 0) {
cerr << "d0 = " << d0 << endl;
subres1 = sg*saved_res*saved_u;
if (deg(subres1) > 1) { // degenerated case
subres2 = subres1;
if (deg(subres1) > 2)
clear(subres1);
else
subres1 = saved_res*lc*power(LeadCoeff(saved_u),d0);
}
power(lc, lc, d0);
mul(res, res, lc);
}
else {
subres1 = sg*res*u;
if (deg(subres1) > 1) { // *very* degenerated case
cerr << "WARNING: This case is not fully implemented!!!\n";
subres2 = subres1;
if (deg(subres1) > 2)
clear(subres1);
else
subres1 = res*lc*power(LeadCoeff(u),d0);
}
clear(res);
}
break;
}
}
rres = res;
}
}
static void Coef(ZZ_p& Coe, vec_ZZ_p& cvec, vec_long& dvec) {
int r = cvec.length();
ZZ_p t;
set(Coe);
for (int l = 3; l <= r-1; ++l) {
power(t, cvec[l-2], dvec[l-3]-dvec[l-1]);
mul(Coe, Coe, t);
}
}
static void Signs(long& sig1, long& sig2, vec_long& dvec) {
int r = dvec.length();
sig1 = 0; sig2 = 0;
for (int l = 3; l <= r-1; ++l) {
sig1 = (sig1 + dvec[l-3]*(dvec[l-2] & 1L)) & 1L;
sig2 = (sig2 + 1L + dvec[l-3] + dvec[l-2]) & 1L;
}
}
// For the moment, I assume that deg(u) > deg(v) and v1 <> 0
void resultantWithSubRes(ZZ_p& rres, ZZ_pX& subres,
const ZZ_pX& u, const ZZ_pX& v)
{
#if 0
if (deg(u) <= NTL_ZZ_pX_GCD_CROSSOVER || deg(v) <= NTL_ZZ_pX_GCD_CROSSOVER) {
PlainResultantWithSubRes(rres, subres1, subres2, u, v);
return;
}
#endif
ZZ_pX u1, v1;
u1 = u;
v1 = v;
ZZ_p t;
ZZ_pX res;
#if 0
clear(subres1);
clear(subres2);
if (deg(u1) == deg(v1)) {
rem(u1, u1, v1);
swap(u1, v1);
if (IsZero(v1)) {
clear(rres);
return;
}
power(t, LeadCoeff(u1), deg(u1) - deg(v1));
mul(res, res, t);
if (deg(u1) & 1)
negate(res, res);
}
else if (deg(u1) < deg(v1)) {
swap(u1, v1);
if (deg(u1) & deg(v1) & 1)
negate(res, res);
}
#endif
assert ((deg(u1) >= deg(v1)) && (v1 != 0));
vec_ZZ_p cvec;
vec_long dvec;
cvec.SetMaxLength(deg(v1)+2);
dvec.SetMaxLength(deg(v1)+2);
append(cvec, LeadCoeff(u1));
append(dvec, deg(u1));
while (deg(u1) > NTL_ZZ_pX_GCD_CROSSOVER && !IsZero(v1)) {
ResHalfGCD(u1, v1, cvec, dvec);
if (!IsZero(v1)) {
append(cvec, LeadCoeff(v1));
append(dvec, deg(v1));
rem(u1, u1, v1);
swap(u1, v1);
}
}
if (IsZero(v1) && deg(u1) > 0) {
clear(rres);
return;
}
if (deg(u1) <= 1) {
// we went all the way...
// Too Bad!!! we have skipped the interesting subresultant.
// ==> in this rare case, we go back to the naive algorithm
// (in principle we should detect exactly which step of the
// recursive call went too fast and to redo just that one... but
// this is too much work for something that nether (?) happens.)
u1 = u; v1 = v;
cvec.SetLength(0);
dvec.SetLength(0);
append(cvec, LeadCoeff(u1));
append(dvec, deg(u1));
}
ZZ_pX UU;
ZZ_p Coe;
long sig1, sig2;
PlainResultant(UU, u1, v1, cvec, dvec);
long r = dvec.length();
assert (cvec.length() == dvec.length());
Coef(Coe, cvec, dvec);
Signs(sig1, sig2, dvec);
power(t, cvec[r-2], dvec[r-3] - dvec[r-2] - 1);
mul(Coe, Coe, t);
mul(subres, Coe, UU);
if ((sig1 + dvec[r-2]*sig2) & 1UL)
negate(subres, subres);
sig1 = (sig1 + dvec[r-3]*(dvec[r-2] & 1L)) & 1L;
sig2 = (sig2 + 1L + dvec[r-3] + dvec[r-2]) & 1L;
power(t, cvec[r-2], dvec[r-2] - dvec[r-1] + 1);
mul(Coe, Coe, t);
mul(res, u1, Coe);
power(t, cvec[r-1], dvec[r-2] - dvec[r-1] - 1);
mul(res, res, t);
if ((sig1 + dvec[r-1]*sig2) & 1L)
negate(res, res);
if (deg(res) >= 1) {
subres = res;
res = 0;
if (deg(subres) >= 2)
subres = 0;
} else if (deg(subres) > 2) {
subres = 0;
} else if (deg(subres) == 2) {
mul(subres, Coe, u1);
power(t, cvec[r-2], 1-(dvec[r-2] - dvec[r-1]));
mul(subres, subres, t);
if ((sig1 + (dvec[r-2]+1)*sig2) & 1L)
negate(subres, subres);
}
rres = coeff(res,0);
}
NTL_END_IMPL
/* Common header file for the rrspace software
This file is part of the rrspace project.
This library is free software; you can redistribute it and/or modify it under
the terms of the GNU Lesser General Public License as published by the Free
Software Foundation; either version 2.1 of the License, or (at your option) any
later version.
This library is distributed in the hope that it will be useful, but WITHOUT ANY
WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
details.
You should have received a copy of the GNU Lesser General Public License along
with this library; if not, write to the Free Software Foundation, Inc., 51
Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
*/
/* This file has been written by Pierrick Gaudry <pierrick.gaudry@loria.fr>
* */
#ifndef NTL_ZZ_pXResultant__H
#define NTL_ZZ_pXResultant__H
#include <NTL/ZZ.h>
#include <NTL/ZZ_p.h>
#include <NTL/ZZ_pX.h>
NTL_OPEN_NNS
void resultantWithSubRes(ZZ_p& rres, ZZ_pX& subres,
const ZZ_pX& u, const ZZ_pX& v);
NTL_CLOSE_NNS
#endif
/* Common header file for the rrspace software
This file is part of the rrspace project.
This library is free software; you can redistribute it and/or modify it under
the terms of the GNU Lesser General Public License as published by the Free
Software Foundation; either version 2.1 of the License, or (at your option) any
later version.
This library is distributed in the hope that it will be useful, but WITHOUT ANY
WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
PARTICULAR PURPOSE. See the GNU Lesser General Public License for more
details.
You should have received a copy of the GNU Lesser General Public License along
with this library; if not, write to the Free Software Foundation, Inc., 51
Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
*/
#include "algos.h"
#include <utility>
#include <vector>
#include "ZZ_pXResultant.h"
using namespace NTL;
using namespace std;
size_t
GetInterpolationDegree(size_t deg_curve, size_t deg_div) {
if ((deg_curve + 1)*deg_curve > 2*deg_div)
return floor((1+sqrt(1+8*deg_div)-1)/2);
return floor((double)deg_div/(double)deg_curve+((double)deg_curve-1)/2);
}
// Returns a vector such that the ind-th element is the pair (i, j) which
// corresponds to the ind-th monomial X^i Y^j such that j < deg_curve
vector<pair<size_t, size_t> >
BuildMonomialVector(size_t deg_curve, size_t deg_max) {
vector<pair<size_t, size_t> > res;
for (size_t d = 0; d <= deg_max; ++d)
for (size_t i = 0; i <= min(d, deg_curve - 1); ++i)
res.push_back(pair<size_t, size_t>(d-i, i));
return res;
}
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());
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);
}
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);
}
auto vec_mons = BuildMonomialVector(D.curve().degree(), deg);
Vec<Vec<ZZ_p> > vec_vec_res;
Vec<ZZ_p> tmp_row;
ZZ_pX tmp_pol;
for (const auto& p : vec_mons) {
tmp_row.SetLength(0);
MulMod(tmp_pol, powers_of_x[p.first], powers_of_g[p.second], modf);
for (size_t i = 0; i < D.degree(); ++i)
tmp_row.append(coeff(tmp_pol, i));
vec_vec_res.append(tmp_row);
}
mat_ZZ_p res;
res.SetDims(vec_mons.size(), D.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);
mat_ZZ_p kerM;
kernel(kerM, M);
mat_ZZ_p randmat;
randmat.SetDims(1, kerM.NumRows());
for (size_t i = 0; i < (size_t)kerM.NumRows(); ++i)
randmat.put(0, i, random_ZZ_p());
mul(kerM, randmat, kerM);
vector<ZZ_pX> interp_pol;
auto vec_mons = BuildMonomialVector(D.curve().degree(), deg_res);
interp_pol.insert(interp_pol.begin(),
min(D.curve().degree(), deg_res + 1),
ZZ_pX(0));
for (size_t i = 0; i < vec_mons.size(); ++i) {
pair<size_t, size_t> p = vec_mons[i];
interp_pol[p.second] += kerM.get(0, i)* (ZZ_pX(1) << p.first);
}
while (IsZero(interp_pol.back()))
interp_pol.pop_back();
interp_pol.shrink_to_fit();
BivPol res(interp_pol);
assert(!res.IsZero());
return res;
}
EffectiveDivisor
PrincipalDivisor(const Curve& C, const BivPol& h) {
size_t deg_res = h.degree() * C.degree(); // Bézout bound
vec_ZZ_p eval_pts;
for (long i = 0; i <= (long)deg_res; ++i)
eval_pts.append(ZZ_p(i));
vec_vec_ZZ_p eval_C, eval_h;
vec_ZZ_p eval_res, eval_subres0, eval_subres1;
for (size_t i = 0; i <= C.degree(); ++i)
eval_C.append(eval(C.coeff(i), eval_pts));
for (size_t i = 0; i <= h.degree_y(); ++i)
eval_h.append(eval(h.coeff(i), eval_pts));
for (size_t i = 0; i <= deg_res; ++i) {
ZZ_pX C_tmp, h_tmp, subres_tmp;
ZZ_p res_tmp;
for (size_t j = 0; j <= C.degree(); ++j)
C_tmp += eval_C[j][i] * (ZZ_pX(1) << j);
for (size_t j = 0; j <= h.degree_y(); ++j)
h_tmp += eval_h[j][i] * (ZZ_pX(1) << j);
resultantWithSubRes(res_tmp, subres_tmp, C_tmp, h_tmp);
assert(deg(subres_tmp) == 1);
eval_res.append(res_tmp);
eval_subres0.append(ConstTerm(subres_tmp));
eval_subres1.append(LeadCoeff(subres_tmp));
}
ZZ_pX res = interpolate(eval_pts, eval_res);
ZZ_pX subres0 = interpolate(eval_pts, eval_subres0);
ZZ_pX subres1 = interpolate(eval_pts, eval_subres1);
return EffectiveDivisor(C, res, -MulMod(InvMod(subres1, res), subres0, res));
}
EffectiveDivisor
PositiveDifference(const EffectiveDivisor& D1, const EffectiveDivisor& D2) {
assert(D1.curve() == D2.curve());
ZZ_pX f;
divide(f, D1.get_f(), GCD(D1.get_f(), D2.get_f()));
return EffectiveDivisor(D1.curve(), f/LeadCoeff(f), D1.get_g() % f);
}
EffectiveDivisor
Sum(const EffectiveDivisor& D1, const EffectiveDivisor& D2) {
assert(D1.curve() == D2.curve());
assert(deg(GCD(D1.get_f(), D2.get_f())) == 0);
if (deg(D1.get_f()) == 0)
return D2;
if (deg(D2.get_f()) == 0)
return D1;
ZZ_pX new_f = D1.get_f()*D2.get_f();
ZZ_pX new_g = (D1.get_g()*D2.get_f()*InvMod(D2.get_f() % D1.get_f(), D1.get_f()) +
D2.get_g()*D1.get_f()*InvMod(D1.get_f() % D2.get_f(), D2.get_f())) %
(D1.get_f()*D2.get_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);
}
EffectiveDivisor
MultiplyByInt(const EffectiveDivisor& D, unsigned int k) {
ZZ_pX _f = D.get_f();
ZZ_pX _g = D.get_g();
assert(IsZero(D.curve().get_pdefpol()->mod_eval(_g, _f)));
// Case where D is the zero divisor
if (deg(_f) == 0)
return D;
int nbsteps;
frexp((double)k, &nbsteps);
// TODO: do all computations mod _f
_g = NewtonHensel(*D.curve().get_pdefpol(), _g, _f, nbsteps);
power(_f, _f, k);
return EffectiveDivisor(D.curve(), _f, _g % _f);
}
Divisor
<