diff --git a/doc/source/nmod_poly.rst b/doc/source/nmod_poly.rst
index 7e3fc174ea..cefcb29189 100644
--- a/doc/source/nmod_poly.rst
+++ b/doc/source/nmod_poly.rst
@@ -1406,99 +1406,99 @@ Multipoint evaluation
--------------------------------------------------------------------------------
-.. function:: void _nmod_poly_evaluate_nmod_vec_iter(nn_ptr ys, nn_srcptr poly, slong len, nn_srcptr xs, slong n, nmod_t mod)
+.. function:: void _nmod_poly_evaluate_nmod_vec_iter(nn_ptr ys, nn_srcptr poly, slong ilen, nn_srcptr xs, slong olen, nmod_t mod)
- Evaluates (``coeffs``, ``len``) at the ``n`` values
+ Evaluates (``coeffs``, ``ilen``) at the ``olen`` values
given in the vector ``xs``, writing the output values
to ``ys``. The values in ``xs`` should be reduced
modulo the modulus.
Uses Horner's method iteratively.
-.. function:: void nmod_poly_evaluate_nmod_vec_iter(nn_ptr ys, const nmod_poly_t poly, nn_srcptr xs, slong n)
+.. function:: void nmod_poly_evaluate_nmod_vec_iter(nn_ptr ys, const nmod_poly_t poly, nn_srcptr xs, slong olen)
- Evaluates ``poly`` at the ``n`` values given in the vector
+ Evaluates ``poly`` at the ``olen`` values given in the vector
``xs``, writing the output values to ``ys``. The values in
``xs`` should be reduced modulo the modulus.
Uses Horner's method iteratively.
-.. function:: void _nmod_poly_evaluate_nmod_vec_fast_precomp(nn_ptr vs, nn_srcptr poly, slong plen, const nn_ptr * tree, slong len, nmod_t mod)
+.. function:: void _nmod_poly_evaluate_nmod_vec_fast_precomp(nn_ptr vs, nn_srcptr poly, slong ilen, const nn_ptr * tree, slong olen, nmod_t mod)
- Evaluates (``poly``, ``plen``) at the ``len`` values given
+ Evaluates (``poly``, ``ilen``) at the ``olen`` values given
by the precomputed subproduct tree ``tree``.
-.. function:: void _nmod_poly_evaluate_nmod_vec_fast(nn_ptr ys, nn_srcptr poly, slong len, nn_srcptr xs, slong n, nmod_t mod)
+.. function:: void _nmod_poly_evaluate_nmod_vec_fast(nn_ptr ys, nn_srcptr poly, slong ilen, nn_srcptr xs, slong olen, nmod_t mod)
- Evaluates (``coeffs``, ``len``) at the ``n`` values
+ Evaluates (``coeffs``, ``ilen``) at the ``olen`` values
given in the vector ``xs``, writing the output values
to ``ys``. The values in ``xs`` should be reduced
modulo the modulus.
Uses fast multipoint evaluation, building a temporary subproduct tree.
-.. function:: void nmod_poly_evaluate_nmod_vec_fast(nn_ptr ys, const nmod_poly_t poly, nn_srcptr xs, slong n)
+.. function:: void nmod_poly_evaluate_nmod_vec_fast(nn_ptr ys, const nmod_poly_t poly, nn_srcptr xs, slong olen)
- Evaluates ``poly`` at the ``n`` values given in the vector
+ Evaluates ``poly`` at the ``olen`` values given in the vector
``xs``, writing the output values to ``ys``. The values in
``xs`` should be reduced modulo the modulus.
Uses fast multipoint evaluation, building a temporary subproduct tree.
-.. function:: void _nmod_poly_evaluate_nmod_vec(nn_ptr ys, nn_srcptr poly, slong len, nn_srcptr xs, slong n, nmod_t mod)
+.. function:: void _nmod_poly_evaluate_nmod_vec(nn_ptr ys, nn_srcptr poly, slong ilen, nn_srcptr xs, slong olen, nmod_t mod)
- Evaluates (``poly``, ``len``) at the ``n`` values
+ Evaluates (``poly``, ``ilen``) at the ``olen`` values
given in the vector ``xs``, writing the output values
to ``ys``. The values in ``xs`` should be reduced
modulo the modulus.
-.. function:: void nmod_poly_evaluate_nmod_vec(nn_ptr ys, const nmod_poly_t poly, nn_srcptr xs, slong n)
+.. function:: void nmod_poly_evaluate_nmod_vec(nn_ptr ys, const nmod_poly_t poly, nn_srcptr xs, slong olen)
- Evaluates ``poly`` at the ``n`` values given in the vector
+ Evaluates ``poly`` at the ``olen`` values given in the vector
``xs``, writing the output values to ``ys``. The values in
``xs`` should be reduced modulo the modulus.
-.. function:: void _nmod_poly_evaluate_geometric_nmod_vec_iter(nn_ptr ys, nn_srcptr coeffs, slong len, ulong r, slong n, nmod_t mod)
+.. function:: void _nmod_poly_evaluate_geometric_nmod_vec_iter(nn_ptr ys, nn_srcptr coeffs, slong ilen, ulong r, slong olen, nmod_t mod)
- Evaluates (``coeffs``, ``len``) at the first ``n`` powers
+ Evaluates (``coeffs``, ``ilen``) at the first ``olen`` powers
of the square of ``r``, writing the output values
to ``ys``. The value of ``r`` should be reduced
modulo the modulus.
Uses Horner's method iteratively.
-.. function:: void nmod_poly_evaluate_geometric_nmod_vec_iter(nn_ptr ys, const nmod_poly_t poly, ulong r, slong n)
+.. function:: void nmod_poly_evaluate_geometric_nmod_vec_iter(nn_ptr ys, const nmod_poly_t poly, ulong r, slong olen)
- Evaluates ``poly`` at the first ``n`` powers
+ Evaluates ``poly`` at the first ``olen`` powers
of the square of ``r``, writing the output values
to ``ys``. The value of ``r`` should be reduced
modulo the modulus.
Uses Horner's method iteratively.
-.. function:: void _nmod_poly_evaluate_geometric_nmod_vec_fast_precomp(nn_ptr vs, nn_srcptr poly, slong plen, const nmod_geometric_progression_t G, slong len, nmod_t mod)
+.. function:: void _nmod_poly_evaluate_geometric_nmod_vec_fast_precomp(nn_ptr vs, nn_srcptr poly, slong ilen, const nmod_geometric_progression_t G, slong olen, nmod_t mod)
- Evaluates (``poly``, ``plen``) at the ``len`` values given
- by the precomputed geometric progression ``G``. The value of
- ``len`` should be less than or equal to the precomputation size parameter ``G->len``.
+ Evaluates (``poly``, ``ilen``) at the first ``olen`` values given by the
+ precomputed geometric progression ``G``, which are the first ``olen`` powers
+ of the square of ``r``. Requires ``olen <= G->len``.
-.. function:: void _nmod_poly_evaluate_geometric_nmod_vec_fast(nn_ptr ys, nn_srcptr coeffs, slong len, ulong r, slong n, nmod_t mod)
+.. function:: void _nmod_poly_evaluate_geometric_nmod_vec_fast(nn_ptr ys, nn_srcptr coeffs, slong ilen, ulong r, slong olen, nmod_t mod)
- Evaluates (``coeffs``, ``len``) at the first ``n`` powers
- of the square of ``r``, writing the output values to ``ys``.
+ Evaluates (``coeffs``, ``ilen``) at the first ``olen`` powers
+ of the square of ``r``, writing the output values to ``ys``.
The value of ``r`` should be reduced modulo the modulus ``mod``
- and of sufficient multiplicative order such that none of
- the first `n` powers of `r^2` is one.
+ and of sufficient multiplicative order such that none of
+ the first ``olen`` powers of `r^2` is one.
Uses fast geometric multipoint evaluation, building a temporary geometric progression precomputation.
-.. function:: void nmod_poly_evaluate_geometric_nmod_vec_fast(nn_ptr ys, const nmod_poly_t poly, ulong r, slong n)
+.. function:: void nmod_poly_evaluate_geometric_nmod_vec_fast(nn_ptr ys, const nmod_poly_t poly, ulong r, slong olen)
- Evaluates ``poly`` at the first ``n`` powers
- of the square of ``r``, writing the output values to ``ys``.
+ Evaluates ``poly`` at the first ``olen`` powers
+ of the square of ``r``, writing the output values to ``ys``.
The value of ``r`` should be reduced modulo the modulus of the polynomial
- and of sufficient multiplicative order such that none of
- the first `n` powers of `r^2` is one.
+ and of sufficient multiplicative order such that none of
+ the first ``olen`` powers of `r^2` is one.
Uses fast geometric multipoint evaluation, building a temporary geometric progression precomputation.
@@ -1506,21 +1506,21 @@ Interpolation
--------------------------------------------------------------------------------
-.. function:: void _nmod_poly_interpolate_nmod_vec(nn_ptr poly, nn_srcptr xs, nn_srcptr ys, slong n, nmod_t mod)
+.. function:: void _nmod_poly_interpolate_nmod_vec(nn_ptr poly, nn_srcptr xs, nn_srcptr ys, slong len, nmod_t mod)
- Sets ``poly`` to the unique polynomial of length at most ``n``
- that interpolates the ``n`` given evaluation points ``xs`` and
+ Sets ``poly`` to the unique polynomial of length at most ``len``
+ that interpolates the ``len`` given evaluation points ``xs`` and
values ``ys``. If the interpolating polynomial is shorter than
- length ``n``, the leading coefficients are set to zero.
+ length ``len``, the leading coefficients are set to zero.
The values in ``xs`` and ``ys`` should be reduced modulo the
modulus, and all ``xs`` must be distinct. Aliasing between
``poly`` and ``xs`` or ``ys`` is not allowed.
-.. function:: void nmod_poly_interpolate_nmod_vec(nmod_poly_t poly, nn_srcptr xs, nn_srcptr ys, slong n)
+.. function:: void nmod_poly_interpolate_nmod_vec(nmod_poly_t poly, nn_srcptr xs, nn_srcptr ys, slong len)
- Sets ``poly`` to the unique polynomial of length ``n`` that
- interpolates the ``n`` given evaluation points ``xs`` and
+ Sets ``poly`` to the unique polynomial of length ``len`` that
+ interpolates the ``len`` given evaluation points ``xs`` and
values ``ys``. The values in ``xs`` and ``ys`` should be
reduced modulo the modulus, and all ``xs`` must be distinct.
@@ -1539,34 +1539,34 @@ Interpolation
interpolation weights ``weights`` corresponding to the
roots.
-.. function:: void _nmod_poly_interpolate_nmod_vec_fast(nn_ptr poly, nn_srcptr xs, nn_srcptr ys, slong n, nmod_t mod)
+.. function:: void _nmod_poly_interpolate_nmod_vec_fast(nn_ptr poly, nn_srcptr xs, nn_srcptr ys, slong len, nmod_t mod)
Performs interpolation using the fast Lagrange interpolation
algorithm, generating a temporary subproduct tree.
-.. function:: void nmod_poly_interpolate_nmod_vec_fast(nmod_poly_t poly, nn_srcptr xs, nn_srcptr ys, slong n)
+.. function:: void nmod_poly_interpolate_nmod_vec_fast(nmod_poly_t poly, nn_srcptr xs, nn_srcptr ys, slong len)
Performs interpolation using the fast Lagrange interpolation algorithm,
generating a temporary subproduct tree.
-.. function:: void _nmod_poly_interpolate_nmod_vec_newton(nn_ptr poly, nn_srcptr xs, nn_srcptr ys, slong n, nmod_t mod)
+.. function:: void _nmod_poly_interpolate_nmod_vec_newton(nn_ptr poly, nn_srcptr xs, nn_srcptr ys, slong len, nmod_t mod)
Forms the interpolating polynomial in the Newton basis using
the method of divided differences and then converts it to
monomial form.
-.. function:: void nmod_poly_interpolate_nmod_vec_newton(nmod_poly_t poly, nn_srcptr xs, nn_srcptr ys, slong n)
+.. function:: void nmod_poly_interpolate_nmod_vec_newton(nmod_poly_t poly, nn_srcptr xs, nn_srcptr ys, slong len)
Forms the interpolating polynomial in the Newton basis using
the method of divided differences and then converts it to
monomial form.
-.. function:: void _nmod_poly_interpolate_nmod_vec_barycentric(nn_ptr poly, nn_srcptr xs, nn_srcptr ys, slong n, nmod_t mod)
+.. function:: void _nmod_poly_interpolate_nmod_vec_barycentric(nn_ptr poly, nn_srcptr xs, nn_srcptr ys, slong len, nmod_t mod)
Forms the interpolating polynomial using a naive implementation
of the barycentric form of Lagrange interpolation.
-.. function:: void nmod_poly_interpolate_nmod_vec_barycentric(nmod_poly_t poly, nn_srcptr xs, nn_srcptr ys, slong n)
+.. function:: void nmod_poly_interpolate_nmod_vec_barycentric(nmod_poly_t poly, nn_srcptr xs, nn_srcptr ys, slong len)
Forms the interpolating polynomial using a naive implementation
of the barycentric form of Lagrange interpolation.
@@ -1576,26 +1576,26 @@ Interpolation
Performs interpolation using the geometric progression precomputation ``G``.
- Sets ``poly`` to the unique polynomial of length at most ``len``
- that interpolates according to the parameter set of ``G``.
- The value of ``len`` should be equal to the precomputation size parameter ``G->len``.
+ Sets ``poly`` to the unique polynomial of length at most ``len`` that
+ interpolates the ``len`` values in ``v`` according to the parameter set of
+ ``G``. Requires ``len <= G->len``.
Uses fast geometric multipoint interpolation using a supplied geometric progression precomputation.
-.. function:: void nmod_poly_interpolate_geometric_nmod_vec_fast(nmod_poly_t poly, ulong r, nn_srcptr ys, slong n)
+.. function:: void nmod_poly_interpolate_geometric_nmod_vec_fast(nmod_poly_t poly, ulong r, nn_srcptr ys, slong len)
- Sets ``poly`` to the unique polynomial of length at most ``n``
- that interpolates the first ``n`` powers of ``r`` and
- values ``ys``.
+ Sets ``poly`` to the unique polynomial of length at most ``len``
+ that interpolates the first ``len`` powers of the square of ``r`` and
+ the ``len`` values in ``ys``.
The values ``ys`` and ``r`` should be reduced modulo the
- modulus, and all ``r`` should be of sufficient order such that
- none of the first `n` powers of `r^2` is one. Aliasing between
+ modulus, and ``r`` should be of sufficient order such that
+ none of the first ``len`` powers of `r^2` is one. Aliasing between
``poly`` and ``ys`` is not allowed.
Uses fast geometric multipoint interpolation, building a temporary geometric progression precomputation.
-
+
Composition
--------------------------------------------------------------------------------
@@ -2586,10 +2586,17 @@ Geometric progression
--------------------------------------------------------------------------------
-.. function:: void nmod_geometric_progression_init(nmod_geometric_progression_t G, ulong r, slong len, nmod_t mod)
+.. function:: void _nmod_geometric_progression_init_function(nmod_geometric_progression_t G, ulong r, slong len, nmod_t mod, ulong function)
+ void nmod_geometric_progression_init(nmod_geometric_progression_t G, ulong r, slong len, nmod_t mod)
Builds a geometric progression multipoint evaluation / interpolation structure.
+ The variant with ``function`` variant builds precomputation for specific
+ functionalities: currently, one should set ``function`` to `1` for
+ evaluation only, to `2` for interpolation only, and to `3` for both
+ evaluation and interpolation. The variant without ``function`` precomputes
+ for both.
+
The set of points used will be `1, r^2, r^4, \ldots, r^{2(len-1)}`.
The value of ``r`` should be reduced modulo the modulus ``mod``
@@ -2598,10 +2605,10 @@ Geometric progression
The value of ``len`` should be both greater than or equal to the number of evaluation points to be
considered, and greater than or equal to the length of the polynomials to be evaluated / interpolated.
- This allocates vectors and polynomials for a total space of `8 len - 1` coefficients.
+ This allocates vectors and polynomials for a total space of `6 len - 1` coefficients.
If the modulus is not prime, this function will work under the additional
- assumption that all the used points `r^{2k}` as well as the axuiliary
+ assumption that all the used points `r^{2k}` as well as the auxiliary
values `r^{2k} - 1` are invertible.
.. function:: void nmod_geometric_progression_clear(nmod_geometric_progression_t G)
diff --git a/src/nmod_poly.h b/src/nmod_poly.h
index cfe9154b5e..29b0e74801 100644
--- a/src/nmod_poly.h
+++ b/src/nmod_poly.h
@@ -60,11 +60,13 @@ nmod_poly_compose_mod_precomp_preinv_arg_t;
typedef struct
{
- nn_ptr x, t, w, y, z; // five vectors of precomputed constants
- nmod_poly_t f, g1, g2; // three precomputed polys
- nmod_t mod;
- slong len; // number of points
-
+ nn_ptr ev_s; /* evaluate: scaling constants */
+ nmod_poly_t ev_f; /* evaluate: polynomial */
+ nn_ptr int_s1, int_s2; /* interpolate: scaling constants */
+ nmod_poly_t int_f; /* interpolate: polynomial */
+ nmod_t mod; /* modulus */
+ slong len; /* number of points */
+ ulong function; /* choice of precomputations */
} nmod_geometric_progression_struct;
typedef nmod_geometric_progression_struct nmod_geometric_progression_t[1];
@@ -463,20 +465,20 @@ void _nmod_poly_powmod_x_fmpz_preinv (nn_ptr res, fmpz_t e, nn_srcptr f, slong l
void nmod_poly_powmod_x_fmpz_preinv(nmod_poly_t res, fmpz_t e, const nmod_poly_t f, const nmod_poly_t finv);
void _nmod_poly_powers_mod_preinv_naive(nn_ptr * res, nn_srcptr f,
- slong flen, slong n, nn_srcptr g, slong glen, nn_srcptr ginv,
- slong ginvlen, const nmod_t mod);
+ slong flen, slong n, nn_srcptr g, slong glen, nn_srcptr ginv,
+ slong ginvlen, const nmod_t mod);
void nmod_poly_powers_mod_naive(nmod_poly_struct * res,
const nmod_poly_t f, slong n, const nmod_poly_t g);
void _nmod_poly_powers_mod_preinv_threaded_pool(nn_ptr * res,
- nn_srcptr f, slong flen, slong n, nn_srcptr g, slong glen,
- nn_srcptr ginv, slong ginvlen, const nmod_t mod,
- thread_pool_handle * threads, slong num_threads);
+ nn_srcptr f, slong flen, slong n, nn_srcptr g, slong glen,
+ nn_srcptr ginv, slong ginvlen, const nmod_t mod,
+ thread_pool_handle * threads, slong num_threads);
void
_nmod_poly_powers_mod_preinv_threaded(nn_ptr * res, nn_srcptr f,
- slong flen, slong n, nn_srcptr g, slong glen,
+ slong flen, slong n, nn_srcptr g, slong glen,
nn_srcptr ginv, slong ginvlen, const nmod_t mod);
void nmod_poly_powers_mod_bsgs(nmod_poly_struct * res,
@@ -552,16 +554,16 @@ ulong _nmod_poly_evaluate_nmod_precomp(nn_srcptr poly, slong len, ulong c, ulong
ulong _nmod_poly_evaluate_nmod_precomp_lazy(nn_srcptr poly, slong len, ulong c, ulong c_precomp, ulong modn);
ulong nmod_poly_evaluate_nmod(const nmod_poly_t poly, ulong c);
-void _nmod_poly_evaluate_nmod_vec(nn_ptr ys, nn_srcptr coeffs, slong len, nn_srcptr xs, slong n, nmod_t mod);
-void nmod_poly_evaluate_nmod_vec(nn_ptr ys, const nmod_poly_t poly, nn_srcptr xs, slong n);
+void _nmod_poly_evaluate_nmod_vec(nn_ptr ys, nn_srcptr coeffs, slong ilen, nn_srcptr xs, slong olen, nmod_t mod);
+void nmod_poly_evaluate_nmod_vec(nn_ptr ys, const nmod_poly_t poly, nn_srcptr xs, slong olen);
-void _nmod_poly_evaluate_nmod_vec_iter(nn_ptr ys, nn_srcptr coeffs, slong len, nn_srcptr xs, slong n, nmod_t mod);
-void nmod_poly_evaluate_nmod_vec_iter(nn_ptr ys, const nmod_poly_t poly, nn_srcptr xs, slong n);
+void _nmod_poly_evaluate_nmod_vec_iter(nn_ptr ys, nn_srcptr coeffs, slong ilen, nn_srcptr xs, slong olen, nmod_t mod);
+void nmod_poly_evaluate_nmod_vec_iter(nn_ptr ys, const nmod_poly_t poly, nn_srcptr xs, slong olen);
-void _nmod_poly_evaluate_nmod_vec_fast_precomp(nn_ptr vs, nn_srcptr poly, slong plen, const nn_ptr * tree, slong len, nmod_t mod);
+void _nmod_poly_evaluate_nmod_vec_fast_precomp(nn_ptr vs, nn_srcptr poly, slong ilen, const nn_ptr * tree, slong olen, nmod_t mod);
-void _nmod_poly_evaluate_nmod_vec_fast(nn_ptr ys, nn_srcptr coeffs, slong len, nn_srcptr xs, slong n, nmod_t mod);
-void nmod_poly_evaluate_nmod_vec_fast(nn_ptr ys, const nmod_poly_t poly, nn_srcptr xs, slong n);
+void _nmod_poly_evaluate_nmod_vec_fast(nn_ptr ys, nn_srcptr coeffs, slong ilen, nn_srcptr xs, slong olen, nmod_t mod);
+void nmod_poly_evaluate_nmod_vec_fast(nn_ptr ys, const nmod_poly_t poly, nn_srcptr xs, slong olen);
void nmod_mat_one_addmul(nmod_mat_t dest, const nmod_mat_t mat, ulong c);
@@ -592,34 +594,38 @@ void _nmod_poly_tree_build(nn_ptr * tree, nn_srcptr roots, slong len, nmod_t mod
/* Geometric evaluation / interpolation *************************************/
-void nmod_geometric_progression_init(nmod_geometric_progression_t G, ulong r, slong len, nmod_t mod);
+void _nmod_geometric_progression_init_function(nmod_geometric_progression_t G,
+ ulong r, slong len, nmod_t mod,
+ ulong function);
+void nmod_geometric_progression_init(nmod_geometric_progression_t G,
+ ulong r, slong len, nmod_t mod);
void nmod_geometric_progression_clear(nmod_geometric_progression_t G);
-void _nmod_poly_evaluate_geometric_nmod_vec_iter(nn_ptr ys, nn_srcptr coeffs, slong len, ulong r, slong n, nmod_t mod);
-void nmod_poly_evaluate_geometric_nmod_vec_iter(nn_ptr ys, const nmod_poly_t poly, ulong r, slong n);
+void _nmod_poly_evaluate_geometric_nmod_vec_iter(nn_ptr ys, nn_srcptr coeffs, slong ilen, ulong r, slong olen, nmod_t mod);
+void nmod_poly_evaluate_geometric_nmod_vec_iter(nn_ptr ys, const nmod_poly_t poly, ulong r, slong olen);
-void _nmod_poly_evaluate_geometric_nmod_vec_fast_precomp(nn_ptr vs, nn_srcptr poly, slong plen, const nmod_geometric_progression_t G, slong len, nmod_t mod);
-void _nmod_poly_evaluate_geometric_nmod_vec_fast(nn_ptr ys, nn_srcptr coeffs, slong len, ulong r, slong n, nmod_t mod);
-void nmod_poly_evaluate_geometric_nmod_vec_fast(nn_ptr ys, const nmod_poly_t poly, ulong r, slong n);
+void _nmod_poly_evaluate_geometric_nmod_vec_fast_precomp(nn_ptr vs, nn_srcptr poly, slong ilen, const nmod_geometric_progression_t G, slong olen, nmod_t mod);
+void _nmod_poly_evaluate_geometric_nmod_vec_fast(nn_ptr ys, nn_srcptr coeffs, slong ilen, ulong r, slong olen, nmod_t mod);
+void nmod_poly_evaluate_geometric_nmod_vec_fast(nn_ptr ys, const nmod_poly_t poly, ulong r, slong olen);
void _nmod_poly_interpolate_geometric_nmod_vec_fast_precomp(nn_ptr poly, nn_srcptr v, const nmod_geometric_progression_t G, slong len, nmod_t mod);
void nmod_poly_interpolate_geometric_nmod_vec_fast_precomp(nmod_poly_t poly, nn_srcptr v, const nmod_geometric_progression_t G, slong len);
-void nmod_poly_interpolate_geometric_nmod_vec_fast(nmod_poly_t poly, ulong r, nn_srcptr ys, slong n);
+void nmod_poly_interpolate_geometric_nmod_vec_fast(nmod_poly_t poly, ulong r, nn_srcptr ys, slong len);
/* Interpolation ************************************************************/
-void _nmod_poly_interpolate_nmod_vec_newton(nn_ptr poly, nn_srcptr xs, nn_srcptr ys, slong n, nmod_t mod);
+void _nmod_poly_interpolate_nmod_vec_newton(nn_ptr poly, nn_srcptr xs, nn_srcptr ys, slong len, nmod_t mod);
void nmod_poly_interpolate_nmod_vec_newton(nmod_poly_t poly, nn_srcptr xs, nn_srcptr ys, slong n);
-void _nmod_poly_interpolate_nmod_vec_barycentric(nn_ptr poly, nn_srcptr xs, nn_srcptr ys, slong n, nmod_t mod);
-void nmod_poly_interpolate_nmod_vec_barycentric(nmod_poly_t poly, nn_srcptr xs, nn_srcptr ys, slong n);
+void _nmod_poly_interpolate_nmod_vec_barycentric(nn_ptr poly, nn_srcptr xs, nn_srcptr ys, slong len, nmod_t mod);
+void nmod_poly_interpolate_nmod_vec_barycentric(nmod_poly_t poly, nn_srcptr xs, nn_srcptr ys, slong len);
-void _nmod_poly_interpolate_nmod_vec(nn_ptr poly, nn_srcptr xs, nn_srcptr ys, slong n, nmod_t mod);
-void nmod_poly_interpolate_nmod_vec(nmod_poly_t poly, nn_srcptr xs, nn_srcptr ys, slong n);
+void _nmod_poly_interpolate_nmod_vec(nn_ptr poly, nn_srcptr xs, nn_srcptr ys, slong len, nmod_t mod);
+void nmod_poly_interpolate_nmod_vec(nmod_poly_t poly, nn_srcptr xs, nn_srcptr ys, slong len);
void _nmod_poly_interpolate_nmod_vec_fast(nn_ptr poly, nn_srcptr xs, nn_srcptr ys, slong len, nmod_t mod);
-void nmod_poly_interpolate_nmod_vec_fast(nmod_poly_t poly, nn_srcptr xs, nn_srcptr ys, slong n);
+void nmod_poly_interpolate_nmod_vec_fast(nmod_poly_t poly, nn_srcptr xs, nn_srcptr ys, slong len);
void _nmod_poly_interpolate_nmod_vec_fast_precomp(nn_ptr poly, nn_srcptr ys,
const nn_ptr * tree, nn_srcptr weights, slong len, nmod_t mod);
@@ -695,7 +701,7 @@ _nmod_poly_compose_mod_brent_kung_vec_preinv(nmod_poly_struct * res,
void nmod_poly_compose_mod_brent_kung_vec_preinv(nmod_poly_struct * res,
const nmod_poly_struct * polys, slong len1, slong n,
const nmod_poly_t g, const nmod_poly_t poly,
- const nmod_poly_t polyinv);
+ const nmod_poly_t polyinv);
void _nmod_poly_compose_mod_brent_kung_vec_preinv_worker(void * arg_ptr);
diff --git a/src/nmod_poly/evaluate_geometric_nmod_vec.c b/src/nmod_poly/evaluate_geometric_nmod_vec.c
index c9117a2a26..b6738c74a0 100644
--- a/src/nmod_poly/evaluate_geometric_nmod_vec.c
+++ b/src/nmod_poly/evaluate_geometric_nmod_vec.c
@@ -21,6 +21,7 @@ void _nmod_poly_evaluate_geometric_nmod_vec_fast_precomp(nn_ptr vs, nn_srcptr po
const nmod_geometric_progression_t G, slong len,
nmod_t mod)
{
+ FLINT_ASSERT(G->function & 1);
FLINT_ASSERT(len <= G->len);
if (plen == 0)
@@ -39,79 +40,63 @@ void _nmod_poly_evaluate_geometric_nmod_vec_fast_precomp(nn_ptr vs, nn_srcptr po
}
}
- /* below are 3 different versions (one requires fft_small) */
- /* TODO once some optimized middle product is written, the "version 1" should probably be discarded */
- /** goal is to compute [rev(p) * G->f]_{plen - 1}^{len} (that is, coeffs [plen - 1, plen - 1 + len))
- * if p has valuation val, define a = p / x**val of length alen = plen - val, and this becomes
- * [rev(a) * (G->f >> x**val)]_{alen - 1}^{len} (that is, coeffs [alen - 1, alen - 1 + len))
- **/
+ /** Formula, based on Bluestein's trick:
+ * poly(q**j) = sum_{i=0}^{n-1} poly[i] q**{i*j} is the coefficient x**{n+j-1} of the product
+ * (sum_{i=0}^{len-1} poly[i] * q**(-i*i / 2) x**{len-i-1}) * (sum_{i=0}^{2*len-2} q**{i*i/2} x**i)
+ * -> the right-hand polynomial is G->ev_f truncated at precision 2*len - 1
+ * -> the inverses q**(-i*i / 2) are the first len values stored in G->ev_s
+ *
+ * Thus, goal is to compute [rev(p) * G->ev_f]_{plen - 1}^{len}, where p is poly scaled by G->ev_s
+ * (the bracket notation means coefficients [plen - 1, plen - 1 + len))
+ * If p has valuation val, define a = p / x**val of length alen = plen - val, and this becomes
+ * [rev(a) * (G->ev_f >> x**val)]_{alen - 1}^{len} (that is, coeffs [alen - 1, alen - 1 + len))
+ */
+
+ /* below are 2 different versions (one requiring fft_small) */
+ /* TODO some optimization needed in mulmid: */
+ /* mulmid is often excellent, but fft_small variant still useful */
const slong alen = plen - val;
- /* version 1: (2025-12-04: fastest in small lengths, waiting for optimized middle product) */
- /** this uses a short product: write rev(a) = x**(alen-1) a(1/x), and write F = (G->f(x) >> x**val) rem x**(alen - 1 + len)
- * rev(a) * F has length L = 2 * alen - 2 + len, reverse it: we get a * rev(F),
- * we want its coefficients from L - 1 - (alen - 1) = alen - 1 + len - 1
- * down to, included, L - 1 - (alen - 1 + len - 1) = alen - 1
- **/
+ /* this uses a middle product to compute [rev(p) * G->ev_f]_{plen - 1}^{len} (i.e. coeffs [plen - 1, plen - 1 + len)) */
#if FLINT_HAVE_FFT_SMALL
- if (2 * (plen - val) - 2 + len <= 192)
-#else
- if (1)
-#endif
+ if (2 * (plen - val) - 2 + len > 192)
{
- slong blen = alen - 1 + len;
- nn_ptr a = _nmod_vec_init(alen);
- nn_ptr b = _nmod_vec_init(blen);
+ /* uses fft_small directly */
+ /* 2025-12-04: fastest in medium and large lengths, like 100 and more */
+ /* 2026-05-03: still useful for some short range where mulmid is not yet good */
+ nn_ptr b = _nmod_vec_init(alen + len - 1);
for (slong i = val; i < plen; i++)
- a[i - val] = nmod_mul(G->x[i], poly[i], mod);
+ b[plen - 1 - i] = nmod_mul(G->ev_s[i], poly[i], mod);
- nn_ptr Frev = _nmod_vec_init(blen);
- _nmod_poly_reverse(Frev, G->f->coeffs + val, blen, blen);
- _nmod_poly_mullow(b, Frev, blen, a, alen, blen, mod);
+ _nmod_poly_mul_mid_default_mpn_ctx(b, alen - 1, alen - 1 + len, G->ev_f->coeffs + val, alen - 1 + len, b, alen, mod);
for (slong i = 0; i < len; i++)
- vs[i] = nmod_mul(G->x[i], b[alen - 1 + len - 1 - i], mod);
+ vs[i] = nmod_mul(G->ev_s[i], b[i], mod);
- _nmod_vec_clear(Frev);
- _nmod_vec_clear(a);
_nmod_vec_clear(b);
}
else
+#endif
{
- /* version 2 */
- /* this uses a middle product to compute [rev(p) * G->f]_{plen - 1}^{len} (i.e. coeffs [plen - 1, plen - 1 + len)) */
-#if FLINT_HAVE_FFT_SMALL
- /* version 2.a uses fft_small directly (2025-12-04: fastest in medium and large lengths, like 100 and more) */
- nn_ptr b = _nmod_vec_init(alen + len - 1);
-
- for (slong i = val; i < plen; i++)
- b[plen - 1 - i] = nmod_mul(G->x[i], poly[i], mod);
-
- _nmod_poly_mul_mid_default_mpn_ctx(b, alen - 1, alen - 1 + len, G->f->coeffs + val, alen - 1 + len, b, alen, mod);
-
- for (slong i = 0; i < len; i++)
- vs[i] = nmod_mul(G->x[i], b[i], mod);
-
- _nmod_vec_clear(b);
-#else
- /* version 2.b uses nmod_poly_mulhigh (2025-12-04: tested correct, but disabled: nmod_poly_mulhigh not yet optimized) */
- /* (currently, with the branching mechanism above, this should never be called) */
+ /* uses nmod_poly_mulmid */
+ /* 2025-12-04: disabled: nmod_poly_mulhigh/mulmid not yet optimized */
+ /* 2026-05-03: best method, much better than fft_small for small parameters, */
+ /* but also sometimes quite slower for some parameter ranges */
nn_ptr a = _nmod_vec_init(alen);
- nn_ptr b = _nmod_vec_init(alen + (alen - 1 + len));
+ nn_ptr b = _nmod_vec_init(alen - 1 + len);
for (slong i = val; i < plen; i++)
- a[plen - 1 - i] = nmod_mul(G->x[i], poly[i], mod);
+ a[plen - 1 - i] = nmod_mul(G->ev_s[i], poly[i], mod);
- _nmod_poly_mulhigh(b, G->f->coeffs + val, alen - 1 + len, a, alen, alen - 1, mod);
+ _nmod_poly_mulmid(b, G->ev_f->coeffs + val, alen - 1 + len, a, alen, alen - 1, alen - 1 + len, mod);
for (slong i = 0; i < len; i++)
- vs[i] = nmod_mul(G->x[i], b[alen - 1 + i], mod);
+ vs[i] = nmod_mul(G->ev_s[i], b[i], mod);
_nmod_vec_clear(a);
_nmod_vec_clear(b);
-#endif
}
}
@@ -121,7 +106,7 @@ void _nmod_poly_evaluate_geometric_nmod_vec_fast(nn_ptr ys, nn_srcptr poly, slon
return;
nmod_geometric_progression_t G;
- nmod_geometric_progression_init(G, r, FLINT_MAX(n, plen), mod);
+ _nmod_geometric_progression_init_function(G, r, FLINT_MAX(n, plen), mod, UWORD(1));
_nmod_poly_evaluate_geometric_nmod_vec_fast_precomp(ys, poly, plen, G, n, mod);
nmod_geometric_progression_clear(G);
}
diff --git a/src/nmod_poly/geometric_progression.c b/src/nmod_poly/geometric_progression.c
index 3537a2c906..139fcde853 100644
--- a/src/nmod_poly/geometric_progression.c
+++ b/src/nmod_poly/geometric_progression.c
@@ -1,6 +1,6 @@
/*
- Copyright (C) 2025, Vincent Neiger, Éric Schost
- Copyright (C) 2025, Mael Hostettler
+ Copyright (C) 2025, Vincent Neiger, Éric Schost, Mael Hostettler
+ Copyright (C) 2026, Vincent Neiger
This file is part of FLINT.
@@ -10,127 +10,253 @@
(at your option) any later version. See .
*/
+#include "nmod.h"
#include "nmod_vec.h"
#include "nmod_poly.h"
-void
-nmod_geometric_progression_init(nmod_geometric_progression_t G, ulong r, slong len, nmod_t mod)
+/* specialized for moduli that support n_mulmod_shoup */
+static
+void _nmod_geometric_progression_evaluate_init_nonfullword(nmod_geometric_progression_t G,
+ ulong r, slong len, nmod_t mod,
+ ulong q, ulong q_pr_quo, ulong q_pr_rem,
+ ulong inv_r,
+ ulong inv_q, ulong inv_q_pr_quo, ulong inv_q_pr_rem)
{
- ulong q, inv_r, inv_q, tmp, qk, inv_qk, qq, s;
- nn_ptr diff, inv_diff, prod_diff;
- slong i;
-
- G->len = len;
- G->mod = mod;
+ /* G->ev_f = sum_{0 <= i < 2*len - 1} q**(i*i/2) * x**i */
+ /* G->ev_s[i] = 1 / q**(i*i/2) = 1 / r**(i*i) */
+ nmod_poly_init2_preinv(G->ev_f, mod.n, mod.ninv, 2*len - 1);
+ G->ev_s = _nmod_vec_init(len);
+
+ G->ev_f->length = 2*len - 1;
+
+ /* precomputations for Shoup multiplication */
+ ulong r_pow_i2 = r;
+ ulong r_pr_quo = n_mulmod_precomp_shoup(r_pow_i2, mod.n);
- nmod_poly_init2(G->f, mod.n, 2*len - 1); G->f->length = 2*len - 1;
- nmod_poly_init2(G->g1, mod.n, len); G->g1->length = len;
- nmod_poly_init2(G->g2, mod.n, len); G->g2->length = len;
- G->g1->coeffs[0] = 1;
- G->g2->coeffs[0] = 1;
- G->f->coeffs[0] = 1;
-
- G->x = _nmod_vec_init(len);
- G->w = _nmod_vec_init(len);
- G->z = _nmod_vec_init(len);
- G->y = _nmod_vec_init(len);
-
- G->x[0] = 1;
- G->y[0] = 1;
- G->w[0] = 1;
- G->z[0] = 1;
-
- q = nmod_mul(r, r, mod);
- inv_r = nmod_inv(r, mod);
- inv_q = nmod_mul(inv_r, inv_r, mod);
-
- tmp = r;
- for (i = 1; i < 2*len - 1; i++)
+ G->ev_f->coeffs[0] = 1;
+ for (slong i = 1; i < 2*len - 1; i++)
{
- G->f->coeffs[i] = nmod_mul(G->f->coeffs[i - 1], tmp, mod);
- tmp = nmod_mul(tmp, q, mod);
+ G->ev_f->coeffs[i] = n_mulmod_shoup(r_pow_i2, G->ev_f->coeffs[i - 1], r_pr_quo, mod.n);
+ n_mulmod_and_precomp_shoup(&r_pow_i2, &r_pr_quo, q, r_pow_i2, q_pr_quo, q_pr_rem, r_pr_quo, mod.n);
}
- // If we had to normalize G->f than that means r is of low order and following
- // inversion will fail
- tmp = inv_r;
- for (i = 1; i < len; i++)
+ r_pr_quo = n_mulmod_precomp_shoup(inv_r, mod.n);
+ G->ev_s[0] = 1;
+ for (slong i = 1; i < len; i++)
{
- G->x[i] = nmod_mul(G->x[i - 1], tmp, mod);
- tmp = nmod_mul(tmp, inv_q, mod);
+ G->ev_s[i] = n_mulmod_shoup(inv_r, G->ev_s[i - 1], r_pr_quo, mod.n);
+ n_mulmod_and_precomp_shoup(&inv_r, &r_pr_quo, inv_q, inv_r, inv_q_pr_quo, inv_q_pr_rem, r_pr_quo, mod.n);
}
-
- inv_diff = _nmod_vec_init(len);
- diff = _nmod_vec_init(len);
- prod_diff = _nmod_vec_init(len);
- inv_diff[0] = 1;
- diff[0] = 1;
- prod_diff[0] = 1;
-
- qk = q; // montgomery inversion
- for (i = 1; i < len; i++)
+}
+
+/* general variant */
+static
+void _nmod_geometric_progression_evaluate_init(nmod_geometric_progression_t G,
+ ulong r, slong len, nmod_t mod,
+ ulong q, ulong inv_r, ulong inv_q)
+{
+ /* G->ev_f = sum_{0 <= i < 2*len - 1} q**(i*i/2) * x**i */
+ /* G->ev_s[i] = 1 / q**(i*i/2) */
+ nmod_poly_init2_preinv(G->ev_f, mod.n, mod.ninv, 2*len - 1);
+ G->ev_s = _nmod_vec_init(len);
+
+ G->ev_f->length = 2*len - 1;
+ G->ev_f->coeffs[0] = 1;
+ for (slong i = 1; i < 2*len - 1; i++)
{
- diff[i] = qk - 1;
- inv_diff[i] = diff[i];
- qk = nmod_mul(qk, q, mod);
- prod_diff[i] = nmod_mul(diff[i], prod_diff[i - 1], mod);
+ G->ev_f->coeffs[i] = nmod_mul(G->ev_f->coeffs[i - 1], r, mod);
+ r = nmod_mul(r, q, mod); /* r**(2*i+1) */
}
-
- tmp = nmod_inv(prod_diff[len-1], mod);
- for (i = len - 1; i > 0; i--)
+
+ G->ev_s[0] = 1;
+ for (slong i = 1; i < len; i++)
{
- inv_diff[i] = nmod_mul(prod_diff[i - 1], tmp, mod);
- tmp = nmod_mul(tmp, diff[i], mod);
+ G->ev_s[i] = nmod_mul(G->ev_s[i - 1], inv_r, mod);
+ inv_r = nmod_mul(inv_r, inv_q, mod);
}
- inv_diff[0] = tmp;
- // end montgomery inversion
+}
+
+static
+void _nmod_geometric_progression_interpolate_init_nonfullword(nmod_geometric_progression_t G,
+ slong len, nmod_t mod,
+ ulong q, ulong q_pr_quo, ulong q_pr_rem,
+ ulong inv_q, ulong inv_q_pr_quo, ulong inv_q_pr_rem)
+{
+ /* quantities for Newton interpolation/evaluation/change-of-basis */
+ /* see [Bostan - Schost, J.Complexity 2005, Section 5.1] */
+ /* write u_i for prod_{1 <= k <= i} (q**k - 1), */
+ /* and q_i for q**(i*(i-1)/2) */
+
+ /* coeff(G->int_f, i) = (-1)**i * q_i / u_i */
+ nmod_poly_init2_preinv(G->int_f, mod.n, mod.ninv, len);
+ G->int_f->length = len;
+
+ /* G->int_s1[i] = 1 / u_i */
+ /* G->int_s2[i] = u_i / q_i */
+ G->int_s1 = _nmod_vec_init(len);
+ G->int_s2 = _nmod_vec_init(len);
- // sets sequences w, y, z and polynomials g1, g2
- qk = 1;
- inv_qk = 1;
- qq = 1;
- s = 1;
+ G->int_f->coeffs[0] = 1;
+ G->int_s2[0] = 1;
- for (i = 1; i < len; i++)
+ const ulong one_precomp = n_mulmod_precomp_shoup(UWORD(1), mod.n);
+ ulong q_pow_i = 1;
+ ulong q_pow_i_pr = one_precomp;
+ ulong inv_q_pow_i = 1;
+ ulong inv_q_pow_i_pr = one_precomp;
+ ulong inv_q_i = 1;
+ ulong inv_q_i_pr = one_precomp;
+ ulong prod_diff = 1;
+
+ for (slong i = 1; i < len; i++)
+ {
+ ulong inv_q_pow_i_pr_rem = n_mulmod_precomp_shoup_rem_from_quo(inv_q_pow_i_pr, mod.n);
+ n_mulmod_and_precomp_shoup(&inv_q_i, &inv_q_i_pr, inv_q_pow_i, inv_q_i, inv_q_pow_i_pr, inv_q_pow_i_pr_rem, inv_q_i_pr, mod.n); /* 1 / q_i */
+ n_mulmod_and_precomp_shoup(&inv_q_pow_i, &inv_q_pow_i_pr, inv_q, inv_q_pow_i, inv_q_pr_quo, inv_q_pr_rem, inv_q_pow_i_pr, mod.n); /* 1 / q**i */
+ G->int_f->coeffs[i] = n_mulmod_shoup(q_pow_i, G->int_f->coeffs[i-1], q_pow_i_pr, mod.n); /* q_i */
+ n_mulmod_and_precomp_shoup(&q_pow_i, &q_pow_i_pr, q, q_pow_i, q_pr_quo, q_pr_rem, q_pow_i_pr, mod.n); /* q**i */
+ G->int_s1[i-1] = q_pow_i - 1; /* temporarily, q**i - 1 */
+ prod_diff = nmod_mul(q_pow_i - 1, prod_diff, mod); /* u_i */
+ G->int_s2[i] = n_mulmod_shoup(inv_q_i, prod_diff, inv_q_i_pr, mod.n); /* u_i / q_i */
+ }
+
+ G->int_s1[len-1] = nmod_inv(prod_diff, mod); /* 1 / u_{len-1} */
+ for (slong i = len - 1; i > 0; i--)
{
- qq = nmod_mul(qq, qk, mod); // prod q^i
- s = nmod_mul(s, inv_qk, mod); // prod 1/q^i
- G->w[i] = nmod_mul(G->w[i - 1], inv_diff[i], mod); // prod 1/(q^i-1)
- tmp = nmod_mul(qq, G->w[i], mod); // prod q^i/(q^i-1)
- G->g2->coeffs[i] = tmp;
-
- if ((i & 1) == 1) /* i is odd */
- {
- G->g1->coeffs[i] = mod.n - tmp;
- G->y[i] = mod.n - prod_diff[i];
- G->z[i] = mod.n - G->w[i];
- }
- else /* i is even */
- {
- G->g1->coeffs[i] = tmp;
- G->y[i] = prod_diff[i];
- G->z[i] = G->w[i];
- }
- G->y[i] = nmod_mul(G->y[i], s, mod);
-
- qk = nmod_mul(qk, q, mod);
- inv_qk = nmod_mul(inv_qk, inv_q, mod);
+ ulong w_i = G->int_s1[i]; /* 1 / u_i */
+ if (i % 2) /* i odd, -q_i / u_i */
+ G->int_f->coeffs[i] = mod.n - nmod_mul(G->int_f->coeffs[i], w_i, mod);
+ else /* i even, q_i / u_i */
+ G->int_f->coeffs[i] = nmod_mul(G->int_f->coeffs[i], w_i, mod);
+ G->int_s1[i-1] = nmod_mul(G->int_s1[i-1], w_i, mod); /* 1 / u_{i-1} */
}
- // similarly, if either g1 or g2 have leading 0 coefficient, something is wrong
+}
+
+static
+void _nmod_geometric_progression_interpolate_init(nmod_geometric_progression_t G,
+ slong len, nmod_t mod,
+ ulong q, ulong inv_q)
+{
+ /* quantities for Newton interpolation/evaluation/change-of-basis */
+ /* see [Bostan - Schost, J.Complexity 2005, Section 5.1] */
+ /* write u_i for prod_{1 <= k <= i} (q**k - 1), */
+ /* and q_i for q**(i*(i-1)/2) */
- _nmod_vec_clear(prod_diff);
- _nmod_vec_clear(inv_diff);
- _nmod_vec_clear(diff);
+ /* coeff(G->int_f, i) = (-1)**i * q_i / u_i */
+ nmod_poly_init2_preinv(G->int_f, mod.n, mod.ninv, len);
+ G->int_f->length = len;
+
+ /* G->int_s1[i] = 1 / u_i */
+ /* G->int_s2[i] = u_i / q_i */
+ G->int_s1 = _nmod_vec_init(len);
+ G->int_s2 = _nmod_vec_init(len);
+
+ G->int_f->coeffs[0] = 1;
+ G->int_s2[0] = 1;
+
+ ulong q_pow_i = 1;
+ ulong inv_q_pow_i = 1;
+ ulong inv_q_i = 1;
+ ulong prod_diff = 1;
+
+ for (slong i = 1; i < len; i++)
+ {
+ inv_q_i = nmod_mul(inv_q_i, inv_q_pow_i, mod); /* 1 / q_i */
+ inv_q_pow_i = nmod_mul(inv_q_pow_i, inv_q, mod); /* 1 / q**i */
+ G->int_f->coeffs[i] = nmod_mul(G->int_f->coeffs[i-1], q_pow_i, mod); /* q_i */
+ q_pow_i = nmod_mul(q_pow_i, q, mod); /* q**i */
+ G->int_s1[i-1] = q_pow_i - 1; /* temporarily, q**i - 1 */
+ prod_diff = nmod_mul(q_pow_i - 1, prod_diff, mod); /* u_i */
+ G->int_s2[i] = nmod_mul(prod_diff, inv_q_i, mod); /* u_i / q_i */
+ }
+
+ G->int_s1[len-1] = nmod_inv(prod_diff, mod); /* 1 / u_{len-1} */
+ for (slong i = len - 1; i > 0; i--)
+ {
+ ulong w_i = G->int_s1[i]; /* 1 / u_i */
+ if (i % 2) /* i odd, - q_i / u_i */
+ G->int_f->coeffs[i] = mod.n - nmod_mul(G->int_f->coeffs[i], w_i, mod);
+ else /* i even, q_i / u_i */
+ G->int_f->coeffs[i] = nmod_mul(G->int_f->coeffs[i], w_i, mod);
+ G->int_s1[i-1] = nmod_mul(G->int_s1[i-1], w_i, mod); /* 1 / u_{i-1} */
+ }
+}
+
+/* initialize for selection of functionalities: */
+/* the lowest 3 bits of `function` act as a mask for the */
+/* three functionalities, in this order: */
+/* evaluate(bit0)+interpolate(bit1)+extrapolate(bit2) */
+void _nmod_geometric_progression_init_function(nmod_geometric_progression_t G,
+ ulong r, slong len, nmod_t mod,
+ ulong function)
+{
+ G->len = len;
+ G->mod = mod;
+ G->function = function;
+
+ if (NMOD_CAN_USE_SHOUP(mod))
+ {
+ const ulong q = nmod_mul(r, r, mod);
+ const ulong inv_r = nmod_inv(r, mod);
+ const ulong inv_q = nmod_mul(inv_r, inv_r, mod);
+
+ ulong q_pr_rem;
+ ulong q_pr_quo;
+ ulong inv_q_pr_rem;
+ ulong inv_q_pr_quo;
+ n_mulmod_precomp_shoup_quo_rem(&q_pr_quo, &q_pr_rem, q, mod.n);
+ n_mulmod_precomp_shoup_quo_rem(&inv_q_pr_quo, &inv_q_pr_rem, inv_q, mod.n);
+
+ if (function & UWORD(1)) /* evaluate */
+ _nmod_geometric_progression_evaluate_init_nonfullword(G, r, len, mod, q, q_pr_quo, q_pr_rem, inv_r, inv_q, inv_q_pr_quo, inv_q_pr_rem);
+ if ((function>>1) & UWORD(1)) /* interpolate */
+ _nmod_geometric_progression_interpolate_init_nonfullword(G, len, mod, q, q_pr_quo, q_pr_rem, inv_q, inv_q_pr_quo, inv_q_pr_rem);
+ /* if ((function>>2) & UWORD(1)) */
+ /* extrapolate */
+ }
+ else
+ {
+ const ulong q = nmod_mul(r, r, mod);
+ const ulong inv_r = nmod_inv(r, mod);
+ const ulong inv_q = nmod_mul(inv_r, inv_r, mod);
+
+ if (function & UWORD(1)) /* evaluate */
+ _nmod_geometric_progression_evaluate_init(G, r, len, mod, q, inv_r, inv_q);
+ if ((function>>1) & UWORD(1)) /* interpolate */
+ _nmod_geometric_progression_interpolate_init(G, len, mod, q, inv_q);
+ /* if ((function>>2) & UWORD(1)) */
+ /* extrapolate */
+ }
+}
+
+/* clear for selection of functionalities (see init) */
+static
+void _nmod_geometric_progression_clear_function(nmod_geometric_progression_t G, ulong function)
+{
+ if (function & UWORD(1)) /* evaluate */
+ {
+ _nmod_vec_clear(G->ev_s);
+ nmod_poly_clear(G->ev_f);
+ }
+ if ((function>>1) & UWORD(1)) /* interpolate */
+ {
+ _nmod_vec_clear(G->int_s1);
+ _nmod_vec_clear(G->int_s2);
+ nmod_poly_clear(G->int_f);
+ }
+ /* if ((function>>2) & UWORD(1)) */
+ /* extrapolate */
+}
+
+/* initialize for all: evaluate+interpolate+extrapolate */
+void nmod_geometric_progression_init(nmod_geometric_progression_t G,
+ ulong r, slong len, nmod_t mod)
+{
+ _nmod_geometric_progression_init_function(G, r, len, mod, UWORD(7));
}
-void
-nmod_geometric_progression_clear(nmod_geometric_progression_t G)
+/* clear for all: evaluate+interpolate+extrapolate */
+void nmod_geometric_progression_clear(nmod_geometric_progression_t G)
{
- nmod_poly_clear(G->f);
- nmod_poly_clear(G->g2);
- nmod_poly_clear(G->g1);
- _nmod_vec_clear(G->x);
- _nmod_vec_clear(G->z);
- _nmod_vec_clear(G->y);
- _nmod_vec_clear(G->w);
+ _nmod_geometric_progression_clear_function(G, G->function);
}
diff --git a/src/nmod_poly/interpolate_geometric_nmod_vec.c b/src/nmod_poly/interpolate_geometric_nmod_vec.c
index 4d8a93028f..9f80619894 100644
--- a/src/nmod_poly/interpolate_geometric_nmod_vec.c
+++ b/src/nmod_poly/interpolate_geometric_nmod_vec.c
@@ -14,16 +14,11 @@
#include "nmod_poly.h"
#include "nmod_vec.h"
-void
-_nmod_poly_interpolate_geometric_nmod_vec_fast_precomp(nn_ptr poly, nn_srcptr v,
- const nmod_geometric_progression_t G, slong len, nmod_t mod)
+void _nmod_poly_interpolate_geometric_nmod_vec_fast_precomp(nn_ptr poly,
+ nn_srcptr v, const nmod_geometric_progression_t G, slong len, nmod_t mod)
{
- slong i, N, f1_len, f2_len, h1_len, h2_min;
- nn_ptr f, h;
-
- N = G->len;
-
- FLINT_ASSERT(len == N);
+ FLINT_ASSERT(len <= G->len);
+ FLINT_ASSERT((G->function >> 1) & 1);
if (len == 1)
{
@@ -31,75 +26,124 @@ _nmod_poly_interpolate_geometric_nmod_vec_fast_precomp(nn_ptr poly, nn_srcptr v,
return;
}
- f = _nmod_vec_init(N);
- h = _nmod_vec_init(N);
-
- for (i = 0; i < N; i++)
- {
- if (v[N - i - 1] != 0)
- {
+ /** step1: Newton interpolation
+ * [Bostan - Schost, J.Complexity 2005, Section 5.1]
+ * -> The coefficients of the interpolant, in the Newton basis associated
+ * to the geometric progression 1, q, q**2, q**3, etc., are obtained as
+ * c_0 / q_0, ..., c_{len-1} / q_{len-1}, where c_0, ..., c_{len-1} are
+ * the first coefficients of the product
+ * (sum_{i=0}^{len-1} v[i]/u_i x**i) (sum_{i=0}^{len-1} (-1)**i q_i/u_i x**i)
+ * where v[i] is the element at index `i` in the input values `v`,
+ * and u_i = prod_{1 <= k <= i} (q**k - 1),
+ * and q_i = q**(i*(i-1)/2)
+ * -> With the precomputed data, these are the `len` coefficients of
+ * f * G->int_f mod x**len
+ * where f = sum_{i=0}^{len-1} v[i] * G->int_s1[i] x**i
+ */
+
+ /* val = valuation of output poly in Newton basis */
+ slong val = 0;
+ for (; val < len; val++)
+ if (v[val] != 0)
break;
- }
- }
-
- f1_len = N - i;
- h1_len = FLINT_MIN(G->g1->length + f1_len - 1, N);
-
- for (i = 0; i < f1_len; i++)
+ if (val == len)
{
- f[i] = nmod_mul(v[i], G->w[i], mod);
+ _nmod_vec_zero(poly, len);
+ return;
}
- _nmod_poly_mullow(h, G->g1->coeffs, G->g1->length, f, f1_len, N, mod);
- while (h1_len > 0 && h[h1_len - 1] == 0)
- {
- h1_len--;
- }
+ slong f_len, h_len;
+ nn_ptr f = _nmod_vec_init(len);
+ nn_ptr h = _nmod_vec_init(len);
- for (i = 0; i < h1_len; i++)
- {
- if (h[i] != 0)
- {
+ /* actual length of f */
+ for (f_len = len; f_len > val; f_len--)
+ if (v[f_len-1] != 0)
+ break;
+ f_len = f_len - val;
+
+ /* f = sum_{i=val}^{len-1} v[i] * G->int_s1[i] x**{i-val} */
+ /* == sum_{i=0}^{f_len-1} v[i+val] * G->int_s1[i+val] x**i */
+ for (slong i = 0; i < f_len; i++)
+ f[i] = nmod_mul(v[i+val], G->int_s1[i+val], mod);
+
+ /* h = (x**val * f) * G->int_f mod x**len */
+ /* == x**val (f * G->int_f mod x**(len-val)) */
+ /* note: len - val is <= G->int_f->length, since G->intf_1 has */
+ /* length G->len >= len (all its coefficients are nonzero) */
+ _nmod_poly_mullow(h+val, G->int_f->coeffs, len - val, f, f_len, len - val, mod);
+
+ /* for Newton interpolation, here we should compute h[i] = h[i]/q_i */
+ /* yet this "/q_i" will simplify with another operation just below, */
+ /* so we just leave h as it is */
+
+ /** step2: Newton basis -> monomial basis
+ * [Bostan - Schost, J.Complexity 2005, Section 5.2]
+ * -> Convert h[i]/q_i to monomial basis, through the transposed
+ * x**len-truncated multiplication of two polynomials
+ * sum_{i=0}^{len-1} uu_i x**i
+ * and sum_{i=0}^{len-1} (-1)**i * (h[i]/q_i)*q_i/uu_i x**i)
+ * where q_i = q**(i*(i-1)/2) as above,
+ * and uu_i = prod_{1 <= k <= i} q**(k-1) / (1 - q**k)
+ * == (-1)**i * q_i / u_i, for u_i as above
+ * (in the paper this is prod q**k / (1 - q**k), is this a typo?)
+ * This gives `len` coefficients that must then be scaled
+ * by (-1)**i * uu_i / q_i == 1 / u_i
+ * -> Transposing the truncated product of poly1,poly2 of degree < len,
+ * mullow_t(res, poly1, poly2, len)
+ * simply amounts to a (non-tranposed) mullow and reversals:
+ * mullow(res, poly1, rev(poly2, len), len)
+ * res = rev(res, len)
+ * -> With the precomputed data, the two polynomials have coefficients
+ * uu_i == G->int_f[i] and (-1)**i * h[i]/uu_i == h[i] * G->int_s2[i]
+ * meaning that we want to compute
+ * mullow_t(res, G->int_f, F, len)
+ * where F = sum_{i=0}^{len-1} h[i] * G->int_s2[i] x**i,
+ * and then scale by 1/u_i == G->int_s1[i]
+ */
+
+ /* valuation of h is at least val, see if it is higher */
+ for (; val < len; val++)
+ if (h[val] != 0)
break;
- }
- }
- h2_min = i;
- for (i = h2_min; i < h1_len; i++)
- {
- f[N - 1 - i] = nmod_mul(h[i], G->y[i], mod);
- }
+ /* actual length of h */
+ for (h_len = len; h_len > val; h_len--)
+ if (h[h_len-1] != 0)
+ break;
- f2_len = N - h2_min;
- _nmod_vec_zero(f, N - h1_len);
- _nmod_poly_mullow(h, G->g2->coeffs, G->g2->length, f, f2_len, N, mod);
+ /* compute reversed and scaled f */
+ for (slong i = 0; i < h_len-val; i++)
+ f[i] = nmod_mul(h[h_len-1-i], G->int_s2[h_len-1-i], mod);
- for (i = 0; i < len; i++)
- {
- poly[i] = nmod_mul(h[N - 1 - i], G->z[i], mod);
- }
+ /* transposed short product */
+ _nmod_poly_mullow(h+len-h_len, f, h_len-val, G->int_f->coeffs, h_len, h_len, mod);
+
+ /* final scaling */
+ _nmod_vec_zero(poly+h_len, len-h_len);
+ for (slong i = 0; i < h_len; i++)
+ poly[i] = nmod_mul(h[len-1-i], G->int_s1[i], mod);
_nmod_vec_clear(f);
_nmod_vec_clear(h);
}
-void
-nmod_poly_interpolate_geometric_nmod_vec_fast_precomp(nmod_poly_t poly, nn_srcptr v,
- const nmod_geometric_progression_t G, slong len)
+void nmod_poly_interpolate_geometric_nmod_vec_fast_precomp(nmod_poly_t poly,
+ nn_srcptr v, const nmod_geometric_progression_t G, slong len)
{
+ FLINT_ASSERT((G->function >> 1) & 1);
+
nmod_poly_fit_length(poly, len);
_nmod_poly_set_length(poly, len);
_nmod_poly_interpolate_geometric_nmod_vec_fast_precomp(poly->coeffs, v, G, len, G->mod);
_nmod_poly_normalise(poly);
}
-void
-nmod_poly_interpolate_geometric_nmod_vec_fast(nmod_poly_t poly, ulong r,
- nn_srcptr ys, slong n)
+void nmod_poly_interpolate_geometric_nmod_vec_fast(nmod_poly_t poly,
+ ulong r, nn_srcptr ys, slong len)
{
nmod_geometric_progression_t G;
- nmod_geometric_progression_init(G, r, n, poly->mod);
-
- nmod_poly_interpolate_geometric_nmod_vec_fast_precomp(poly, ys, G, n);
+ _nmod_geometric_progression_init_function(G, r, len, poly->mod, UWORD(2));
+ nmod_poly_interpolate_geometric_nmod_vec_fast_precomp(poly, ys, G, len);
nmod_geometric_progression_clear(G);
}
diff --git a/src/nmod_poly/profile/p-eval_interp_precomputations.c b/src/nmod_poly/profile/p-eval_interp_precomputations.c
new file mode 100644
index 0000000000..d5571f3479
--- /dev/null
+++ b/src/nmod_poly/profile/p-eval_interp_precomputations.c
@@ -0,0 +1,303 @@
+/*
+ Copyright (C) 2026 Vincent Neiger
+
+ This file is part of FLINT.
+
+ FLINT is free software: you can redistribute it and/or modify it under
+ the terms of the GNU Lesser General Public License (LGPL) as published
+ by the Free Software Foundation; either version 3 of the License, or
+ (at your option) any later version. See .
+*/
+
+#include /* for atoi */
+#include
+
+#include "nmod.h"
+#include "profiler.h"
+#include "ulong_extras.h"
+#include "nmod_poly.h"
+#include "nmod_vec.h"
+
+#define __NB_ITER 10
+
+/*------------------------------------------------------------*/
+/* finds an element of order at least n */
+/* returns 0 if not found */
+/*------------------------------------------------------------*/
+static ulong nmod_find_root(slong n, nmod_t mod)
+{
+ ulong attempts = 0;
+ for (ulong q = 2; q < mod.n; q++)
+ {
+ slong k = 1;
+ slong qk = q;
+ while (qk != 1 && k < n)
+ {
+ qk = nmod_mul(qk, q, mod);
+ k++;
+ }
+ if (qk != 1)
+ {
+ return q;
+ }
+ attempts += 1;
+ if (attempts >= 10)
+ return 0;
+ }
+ return 0;
+}
+
+typedef struct
+{
+ flint_bitcnt_t bits;
+ slong length;
+ slong npoints_precomp;
+} info_t;
+
+/* precomputation for general points (subproduct tree) */
+void sample_general_precomp(void * arg, ulong count)
+{
+ ulong n;
+ nmod_t mod;
+ ulong i;
+
+ info_t * info = (info_t *) arg;
+ flint_bitcnt_t bits = info->bits;
+ slong npoints_precomp = info->npoints_precomp;
+
+ FLINT_TEST_INIT(state);
+
+ nn_ptr pts = _nmod_vec_init(npoints_precomp);
+
+ for (i = 0; i < count; i++)
+ {
+ n = n_randbits(state, bits);
+ if (n == UWORD(0)) n++;
+ nmod_init(&mod, n);
+ _nmod_vec_rand(pts, state, npoints_precomp, mod);
+
+ nn_ptr * tree = _nmod_poly_tree_alloc(npoints_precomp);
+
+ prof_start();
+ for (ulong ii = 0; ii < __NB_ITER; ii++)
+ _nmod_poly_tree_build(tree, pts, npoints_precomp, mod);
+ prof_stop();
+
+ _nmod_poly_tree_free(tree, npoints_precomp);
+ }
+
+ _nmod_vec_clear(pts);
+
+ FLINT_TEST_CLEAR(state);
+}
+
+/* precomputations for geometric progression (eval+interp) */
+void sample_geometric_precomp(void * arg, ulong count)
+{
+ ulong n;
+ nmod_t mod;
+ ulong i;
+
+ info_t * info = (info_t *) arg;
+ flint_bitcnt_t bits = info->bits;
+ slong npoints_precomp = info->npoints_precomp;
+
+ FLINT_TEST_INIT(state);
+
+ for (i = 0; i < count; i++)
+ {
+ n = n_randprime(state, bits, 1);
+ nmod_init(&mod, n);
+ ulong r = nmod_find_root(2*npoints_precomp, mod);
+ if (r == 0)
+ flint_printf("\n...could not find element of suitable order for geometric progression...\n");
+
+ nmod_geometric_progression_t G;
+
+ prof_start();
+ for (ulong ii = 0; ii < __NB_ITER; ii++)
+ {
+ nmod_geometric_progression_init(G, r, npoints_precomp, mod);
+ nmod_geometric_progression_clear(G);
+ }
+ prof_stop();
+
+ }
+
+ FLINT_TEST_CLEAR(state);
+}
+
+/* precomputations for geometric progression (eval only) */
+void sample_geometric_precomp_eval(void * arg, ulong count)
+{
+ ulong n;
+ nmod_t mod;
+ ulong i;
+
+ info_t * info = (info_t *) arg;
+ flint_bitcnt_t bits = info->bits;
+ slong npoints_precomp = info->npoints_precomp;
+
+ FLINT_TEST_INIT(state);
+
+ for (i = 0; i < count; i++)
+ {
+ n = n_randprime(state, bits, 1);
+ nmod_init(&mod, n);
+ ulong r = nmod_find_root(2*npoints_precomp, mod);
+ if (r == 0)
+ flint_printf("\n...could not find element of suitable order for geometric progression...\n");
+
+ nmod_geometric_progression_t G;
+ G->len = npoints_precomp;
+ G->mod = mod;
+
+ prof_start();
+ for (ulong ii = 0; ii < __NB_ITER; ii++)
+ {
+ _nmod_geometric_progression_init_function(G, r, npoints_precomp, mod, UWORD(1));
+ nmod_geometric_progression_clear(G);
+ }
+ prof_stop();
+ }
+
+ FLINT_TEST_CLEAR(state);
+}
+
+/* precomputations for geometric progression (interp only) */
+void sample_geometric_precomp_interp(void * arg, ulong count)
+{
+ ulong n;
+ nmod_t mod;
+ ulong i;
+
+ info_t * info = (info_t *) arg;
+ flint_bitcnt_t bits = info->bits;
+ slong npoints_precomp = info->npoints_precomp;
+
+ FLINT_TEST_INIT(state);
+
+ for (i = 0; i < count; i++)
+ {
+ n = n_randprime(state, bits, 1);
+ nmod_init(&mod, n);
+ ulong r = nmod_find_root(2*npoints_precomp, mod);
+ if (r == 0)
+ flint_printf("\n...could not find element of suitable order for geometric progression...\n");
+
+ nmod_geometric_progression_t G;
+ G->len = npoints_precomp;
+ G->mod = mod;
+
+ prof_start();
+ for (ulong ii = 0; ii < __NB_ITER; ii++)
+ {
+ _nmod_geometric_progression_init_function(G, r, npoints_precomp, mod, UWORD(2));
+ nmod_geometric_progression_clear(G);
+ }
+ prof_stop();
+ }
+
+ FLINT_TEST_CLEAR(state);
+}
+
+int main(int argc, char * argv[])
+{
+ if (argc > 1 && (strcmp(argv[1], "-h") == 0 || strcmp(argv[1], "--help") == 0))
+ {
+ flint_printf("Usage: %s -h for this help, or\n"
+ " %s [func] [short]\n"
+ " Optional arguments (if one is provided, previous ones must as well)\n"
+ " [func] is optional (default 0)\n"
+ " 0 -> all\n"
+ " 1 -> general points only\n"
+ " 2 -> geometric points only\n"
+ " [short] is optional (default 0)\n"
+ " 0 -> short bench, 1 -> full bench\n",
+ argv[0], argv[0]);
+ return 0;
+ }
+
+ const slong func_bench = (argc >= 2) ? atoi(argv[1]) : 0;
+
+ /* number of lens differs: long test / short test */
+ const slong nb_lens = (argc >= 3 && atoi(argv[2]) != 0) ? 25 : 21;
+ slong lengths[] = {1, 2, 3, 4, 6,
+ 8, 10, 12, 16, 20,
+ 30, 45, 70, 100, 200,
+ 400, 800, 1600, 3200, 6400,
+ 12800, 25600, 51200, 102400, 204800};
+
+ double tmp;
+ double time_general;
+ double time_geometric_all;
+ double time_geometric_eval;
+ double time_geometric_interp;
+
+ info_t info;
+ flint_bitcnt_t i;
+
+ flint_printf("unit: measurements in ms\n");
+
+ for (i = 63; i <= FLINT_BITS; i++)
+ {
+ info.bits = i;
+
+ printf("==== nbits = %ld====\n", i);
+
+ if (func_bench == 0) /* bench all */
+ flint_printf("len\tpoints | general | geom | g-eval | g-interp \n");
+ else if (func_bench == 1) /* general only */
+ flint_printf("len\tpoints | general\n");
+ else if (func_bench == 2) /* geometric only */
+ flint_printf("len\tpoints | geom | g-eval | g-interp \n");
+
+ for (int len = 0; len < nb_lens; ++len)
+ {
+ /* time in microsec */
+ const double fac = 1. * __NB_ITER;
+ info.npoints_precomp = lengths[len];
+ info.length = lengths[len];
+
+ if (func_bench == 0 || func_bench == 1)
+ prof_repeat(&time_general, &tmp, sample_general_precomp, (void *) &info);
+
+ if (func_bench == 0 || func_bench == 2)
+ {
+ prof_repeat(&time_geometric_all, &tmp, sample_geometric_precomp, (void *) &info);
+ prof_repeat(&time_geometric_eval, &tmp, sample_geometric_precomp_eval, (void *) &info);
+ prof_repeat(&time_geometric_interp, &tmp, sample_geometric_precomp_interp, (void *) &info);
+ }
+
+ if (func_bench == 0)
+ {
+ flint_printf("%ld\t%7ld| %.1e | %.1e | %.1e | %.1e\n",
+ info.length, info.npoints_precomp,
+ time_general/fac,
+ time_geometric_all/fac,
+ time_geometric_eval/fac,
+ time_geometric_interp/fac);
+ }
+ else if (func_bench == 1)
+ {
+ flint_printf("%ld\t%7ld| %.1e\n",
+ info.length, info.npoints_precomp,
+ time_general/fac);
+ }
+ else if (func_bench == 2)
+ {
+ flint_printf("%ld\t%7ld| %.1e | %.1e | %.1e\n",
+ info.length, info.npoints_precomp,
+ time_geometric_all/fac,
+ time_geometric_eval/fac,
+ time_geometric_interp/fac);
+ }
+ }
+
+ flint_printf("\n");
+ }
+
+ return 0;
+}
+
+#undef __NB_ITER
diff --git a/src/nmod_poly/profile/p-evaluate_nmod_vec.c b/src/nmod_poly/profile/p-evaluate_nmod_vec.c
index 9440bf0807..9fbeb63804 100644
--- a/src/nmod_poly/profile/p-evaluate_nmod_vec.c
+++ b/src/nmod_poly/profile/p-evaluate_nmod_vec.c
@@ -9,13 +9,14 @@
(at your option) any later version. See .
*/
-#include "stdlib.h" /* for atoi */
+#include /* for atoi */
+#include
+
#include "nmod.h"
#include "profiler.h"
#include "ulong_extras.h"
#include "nmod_poly.h"
#include "nmod_vec.h"
-#include
#define __NB_ITER 10
@@ -23,7 +24,7 @@
/* finds an element of order at least n */
/* returns 0 if not found */
/*------------------------------------------------------------*/
-static long nmod_find_root(slong n, nmod_t mod)
+static ulong nmod_find_root(slong n, nmod_t mod)
{
ulong attempts = 0;
for (ulong q = 2; q < mod.n; q++)
@@ -367,6 +368,46 @@ void sample_nmod_vec_geom_fast_onlyprecomp(void * arg, ulong count)
FLINT_TEST_CLEAR(state);
}
+/* multiplying two polynomials */
+void sample_nmod_poly_mul(void * arg, ulong count)
+{
+ ulong n;
+ nmod_t mod;
+ ulong i;
+
+ info_t * info = (info_t *) arg;
+ flint_bitcnt_t bits = info->bits;
+ slong length = info->length;
+
+ FLINT_TEST_INIT(state);
+
+ for (i = 0; i < count; i++)
+ {
+ n = n_randprime(state, bits, 1);
+ nmod_init(&mod, n);
+ nmod_poly_t poly1;
+ nmod_poly_t poly2;
+ nmod_poly_t poly3;
+ nmod_poly_init(poly1, n);
+ nmod_poly_init(poly2, n);
+ nmod_poly_init(poly3, n);
+ nmod_poly_randtest_monic(poly1, state, 2*length);
+ nmod_poly_randtest_monic(poly2, state, length);
+
+ prof_start();
+ for (ulong ii = 0; ii < __NB_ITER; ii++)
+ /* nmod_poly_mul(poly3, poly1, poly2); */
+ nmod_poly_mulmid(poly3, poly1, poly2, length, 2*length);
+ prof_stop();
+
+ nmod_poly_clear(poly1);
+ nmod_poly_clear(poly2);
+ nmod_poly_clear(poly3);
+ }
+
+ FLINT_TEST_CLEAR(state);
+}
+
int main(int argc, char * argv[])
{
if (argc > 1 && (strcmp(argv[1], "-h") == 0 || strcmp(argv[1], "--help") == 0))
@@ -390,9 +431,11 @@ int main(int argc, char * argv[])
const slong pre_fact = (argc >= 3) ? atoi(argv[2]) : 1;
- flint_printf("[pre_fact provided] -> will build geometric progression with `pre_fact * points` points\n");
- flint_printf("[pre_fact provided] -> note: time for general points w/ tree is for tree with `points` points\n");
- flint_printf("[pre_fact provided] -> note: time for general points precomputation is for tree with `pre_fact * points` points\n");
+ if (pre_fact != 1)
+ {
+ flint_printf("[pre_fact provided] -> will build geometric progression with `%ld * points` points\n", pre_fact);
+ flint_printf("[pre_fact provided] -> for general points (except precomp only), tree has `points` points\n");
+ }
/* number of lens differs: long test / short test */
const slong nb_lens = (argc >= 4 && atoi(argv[3]) != 0) ? 25 : 21;
@@ -404,30 +447,30 @@ int main(int argc, char * argv[])
/* specific lengths for finding out fft_small threshold */
/* slong lengths[] = { */
/* 30, */
- /* 32, */
- /* 34, */
- /* 36, */
/* 38, */
- /* 40, */
- /* 42, */
- /* 44, */
/* 46, */
- /* 48, */
- /* 50, */
- /* 52, */
/* 54, */
- /* 56, */
- /* 58, */
- /* 60, */
/* 62, */
- /* 64, */
- /* 66, */
- /* 68, */
/* 70, */
- /* 72, */
- /* 74, */
- /* 76, */
/* 78, */
+ /* 86, */
+ /* 94, */
+ /* 102, */
+ /* 110, */
+ /* 118, */
+ /* 126, */
+ /* 134, */
+ /* 142, */
+ /* 150, */
+ /* 158, */
+ /* 166, */
+ /* 174, */
+ /* 182, */
+ /* 190, */
+ /* 198, */
+ /* 206, */
+ /* 214, */
+ /* 222, */
/* }; */
double tmp;
@@ -439,6 +482,7 @@ int main(int argc, char * argv[])
double time_geom_fast;
double time_geom_fast_precomp;
double time_geom_fast_onlyprecomp;
+ double time_poly_mul;
info_t info;
flint_bitcnt_t i;
@@ -455,20 +499,23 @@ int main(int argc, char * argv[])
{
if (func_bench == 0) /* bench all */
{
- flint_printf("==== nb eval points == %d * (poly length) ====\n", npoints_factor);
- flint_printf("len\tpoints | GENERAL POINTS | GEOMETRIC PROGRESSION \n");
+ flint_printf("==== nb eval points : %d * (poly length) ====\n", npoints_factor);
+ flint_printf("==== nb precomp points : %d * (poly length) ====\n", pre_fact * npoints_factor);
+ flint_printf("len\tpoints | GENERAL POINTS | GEOMETRIC PROGRESSION | POLY_MUL\n");
flint_printf("len\tpoints |iter\tfast\tw/ tree\ttree |iter\tfast\tw/ prec\tprecomp\n");
}
else if (func_bench == 1) /* general only */
{
- flint_printf("==== nb eval points == %d * (poly length) ====\n", npoints_factor);
+ flint_printf("==== nb eval points : %d * (poly length) ====\n", npoints_factor);
+ flint_printf("==== nb precomp points : %d * (poly length) ====\n", pre_fact * npoints_factor);
flint_printf("len\tpoints | GENERAL POINTS |\n");
flint_printf("len\tpoints |iter\tfast\tw/ tree\ttree |\n");
}
else if (func_bench == 2) /* geometric only */
{
- flint_printf("==== nb eval points == %d * (poly length) ====\n", npoints_factor);
- flint_printf("len\tpoints | GEOMETRIC PROGRESSION \n");
+ flint_printf("==== nb eval points : %d * (poly length) ====\n", npoints_factor);
+ flint_printf("==== nb precomp points : %d * (poly length) ====\n", pre_fact * npoints_factor);
+ flint_printf("len\tpoints | GEOMETRIC PROGRESSION | POLY_MUL\n");
flint_printf("len\tpoints |iter\tfast\tw/ prec\tprecomp\n");
}
@@ -510,14 +557,17 @@ int main(int argc, char * argv[])
prof_repeat(&time_geom_fast_precomp, &tmp, sample_nmod_vec_geom_fast_precomp, (void *) &info);
prof_repeat(&time_geom_fast_onlyprecomp, &tmp, sample_nmod_vec_geom_fast_onlyprecomp, (void *) &info);
+
+ prof_repeat(&time_poly_mul, &tmp, sample_nmod_poly_mul, (void *) &info);
}
if (func_bench == 0)
{
- flint_printf("%ld\t%7ld|%.1e\t%.1e\t%.1e\t%.1e|%.1e\t%.1e\t%.1e\t%.1e\n",
+ flint_printf("%ld\t%7ld|%.1e\t%.1e\t%.1e\t%.1e|%.1e\t%.1e\t%.1e\t%.1e\t|%.1e\n",
info.length, info.npoints,
time_iter/fac, time_fast/fac, time_fast_precomp/fac, time_fast_onlyprecomp/fac,
- time_geom_iter/fac, time_geom_fast/fac, time_geom_fast_precomp/fac, time_geom_fast_onlyprecomp/fac);
+ time_geom_iter/fac, time_geom_fast/fac, time_geom_fast_precomp/fac, time_geom_fast_onlyprecomp/fac,
+ time_poly_mul/fac);
}
else if (func_bench == 1)
{
@@ -527,9 +577,10 @@ int main(int argc, char * argv[])
}
else if (func_bench == 2)
{
- flint_printf("%ld\t%7ld|%.1e\t%.1e\t%.1e\t%.1e\n",
+ flint_printf("%ld\t%7ld|%.1e\t%.1e\t%.1e\t%.1e\t|%.1e\n",
info.length, info.npoints,
- time_geom_iter/fac, time_geom_fast/fac, time_geom_fast_precomp/fac, time_geom_fast_onlyprecomp/fac);
+ time_geom_iter/fac, time_geom_fast/fac, time_geom_fast_precomp/fac, time_geom_fast_onlyprecomp/fac,
+ time_poly_mul/fac);
}
}
}
diff --git a/src/nmod_poly/profile/p-interpolate_nmod_vec.c b/src/nmod_poly/profile/p-interpolate_nmod_vec.c
new file mode 100644
index 0000000000..7ff726b26b
--- /dev/null
+++ b/src/nmod_poly/profile/p-interpolate_nmod_vec.c
@@ -0,0 +1,592 @@
+/*
+ Copyright (C) 2026 Vincent Neiger
+
+ This file is part of FLINT.
+
+ FLINT is free software: you can redistribute it and/or modify it under
+ the terms of the GNU Lesser General Public License (LGPL) as published
+ by the Free Software Foundation; either version 3 of the License, or
+ (at your option) any later version. See .
+*/
+
+#include /* for atoi */
+#include
+
+#include "nmod.h"
+#include "profiler.h"
+#include "ulong_extras.h"
+#include "nmod_poly.h"
+#include "nmod_vec.h"
+
+#define __NB_ITER 10
+
+/*------------------------------------------------------------*/
+/* finds an element of order at least n */
+/* returns 0 if not found */
+/*------------------------------------------------------------*/
+static ulong nmod_find_root(slong n, nmod_t mod)
+{
+ ulong attempts = 0;
+ for (ulong q = 2; q < mod.n; q++)
+ {
+ slong k = 1;
+ slong qk = q;
+ while (qk != 1 && k < n)
+ {
+ qk = nmod_mul(qk, q, mod);
+ k++;
+ }
+ if (qk != 1)
+ {
+ return q;
+ }
+ attempts += 1;
+ if (attempts >= 10)
+ return 0;
+ }
+ return 0;
+}
+
+typedef struct
+{
+ flint_bitcnt_t bits;
+ slong length;
+ slong npoints_precomp;
+} info_t;
+
+/* Newton */
+void sample_nmod_vec_newton(void * arg, ulong count)
+{
+ ulong n;
+ nmod_t mod;
+ ulong i;
+
+ info_t * info = (info_t *) arg;
+ flint_bitcnt_t bits = info->bits;
+ slong length = info->length;
+
+ FLINT_TEST_INIT(state);
+
+ nn_ptr pts = _nmod_vec_init(length);
+ nn_ptr vals = _nmod_vec_init(length);
+ nn_ptr poly = _nmod_vec_init(length);
+
+ for (i = 0; i < count; i++)
+ {
+ n = n_randprime(state, bits, 0);
+ if (n == UWORD(0)) n++;
+ nmod_init(&mod, n);
+ _nmod_vec_rand(vals, state, length, mod);
+ _nmod_vec_rand(pts, state, length, mod);
+
+ prof_start();
+ for (ulong ii = 0; ii < __NB_ITER; ii++)
+ _nmod_poly_interpolate_nmod_vec_newton(poly, pts, vals, length, mod);
+ prof_stop();
+ }
+
+ _nmod_vec_clear(pts);
+ _nmod_vec_clear(vals);
+ _nmod_vec_clear(poly);
+
+ FLINT_TEST_CLEAR(state);
+}
+
+/* barycentric */
+void sample_nmod_vec_barycentric(void * arg, ulong count)
+{
+ ulong n;
+ nmod_t mod;
+ ulong i;
+
+ info_t * info = (info_t *) arg;
+ flint_bitcnt_t bits = info->bits;
+ slong length = info->length;
+
+ FLINT_TEST_INIT(state);
+
+ nn_ptr pts = _nmod_vec_init(length);
+ nn_ptr vals = _nmod_vec_init(length);
+ nn_ptr poly = _nmod_vec_init(length);
+
+ for (i = 0; i < count; i++)
+ {
+ n = n_randprime(state, bits, 0);
+ if (n == UWORD(0)) n++;
+ nmod_init(&mod, n);
+ _nmod_vec_rand(vals, state, length, mod);
+ _nmod_vec_rand(pts, state, length, mod);
+
+ prof_start();
+ for (ulong ii = 0; ii < __NB_ITER; ii++)
+ _nmod_poly_interpolate_nmod_vec_barycentric(poly, pts, vals, length, mod);
+ prof_stop();
+ }
+
+ _nmod_vec_clear(pts);
+ _nmod_vec_clear(vals);
+ _nmod_vec_clear(poly);
+
+ FLINT_TEST_CLEAR(state);
+}
+
+/* fast */
+void sample_nmod_vec_fast(void * arg, ulong count)
+{
+ ulong n;
+ nmod_t mod;
+ ulong i;
+
+ info_t * info = (info_t *) arg;
+ flint_bitcnt_t bits = info->bits;
+ slong length = info->length;
+
+ FLINT_TEST_INIT(state);
+
+ nn_ptr pts = _nmod_vec_init(length);
+ nn_ptr vals = _nmod_vec_init(length);
+ nn_ptr poly = _nmod_vec_init(length);
+
+ for (i = 0; i < count; i++)
+ {
+ n = n_randprime(state, bits, 0);
+ if (n == UWORD(0)) n++;
+ nmod_init(&mod, n);
+ _nmod_vec_rand(vals, state, length, mod);
+ _nmod_vec_rand(pts, state, length, mod);
+
+ prof_start();
+ for (ulong ii = 0; ii < __NB_ITER; ii++)
+ _nmod_poly_interpolate_nmod_vec_fast(poly, pts, vals, length, mod);
+ prof_stop();
+ }
+
+ _nmod_vec_clear(pts);
+ _nmod_vec_clear(vals);
+ _nmod_vec_clear(poly);
+
+ FLINT_TEST_CLEAR(state);
+}
+
+/* fast, tree, not counting precomputations */
+void sample_nmod_vec_fast_precomp(void * arg, ulong count)
+{
+ ulong n;
+ nmod_t mod;
+ ulong i;
+
+ info_t * info = (info_t *) arg;
+ flint_bitcnt_t bits = info->bits;
+ slong length = info->length;
+
+ FLINT_TEST_INIT(state);
+
+ nn_ptr pts = _nmod_vec_init(length);
+ nn_ptr vals = _nmod_vec_init(length);
+ nn_ptr poly = _nmod_vec_init(length);
+
+ for (i = 0; i < count; i++)
+ {
+ n = n_randprime(state, bits, 0);
+ if (n == UWORD(0)) n++;
+ nmod_init(&mod, n);
+ _nmod_vec_rand(vals, state, length, mod);
+ _nmod_vec_rand(pts, state, length, mod);
+
+ nn_ptr * tree = _nmod_poly_tree_alloc(length);
+ _nmod_poly_tree_build(tree, pts, length, mod);
+
+ nn_ptr w = _nmod_vec_init(length);
+ _nmod_poly_interpolation_weights(w, tree, length, mod);
+
+ prof_start();
+ for (ulong ii = 0; ii < __NB_ITER; ii++)
+ _nmod_poly_interpolate_nmod_vec_fast_precomp(poly, vals, tree, w, length, mod);
+ prof_stop();
+
+ _nmod_poly_tree_free(tree, length);
+ flint_free(w);
+ }
+
+ _nmod_vec_clear(pts);
+ _nmod_vec_clear(vals);
+ _nmod_vec_clear(poly);
+
+ FLINT_TEST_CLEAR(state);
+}
+
+/* tree, counting only precomputations */
+void sample_nmod_vec_fast_onlyprecomp(void * arg, ulong count)
+{
+ ulong n;
+ nmod_t mod;
+ ulong i;
+
+ info_t * info = (info_t *) arg;
+ flint_bitcnt_t bits = info->bits;
+ slong length = info->length;
+
+ FLINT_TEST_INIT(state);
+
+ nn_ptr pts = _nmod_vec_init(length);
+
+ for (i = 0; i < count; i++)
+ {
+ n = n_randprime(state, bits, 0);
+ if (n == UWORD(0)) n++;
+ nmod_init(&mod, n);
+ _nmod_vec_rand(pts, state, length, mod);
+
+ prof_start();
+ for (ulong ii = 0; ii < __NB_ITER; ii++)
+ {
+ nn_ptr * tree = _nmod_poly_tree_alloc(length);
+ _nmod_poly_tree_build(tree, pts, length, mod);
+
+ nn_ptr w = _nmod_vec_init(length);
+ _nmod_poly_interpolation_weights(w, tree, length, mod);
+
+ _nmod_poly_tree_free(tree, length);
+ flint_free(w);
+ }
+ prof_stop();
+ }
+
+ _nmod_vec_clear(pts);
+
+ FLINT_TEST_CLEAR(state);
+}
+
+/* geometric: fast, with precomp (by default, on `length` points) */
+void sample_nmod_vec_geom_fast(void * arg, ulong count)
+{
+ ulong n;
+ nmod_t mod;
+ ulong i;
+
+ info_t * info = (info_t *) arg;
+ flint_bitcnt_t bits = info->bits;
+ slong length = info->length;
+
+ FLINT_TEST_INIT(state);
+
+ nn_ptr vals = _nmod_vec_init(length);
+
+ for (i = 0; i < count; i++)
+ {
+ n = n_randprime(state, bits, 0);
+ nmod_init(&mod, n);
+ nmod_poly_t poly;
+ nmod_poly_init(poly, n);
+
+ _nmod_vec_rand(vals, state, length, mod);
+ ulong r = nmod_find_root(2*length, mod);
+ if (r == 0)
+ flint_printf("\n...could not find element of suitable order for geometric progression...\n");
+
+ prof_start();
+ for (ulong ii = 0; ii < __NB_ITER; ii++)
+ nmod_poly_interpolate_geometric_nmod_vec_fast(poly, r, vals, length);
+ prof_stop();
+ nmod_poly_clear(poly);
+ }
+
+ _nmod_vec_clear(vals);
+
+ FLINT_TEST_CLEAR(state);
+}
+
+/* geometric: fast, not counting precomputations (which is on npoints_precomp points) */
+void sample_nmod_vec_geom_fast_precomp(void * arg, ulong count)
+{
+ ulong n;
+ nmod_t mod;
+ ulong i;
+
+ info_t * info = (info_t *) arg;
+ flint_bitcnt_t bits = info->bits;
+ slong length = info->length;
+ slong npoints_precomp = info->npoints_precomp;
+
+ FLINT_TEST_INIT(state);
+
+ nn_ptr vals = _nmod_vec_init(length);
+
+ for (i = 0; i < count; i++)
+ {
+ n = n_randprime(state, bits, 1);
+ nmod_init(&mod, n);
+ nn_ptr poly = _nmod_vec_init(length);
+ _nmod_vec_rand(vals, state, length, mod);
+ ulong r = nmod_find_root(2*npoints_precomp, mod);
+ if (r == 0)
+ flint_printf("\n...could not find element of suitable order for geometric progression...\n");
+
+ nmod_geometric_progression_t G;
+ nmod_geometric_progression_init(G, r, npoints_precomp, mod);
+
+ prof_start();
+ for (ulong ii = 0; ii < __NB_ITER; ii++)
+ _nmod_poly_interpolate_geometric_nmod_vec_fast_precomp(poly, vals, G, length, mod);
+ prof_stop();
+
+ nmod_geometric_progression_clear(G);
+ _nmod_vec_clear(poly);
+ }
+
+ _nmod_vec_clear(vals);
+
+ FLINT_TEST_CLEAR(state);
+}
+
+/* geometric: fast, counting only precomputations (on npoints_precomp points) */
+void sample_nmod_vec_geom_fast_onlyprecomp(void * arg, ulong count)
+{
+ ulong n;
+ nmod_t mod;
+ ulong i;
+
+ info_t * info = (info_t *) arg;
+ flint_bitcnt_t bits = info->bits;
+ slong npoints_precomp = info->npoints_precomp;
+
+ FLINT_TEST_INIT(state);
+
+ for (i = 0; i < count; i++)
+ {
+ n = n_randprime(state, bits, 1);
+ nmod_init(&mod, n);
+ ulong r = nmod_find_root(2*npoints_precomp, mod);
+ if (r == 0)
+ flint_printf("\n...could not find element of suitable order for geometric progression...\n");
+
+ nmod_geometric_progression_t G;
+
+ prof_start();
+ for (ulong ii = 0; ii < __NB_ITER; ii++)
+ _nmod_geometric_progression_init_function(G, r, npoints_precomp, mod, UWORD(2));
+ prof_stop();
+
+ nmod_geometric_progression_clear(G);
+ }
+
+ FLINT_TEST_CLEAR(state);
+}
+
+/* multiplying two polynomials */
+void sample_nmod_poly_mul(void * arg, ulong count)
+{
+ ulong n;
+ nmod_t mod;
+ ulong i;
+
+ info_t * info = (info_t *) arg;
+ flint_bitcnt_t bits = info->bits;
+ slong length = info->length;
+
+ FLINT_TEST_INIT(state);
+
+ for (i = 0; i < count; i++)
+ {
+ n = n_randprime(state, bits, 1);
+ nmod_init(&mod, n);
+ nmod_poly_t poly1;
+ nmod_poly_t poly2;
+ nmod_poly_t poly3;
+ nmod_poly_init(poly1, n);
+ nmod_poly_init(poly2, n);
+ nmod_poly_init(poly3, n);
+ nmod_poly_randtest_monic(poly1, state, length);
+ nmod_poly_randtest_monic(poly2, state, length);
+
+ prof_start();
+ for (ulong ii = 0; ii < __NB_ITER; ii++)
+ nmod_poly_mul(poly3, poly1, poly2);
+ /* nmod_poly_mulmid(poly3, poly1, poly2, length, 2*length); */
+ prof_stop();
+
+ nmod_poly_clear(poly1);
+ nmod_poly_clear(poly2);
+ nmod_poly_clear(poly3);
+ }
+
+ FLINT_TEST_CLEAR(state);
+}
+
+int main(int argc, char * argv[])
+{
+ if (argc > 1 && (strcmp(argv[1], "-h") == 0 || strcmp(argv[1], "--help") == 0))
+ {
+ flint_printf("Usage: %s -h for this help, or\n"
+ " %s [func] [pre_fact] [short]\n"
+ " Optional arguments (if one is provided, previous ones must as well)\n"
+ " [func] is optional (default 0)\n"
+ " 0 -> all\n"
+ " 1 -> general points only\n"
+ " 2 -> geometric points only\n"
+ " [pre_fact] is optional (default 1)\n"
+ " positive integer; precomputation will be on pre_fact * npoints\n"
+ " [short] is optional (default 0)\n"
+ " 0 -> short bench, 1 -> full bench\n",
+ argv[0], argv[0]);
+ return 0;
+ }
+
+ const slong func_bench = (argc >= 2) ? atoi(argv[1]) : 0;
+
+ const slong pre_fact = (argc >= 3) ? atoi(argv[2]) : 1;
+
+ if (pre_fact != 1)
+ {
+ flint_printf("[pre_fact provided] -> will build geometric progression with `%ld * points` points\n", pre_fact);
+ flint_printf("[pre_fact provided] -> for general points (except precomp only), tree has `points` points\n");
+ }
+
+ /* number of lens differs: long test / short test */
+ const slong nb_lens = (argc >= 4 && atoi(argv[3]) != 0) ? 25 : 21;
+ slong lengths[] = {1, 2, 3, 4, 6,
+ 8, 10, 12, 16, 20,
+ 30, 45, 70, 100, 200,
+ 400, 800, 1600, 3200, 6400,
+ 12800, 25600, 51200, 102400, 204800};
+ /* specific lengths for finding out fft_small threshold */
+ /* slong lengths[] = { */
+ /* 30, */
+ /* 32, */
+ /* 34, */
+ /* 36, */
+ /* 38, */
+ /* 40, */
+ /* 42, */
+ /* 44, */
+ /* 46, */
+ /* 48, */
+ /* 50, */
+ /* 52, */
+ /* 54, */
+ /* 56, */
+ /* 58, */
+ /* 60, */
+ /* 62, */
+ /* 64, */
+ /* 66, */
+ /* 68, */
+ /* 70, */
+ /* 72, */
+ /* 74, */
+ /* 76, */
+ /* 78, */
+ /* }; */
+
+ double tmp;
+ double time_newton;
+ double time_barycentric;
+ double time_fast;
+ double time_fast_precomp;
+ double time_fast_onlyprecomp;
+ double time_geom_fast;
+ double time_geom_fast_precomp;
+ double time_geom_fast_onlyprecomp;
+ double time_poly_mul;
+
+ info_t info;
+ flint_bitcnt_t i;
+
+ flint_printf("unit: measurements in ms\n");
+
+ for (i = 63; i <= FLINT_BITS; i++)
+ {
+ info.bits = i;
+
+ printf("==== nbits = %ld====\n", i);
+
+ if (func_bench == 0) /* bench all */
+ {
+ flint_printf("==== nb precomp points : %d * (poly length) ====\n", pre_fact);
+ flint_printf("len\tpoints | GENERAL POINTS | GEOMETRIC PROGRESSION | POLY_MUL\n");
+ flint_printf("len\tpoints |newton\tbaryc\tfast\tw/ tree\ttree\t| fast\tw/ prec\tprecomp\n");
+ }
+ else if (func_bench == 1) /* general only */
+ {
+ flint_printf("==== nb precomp points : %d * (poly length) ====\n", pre_fact);
+ flint_printf("len\tpoints | GENERAL POINTS |\n");
+ flint_printf("len\tpoints |newton\tbaryc\tfast\tw/ tree\ttree\t|\n");
+ }
+ else if (func_bench == 2) /* geometric only */
+ {
+ flint_printf("==== nb precomp points : %d * (poly length) ====\n", pre_fact);
+ flint_printf("len\tpoints | GEOMETRIC PROGRESSION | POLY_MUL\n");
+ flint_printf("len\tpoints |fast\tw/ prec\tprecomp\n");
+ }
+
+ for (int len = 0; len < nb_lens; ++len)
+ {
+ /* cycles / limb (if CLOCK_SCALE_FACTOR set correctly) */
+ /* const double fac = npoints[len] * (double)FLINT_CLOCK_SCALE_FACTOR; */
+ /* time in s */
+ /* const double fac = 1000000. * __NB_ITER; */
+ /* time in ms */
+ const double fac = 1. * __NB_ITER;
+ info.npoints_precomp = pre_fact * lengths[len];
+ info.length = lengths[len];
+
+ if (func_bench == 0 || func_bench == 1)
+ {
+ if (info.length <= 512)
+ prof_repeat(&time_newton, &tmp, sample_nmod_vec_newton, (void *) &info);
+ else
+ time_newton = 0.0;
+
+ if (info.length <= 1500)
+ prof_repeat(&time_barycentric, &tmp, sample_nmod_vec_barycentric, (void *) &info);
+ else
+ time_barycentric = 0.0;
+
+ prof_repeat(&time_fast, &tmp, sample_nmod_vec_fast, (void *) &info);
+
+ prof_repeat(&time_fast_precomp, &tmp, sample_nmod_vec_fast_precomp, (void *) &info);
+
+ prof_repeat(&time_fast_onlyprecomp, &tmp, sample_nmod_vec_fast_onlyprecomp, (void *) &info);
+ }
+
+ if (func_bench == 0 || func_bench == 2)
+ {
+ prof_repeat(&time_geom_fast, &tmp, sample_nmod_vec_geom_fast, (void *) &info);
+
+ prof_repeat(&time_geom_fast_precomp, &tmp, sample_nmod_vec_geom_fast_precomp, (void *) &info);
+
+ prof_repeat(&time_geom_fast_onlyprecomp, &tmp, sample_nmod_vec_geom_fast_onlyprecomp, (void *) &info);
+
+ prof_repeat(&time_poly_mul, &tmp, sample_nmod_poly_mul, (void *) &info);
+ }
+
+ if (func_bench == 0)
+ {
+ flint_printf("%ld\t%7ld|%.1e %.1e %.1e %.1e %.1e |%.1e %.1e %.1e |%.1e\n",
+ info.length, info.npoints_precomp,
+ time_newton/fac, time_barycentric/fac, time_fast/fac, time_fast_precomp/fac, time_fast_onlyprecomp/fac,
+ time_geom_fast/fac, time_geom_fast_precomp/fac, time_geom_fast_onlyprecomp/fac,
+ time_poly_mul/fac);
+ }
+ else if (func_bench == 1)
+ {
+ flint_printf("%ld\t%7ld|%.1e %.1e %.1e %.1e %.1e\n",
+ info.length, info.npoints_precomp,
+ time_newton/fac, time_barycentric/fac, time_fast/fac, time_fast_precomp/fac, time_fast_onlyprecomp/fac);
+ }
+ else if (func_bench == 2)
+ {
+ flint_printf("%ld\t%7ld|%.1e %.1e %.1e |%.1e\n",
+ info.length, info.npoints_precomp,
+ time_geom_fast/fac, time_geom_fast_precomp/fac, time_geom_fast_onlyprecomp/fac,
+ time_poly_mul/fac);
+ }
+ }
+
+ flint_printf("\n");
+ }
+
+ return 0;
+}
+
+#undef __NB_ITER
diff --git a/src/nmod_poly/test/t-interpolate_geometric_nmod_vec_fast.c b/src/nmod_poly/test/t-interpolate_geometric_nmod_vec_fast.c
index 1192ce2549..d461e6d612 100644
--- a/src/nmod_poly/test/t-interpolate_geometric_nmod_vec_fast.c
+++ b/src/nmod_poly/test/t-interpolate_geometric_nmod_vec_fast.c
@@ -23,41 +23,98 @@ TEST_FUNCTION_START(nmod_poly_interpolate_geometric_nmod_vec_fast, state)
nmod_poly_t P, Q;
nn_ptr y;
ulong mod, r;
- slong n, npoints;
+ slong len, npoints;
- npoints = 1 + (i == 0 ? 1 : n_randint(state, 100));
- n = 1 + n_randint(state, npoints);
+ /* npoints = 1 + (i == 0 ? 1 : n_randint(state, 200)); */
+ npoints = 1 + n_randint(state, 200);
+ len = 1 + n_randint(state, npoints);
do
{
mod = n_randtest_prime(state, 1);
}
- while (mod <= 2*FLINT_MAX(npoints, n) + 1); // minimum limit for maximum order r
-
- nmod_poly_init(P, mod);
- nmod_poly_init(Q, mod);
+ while (mod <= (ulong)(2*FLINT_MAX(npoints, len) + 1)); // minimum limit for maximum order r
r = n_primitive_root_prime(mod);
y = _nmod_vec_init(npoints);
- nmod_poly_randtest(P, state, n);
+ /* use full `npoints` points */
+ {
+ nmod_poly_init(P, mod);
+ nmod_poly_init(Q, mod);
+ nmod_poly_randtest(P, state, len);
+
+ nmod_poly_evaluate_geometric_nmod_vec_fast(y, P, r, npoints);
+ nmod_poly_interpolate_geometric_nmod_vec_fast(Q, r, y, npoints);
+ result = nmod_poly_equal(P, Q);
+ if (!result)
+ {
+ flint_printf("FAIL (all points):\n");
+ flint_printf("mod=%wu, len=%wd, npoints=%wd\n\n", mod, len, npoints);
+ nmod_poly_print_pretty(P, "x"), flint_printf("\n\n");
+ nmod_poly_print_pretty(Q, "x"), flint_printf("\n\n");
+ fflush(stdout);
+ flint_abort();
+ }
- nmod_poly_evaluate_geometric_nmod_vec_fast(y, P, r, npoints);
- nmod_poly_interpolate_geometric_nmod_vec_fast(Q, r, y, npoints);
+ nmod_poly_clear(P);
+ nmod_poly_clear(Q);
+ }
- result = nmod_poly_equal(P, Q);
- if (!result)
+ /* use only `len` points */
{
- flint_printf("FAIL:\n");
- flint_printf("mod=%wu, n=%wd, npoints=%wd\n\n", mod, n, npoints);
- nmod_poly_print_pretty(P, "x"), flint_printf("\n\n");
- nmod_poly_print_pretty(Q, "x"), flint_printf("\n\n");
- fflush(stdout);
- flint_abort();
+ nmod_poly_init(P, mod);
+ nmod_poly_init(Q, mod);
+ nmod_poly_randtest(P, state, len);
+
+ nmod_poly_evaluate_geometric_nmod_vec_fast(y, P, r, npoints);
+ nmod_poly_interpolate_geometric_nmod_vec_fast(Q, r, y, len);
+ result = nmod_poly_equal(P, Q);
+ if (!result)
+ {
+ flint_printf("FAIL (`len` points):\n");
+ flint_printf("mod=%wu, len=%wd, npoints=%wd\n\n", mod, len, npoints);
+ nmod_poly_print_pretty(P, "x"), flint_printf("\n\n");
+ nmod_poly_print_pretty(Q, "x"), flint_printf("\n\n");
+ fflush(stdout);
+ flint_abort();
+ }
+
+ nmod_poly_clear(P);
+ nmod_poly_clear(Q);
+ }
+
+ /* use variant with given precomputation */
+ {
+ nmod_t mod2;
+ nmod_init(&mod2, mod);
+ nmod_geometric_progression_t G;
+ nmod_geometric_progression_init(G, r, npoints, mod2);
+ nmod_poly_init(P, mod);
+ nmod_poly_init(Q, mod);
+ nmod_poly_randtest(P, state, len);
+
+ _nmod_poly_evaluate_geometric_nmod_vec_fast_precomp(y, P->coeffs, P->length, G, npoints, mod2);
+ nmod_poly_fit_length(Q, len);
+ _nmod_poly_interpolate_geometric_nmod_vec_fast_precomp(Q->coeffs, y, G, len, mod2);
+ _nmod_poly_set_length(Q, len);
+ _nmod_poly_normalise(Q);
+ result = nmod_poly_equal(P, Q);
+ if (!result)
+ {
+ flint_printf("FAIL (precomp):\n");
+ flint_printf("mod=%wu, len=%wd, npoints=%wd\n\n", mod, len, npoints);
+ nmod_poly_print_pretty(P, "x"), flint_printf("\n\n");
+ nmod_poly_print_pretty(Q, "x"), flint_printf("\n\n");
+ fflush(stdout);
+ flint_abort();
+ }
+
+ nmod_poly_clear(P);
+ nmod_poly_clear(Q);
+ nmod_geometric_progression_clear(G);
}
- nmod_poly_clear(P);
- nmod_poly_clear(Q);
- _nmod_vec_clear(y);
+ flint_free(y);
}
TEST_FUNCTION_END(state);