Skip to content
This repository has been archived by the owner on Feb 1, 2023. It is now read-only.

Commit

Permalink
Trac #32841: drop zn_poly interval_products() implementation.
Browse files Browse the repository at this point in the history 
Our interval_products_wrapper() wrapped three different implementations:

  * ntl_interval_products (ZZ_p version)
  * ntl_interval_products (zz_p version)
  * zn_poly_interval_products

But "the zn_poly version occasionally fails," and when that happens,
the NTL implementation is used anyway. Moreover the wrapper takes a
"force_ntl" parameter that allows the user to skip the zn_poly attempt
entirely.

This is the only place that zn_poly is actually used in sagelib
anymore, so this commit removes the underlying zn_poly implementation,
and eliminates the "force_ntl" parameter that was propagated
throughout the "hypellfrob" code. Now NTL is always used, as if
force_ntl was enabled.
  • Loading branch information
@orlitzky @dimpase
orlitzky authored and dimpase committed Dec 7, 2022
1 parent c938e59 commit abc091e
Show file tree
Hide file tree
Showing 6 changed files with 18 additions and 795 deletions.
1 change: 0 additions & 1 deletion .gitignore
View file
Original file line number Diff line number Diff line change
Expand Up @@ -157,7 +157,6 @@ __pycache__/
!/src/sage/rings/polynomial/weil/power_sums.c
!/src/sage/schemes/hyperelliptic_curves/hypellfrob/hypellfrob.cpp
!/src/sage/schemes/hyperelliptic_curves/hypellfrob/recurrences_ntl.cpp
!/src/sage/schemes/hyperelliptic_curves/hypellfrob/recurrences_zn_poly.cpp
!/src/sage/stats/distributions/dgs_bern.c
!/src/sage/stats/distributions/dgs_gauss_dp.c
!/src/sage/stats/distributions/dgs_gauss_mp.c
Expand Down
10 changes: 5 additions & 5 deletions src/sage/schemes/hyperelliptic_curves/hypellfrob.pyx
View file
Original file line number Diff line number Diff line change
@@ -1,8 +1,8 @@
# distutils: language = c++
# distutils: sources = sage/schemes/hyperelliptic_curves/hypellfrob/hypellfrob.cpp sage/schemes/hyperelliptic_curves/hypellfrob/recurrences_ntl.cpp sage/schemes/hyperelliptic_curves/hypellfrob/recurrences_zn_poly.cpp
# distutils: depends = sage/schemes/hyperelliptic_curves/hypellfrob/hypellfrob.h sage/schemes/hyperelliptic_curves/hypellfrob/recurrences_ntl.h sage/schemes/hyperelliptic_curves/hypellfrob/recurrences_zn_poly.h
# distutils: sources = sage/schemes/hyperelliptic_curves/hypellfrob/hypellfrob.cpp sage/schemes/hyperelliptic_curves/hypellfrob/recurrences_ntl.cpp
# distutils: depends = sage/schemes/hyperelliptic_curves/hypellfrob/hypellfrob.h sage/schemes/hyperelliptic_curves/hypellfrob/recurrences_ntl.h
# distutils: include_dirs = sage/libs/ntl/ sage/schemes/hyperelliptic_curves/hypellfrob/ NTL_INCDIR
# distutils: libraries = gmp NTL_LIBRARIES zn_poly
# distutils: libraries = gmp NTL_LIBRARIES
# distutils: extra_compile_args = NTL_CFLAGS
# distutils: library_dirs = NTL_LIBDIR
# distutils: extra_link_args = NTL_LIBEXTRA
Expand Down Expand Up @@ -59,7 +59,7 @@ cdef extern from "hypellfrob.h":
void interval_products_wrapper \
"hypellfrob::hypellfrob_interval_products_wrapper" \
(mat_ZZ_p_c &output, const mat_ZZ_p_c &M0, const mat_ZZ_p_c &M1,
const vector[ZZ_c] target, int force_ntl)
const vector[ZZ_c] target)


def interval_products(M0, M1, target):
Expand Down Expand Up @@ -146,7 +146,7 @@ def interval_products(M0, M1, target):
numintervals = len(target)/2
out.SetDims(dim, dim*numintervals)

interval_products_wrapper(out, mm0, mm1, targ, 1)
interval_products_wrapper(out, mm0, mm1, targ)
sig_off()

R = M0.matrix_space()
Expand Down
75 changes: 9 additions & 66 deletions src/sage/schemes/hyperelliptic_curves/hypellfrob/hypellfrob.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,7 +26,6 @@

#include "hypellfrob.h"
#include "recurrences_ntl.h"
#include "recurrences_zn_poly.h"

#include <cassert>

Expand Down Expand Up @@ -75,76 +74,21 @@ The intervals are supplied in "target", simply as the list
There are three possible underlying implementations:
* ntl_interval_products (ZZ_p version),
* ntl_interval_products (zz_p version),
* zn_poly_interval_products.
This function is a wrapper which takes ZZ_p input, calls one of the three
This function is a wrapper which takes ZZ_p input, calls one of the two
above implementations depending on the size of the current ZZ_p modulus, and
produces output in ZZ_p format.
If the force_ntl flag is set, it will never use the zn_poly version.
Note that the zn_poly version occasionally fails; this happens more frequently
for smaller p, but is extremely rare for larger p. This wrapper detects this
and falls back on the zz_p/ZZ_p versions, which should never fail.
PRECONDITIONS:
Let d = b[n-1] - a[0]. Then 2, 3, ... 1 + floor(sqrt(d)) must all be
invertible.
*/
void interval_products_wrapper(vector<mat_ZZ_p>& output,
const mat_ZZ_p& M0, const mat_ZZ_p& M1,
const vector<ZZ>& target,
int force_ntl = 0)
const vector<ZZ>& target)
{
const ZZ& modulus = ZZ_p::modulus();

if (!force_ntl && modulus <= (1UL << (ULONG_BITS - 1)) - 1)
{
// Small modulus; let's try using zn_poly if we're allowed.
// (we don't allow the full ULONG_BITS bits, because I'm worried I
// sometimes use NTL code which requires longs rather than ulongs.)
zn_mod_t mod;
zn_mod_init(mod, to_ulong(modulus));

int r = M0.NumRows();

// convert to zn_mod format
vector<vector<ulong> > M0_temp(r, vector<ulong>(r));
vector<vector<ulong> > M1_temp(r, vector<ulong>(r));
vector<vector<vector<ulong> > > output_temp(target.size()/2, M0_temp);

for (int x = 0; x < r; x++)
for (int y = 0; y < r; y++)
{
M0_temp[y][x] = to_ulong(rep(M0[y][x]));
M1_temp[y][x] = to_ulong(rep(M1[y][x]));
}

int success = zn_poly_interval_products(output_temp, M0_temp, M1_temp,
target, mod);
zn_mod_clear(mod);

// note: if we failed, we fall back on the ZZ_p or zz_p versions below
if (success)
{
// convert back to ZZ_p format
output.clear();
mat_ZZ_p temp;
temp.SetDims(r, r);
for (size_t i = 0; i < target.size()/2; i++)
{
for (int x = 0; x < r; x++)
for (int y = 0; y < r; y++)
temp[y][x] = output_temp[i][y][x];
output.push_back(temp);
}

return;
}

// failed
}

if (!modulus.SinglePrecision())
{
// Large modulus.... go straight to the ZZ_p version
Expand Down Expand Up @@ -183,11 +127,10 @@ void interval_products_wrapper(vector<mat_ZZ_p>& output,

void hypellfrob_interval_products_wrapper(mat_ZZ_p& output,
const mat_ZZ_p& M0, const mat_ZZ_p& M1,
const vector<ZZ>& target,
int force_ntl = 0)
const vector<ZZ>& target)
{
vector<mat_ZZ_p> mat_vector;
interval_products_wrapper(mat_vector, M0, M1, target, force_ntl);
interval_products_wrapper(mat_vector, M0, M1, target);
int r = M0.NumRows();
output.SetDims(r, r * mat_vector.size());

Expand Down Expand Up @@ -328,7 +271,7 @@ void padic_invert_matrix(mat_ZZ_p& B, const mat_ZZ_p& A, const ZZ& p, int N)
/*
The main function exported from this module. See hypellfrob.h for information.
*/
int matrix(mat_ZZ& output, const ZZ& p, int N, const ZZX& Q, int force_ntl)
int matrix(mat_ZZ& output, const ZZ& p, int N, const ZZX& Q)
{
// check input conditions
if (N < 1 || p < 3)
Expand Down Expand Up @@ -522,8 +465,8 @@ int matrix(mat_ZZ& output, const ZZ& p, int N, const ZZX& Q, int force_ntl)
vector<mat_ZZ_p> MH, DH;
MH.reserve(L);
DH.reserve(L);
interval_products_wrapper(MH, MH0[j], MH1[j], s, force_ntl);
interval_products_wrapper(DH, DH0[j], DH1[j], s, force_ntl);
interval_products_wrapper(MH, MH0[j], MH1[j], s);
interval_products_wrapper(DH, DH0[j], DH1[j], s);

// Use the vandermonde matrix to extend all the way up to L if necessary.
if (L > N)
Expand Down Expand Up @@ -710,8 +653,8 @@ int matrix(mat_ZZ& output, const ZZ& p, int N, const ZZX& Q, int force_ntl)
}
modulusV.restore();
vector<mat_ZZ_p> MV, DV;
interval_products_wrapper(MV, MV0, MV1, s, force_ntl);
interval_products_wrapper(DV, DV0, DV1, s, force_ntl);
interval_products_wrapper(MV, MV0, MV1, s);
interval_products_wrapper(DV, DV0, DV1, s);

// Divide out each MV[j] by DV[j]. Note that for 1 <= j < N, both of them
// should be divisible by p too, so we take that into account here.
Expand Down
17 changes: 4 additions & 13 deletions src/sage/schemes/hyperelliptic_curves/hypellfrob/hypellfrob.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -44,27 +44,19 @@ The intervals are supplied in "target", simply as the list
There are three possible underlying implementations:
* ntl_interval_products (ZZ_p version),
* ntl_interval_products (zz_p version),
* zn_poly_interval_products.
This function is a wrapper which takes ZZ_p input, calls one of the three
* ntl_interval_products (zz_p version)
This function is a wrapper which takes ZZ_p input, calls one of the two
above implementations depending on the size of the current ZZ_p modulus, and
produces output in ZZ_p format.
If the force_ntl flag is set, it will never use the zn_poly version.
Note that the zn_poly version occasionally fails; this happens more frequently
for smaller p, but is extremely rare for larger p. This wrapper detects this
and falls back on the zz_p/ZZ_p versions, which should never fail.
PRECONDITIONS:
Let d = b[n-1] - a[0]. Then 2, 3, ... 1 + floor(sqrt(d)) must all be
invertible.
*/
void hypellfrob_interval_products_wrapper(NTL::mat_ZZ_p& output,
const NTL::mat_ZZ_p& M0, const NTL::mat_ZZ_p& M1,
const std::vector<NTL::ZZ>& target,
int force_ntl);
const std::vector<NTL::ZZ>& target);

/*
Computes frobenius matrix for given p, to precision p^N, for the
Expand All @@ -82,8 +74,7 @@ RETURN VALUE:
will not check that p is prime. That's up to you.)
*/
int matrix(NTL::mat_ZZ& output, const NTL::ZZ& p, int N, const NTL::ZZX& Q,
int force_ntl = 0);
int matrix(NTL::mat_ZZ& output, const NTL::ZZ& p, int N, const NTL::ZZX& Q);


};
Expand Down
Loading

0 comments on commit abc091e

Please sign in to comment.