From 57d6963c122909ed61e7c82fb33b2f22627edfcc Mon Sep 17 00:00:00 2001 From: Vincent Neiger Date: Fri, 1 May 2026 16:12:06 +0200 Subject: [PATCH 01/13] add profile for precomputation of evaluate/interpolate + clean/improve geometric init --- src/nmod_poly.h | 8 +- src/nmod_poly/evaluate_geometric_nmod_vec.c | 38 +-- src/nmod_poly/geometric_progression.c | 203 ++++++++------- .../interpolate_geometric_nmod_vec.c | 14 +- .../profile/p-eval_interp_precomputations.c | 245 ++++++++++++++++++ 5 files changed, 382 insertions(+), 126 deletions(-) create mode 100644 src/nmod_poly/profile/p-eval_interp_precomputations.c diff --git a/src/nmod_poly.h b/src/nmod_poly.h index cfe9154b5e..978d914a2e 100644 --- a/src/nmod_poly.h +++ b/src/nmod_poly.h @@ -60,13 +60,15 @@ 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 + nn_ptr ev_s; /* evaluate: scaling constants */ + nmod_poly_t ev_f; /* evaluate: polynomial */ + nn_ptr int_s1, int_s2, int_s3; /* interpolate: scaling constants */ + nmod_poly_t int_f1, int_f2; /* interpolate: polynomials */ nmod_t mod; slong len; // number of points - } nmod_geometric_progression_struct; + typedef nmod_geometric_progression_struct nmod_geometric_progression_t[1]; /* Memory management ********************************************************/ diff --git a/src/nmod_poly/evaluate_geometric_nmod_vec.c b/src/nmod_poly/evaluate_geometric_nmod_vec.c index c9117a2a26..a1c42f0f41 100644 --- a/src/nmod_poly/evaluate_geometric_nmod_vec.c +++ b/src/nmod_poly/evaluate_geometric_nmod_vec.c @@ -39,17 +39,25 @@ void _nmod_poly_evaluate_geometric_nmod_vec_fast_precomp(nn_ptr vs, nn_srcptr po } } + /** 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 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)) - **/ 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) + /** this uses a short product: write rev(a) = x**(alen-1) a(1/x), and write F = (G->ev_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 @@ -65,14 +73,14 @@ void _nmod_poly_evaluate_geometric_nmod_vec_fast_precomp(nn_ptr vs, nn_srcptr po nn_ptr b = _nmod_vec_init(blen); for (slong i = val; i < plen; i++) - a[i - val] = nmod_mul(G->x[i], poly[i], mod); + a[i - val] = 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_reverse(Frev, G->ev_f->coeffs + val, blen, blen); _nmod_poly_mullow(b, Frev, blen, a, alen, blen, 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[alen - 1 + len - 1 - i], mod); _nmod_vec_clear(Frev); _nmod_vec_clear(a); @@ -81,18 +89,18 @@ void _nmod_poly_evaluate_geometric_nmod_vec_fast_precomp(nn_ptr vs, nn_srcptr po else { /* version 2 */ - /* this uses a middle product to compute [rev(p) * G->f]_{plen - 1}^{len} (i.e. coeffs [plen - 1, plen - 1 + len)) */ + /* 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 /* 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); + b[plen - 1 - i] = nmod_mul(G->ev_s[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); + _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[i], mod); + vs[i] = nmod_mul(G->ev_s[i], b[i], mod); _nmod_vec_clear(b); #else @@ -102,12 +110,12 @@ void _nmod_poly_evaluate_geometric_nmod_vec_fast_precomp(nn_ptr vs, nn_srcptr po nn_ptr b = _nmod_vec_init(alen + (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_mulhigh(b, G->ev_f->coeffs + val, alen - 1 + len, a, alen, alen - 1, 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[alen - 1 + i], mod); _nmod_vec_clear(a); _nmod_vec_clear(b); diff --git a/src/nmod_poly/geometric_progression.c b/src/nmod_poly/geometric_progression.c index 3537a2c906..fc874f1443 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. @@ -13,124 +13,123 @@ #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) +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) { - 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; - - 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 = 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++) { - G->f->coeffs[i] = nmod_mul(G->f->coeffs[i - 1], tmp, mod); - tmp = nmod_mul(tmp, q, mod); + G->ev_f->coeffs[i] = nmod_mul(G->ev_f->coeffs[i - 1], r, mod); + r = nmod_mul(r, q, mod); } - // 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++) - { - G->x[i] = nmod_mul(G->x[i - 1], tmp, mod); - tmp = nmod_mul(tmp, inv_q, mod); - } - - 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++) + G->ev_s[0] = 1; + for (slong i = 1; i < len; 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_s[i] = nmod_mul(G->ev_s[i - 1], inv_r, mod); + inv_r = nmod_mul(inv_r, inv_q, mod); } - - tmp = nmod_inv(prod_diff[len-1], mod); - for (i = len - 1; i > 0; i--) +} + +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) */ + + /* coeff(G->int_f1, i) = (-1)**i * q_i / u_i */ + /* coeff(G->int_f2, i) = q_i / u_i, i.e. G->int_f2 == G->int_f1(-x) */ + nmod_poly_init2_preinv(G->int_f1, mod.n, mod.ninv, len); + nmod_poly_init2_preinv(G->int_f2, mod.n, mod.ninv, len); + G->int_f1->length = len; + G->int_f2->length = len; + + /* G->int_s2[i] = (-1)**i * u_i / q_i = inverse of coeff(G->int_f1, i) */ + /* G->int_s3[i] = (-1)**i / u_i */ + /* G->int_s1[i] = 1 / u_i */ + G->int_s1 = _nmod_vec_init(len); + G->int_s2 = _nmod_vec_init(len); + G->int_s3 = _nmod_vec_init(len); + + G->int_f1->coeffs[0] = 1; + G->int_f2->coeffs[0] = 1; + G->int_s1[0] = 1; + G->int_s2[0] = 1; + G->int_s3[0] = 1; + + ulong q_pow_i = 1; + ulong inv_q_pow_i = 1; + ulong inv_q_i = 1; + ulong prod_diff = 1; + + /* TODO improve with Shoup */ + 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); + 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_f2->coeffs[i] = nmod_mul(G->int_f2->coeffs[i-1], q_pow_i, mod); /* q_i */ + q_pow_i = nmod_mul(q_pow_i, q, mod); /* q**i */ + G->int_s3[i] = q_pow_i - 1; /* temporarily, G->int_s3[i] = q**i - 1 */ + prod_diff = nmod_mul(q_pow_i - 1, prod_diff, mod); /* u_i */ + if (i % 2) /* i is odd */ + G->int_s2[i] = nmod_mul(mod.n - prod_diff, inv_q_i, mod); /* (-1)**i * u_i / q_i */ + else /* i is even */ + G->int_s2[i] = nmod_mul(prod_diff, inv_q_i, mod); /* (-1)**i * u_i / q_i */ } - inv_diff[0] = tmp; - // end montgomery inversion - - // sets sequences w, y, z and polynomials g1, g2 - qk = 1; - inv_qk = 1; - qq = 1; - s = 1; - for (i = 1; i < len; 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 */ + ulong w_i = G->int_s1[i]; /* 1 / u_i */ + ulong tmp = nmod_mul(G->int_f2->coeffs[i], w_i, mod); + G->int_f2->coeffs[i] = tmp; /* q_i / u_i */ + G->int_s1[i-1] = nmod_mul(G->int_s3[i], w_i, mod); /* 1 / u_{i-1} */ + if (i % 2) /* 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]; + G->int_s3[i] = mod.n - w_i; /* (-1)**i * G->int_s1[i] */ + G->int_f1->coeffs[i] = mod.n - tmp; /* (-1)**i * G->int_f2[i] */ } - else /* i is even */ + else /* i is even */ { - G->g1->coeffs[i] = tmp; - G->y[i] = prod_diff[i]; - G->z[i] = G->w[i]; + G->int_s3[i] = w_i; /* (-1)**i * G->int_s1[i] */ + G->int_f1->coeffs[i] = tmp; /* (-1)**i * G->int_f2[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); } - // similarly, if either g1 or g2 have leading 0 coefficient, something is wrong +} + +void +nmod_geometric_progression_init(nmod_geometric_progression_t G, ulong r, slong len, nmod_t 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); + + G->len = len; + G->mod = mod; + + _nmod_geometric_progression_evaluate_init(G, r, len, mod, q, inv_r, inv_q); - _nmod_vec_clear(prod_diff); - _nmod_vec_clear(inv_diff); - _nmod_vec_clear(diff); + _nmod_geometric_progression_interpolate_init(G, len, mod, q, inv_q); } -void +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_vec_clear(G->ev_s); + nmod_poly_clear(G->ev_f); + _nmod_vec_clear(G->int_s1); + _nmod_vec_clear(G->int_s2); + _nmod_vec_clear(G->int_s3); + nmod_poly_clear(G->int_f1); + nmod_poly_clear(G->int_f2); } diff --git a/src/nmod_poly/interpolate_geometric_nmod_vec.c b/src/nmod_poly/interpolate_geometric_nmod_vec.c index 4d8a93028f..32368d6b5d 100644 --- a/src/nmod_poly/interpolate_geometric_nmod_vec.c +++ b/src/nmod_poly/interpolate_geometric_nmod_vec.c @@ -31,6 +31,8 @@ _nmod_poly_interpolate_geometric_nmod_vec_fast_precomp(nn_ptr poly, nn_srcptr v, return; } + /* TODO what is f? */ + /* TODO what is h? */ f = _nmod_vec_init(N); h = _nmod_vec_init(N); @@ -43,13 +45,13 @@ _nmod_poly_interpolate_geometric_nmod_vec_fast_precomp(nn_ptr poly, nn_srcptr v, } f1_len = N - i; - h1_len = FLINT_MIN(G->g1->length + f1_len - 1, N); + h1_len = FLINT_MIN(G->int_f1->length + f1_len - 1, N); for (i = 0; i < f1_len; i++) { - f[i] = nmod_mul(v[i], G->w[i], mod); + f[i] = nmod_mul(v[i], G->int_s1[i], mod); } - _nmod_poly_mullow(h, G->g1->coeffs, G->g1->length, f, f1_len, N, mod); + _nmod_poly_mullow(h, G->int_f1->coeffs, G->int_f1->length, f, f1_len, N, mod); while (h1_len > 0 && h[h1_len - 1] == 0) { @@ -67,16 +69,16 @@ _nmod_poly_interpolate_geometric_nmod_vec_fast_precomp(nn_ptr poly, nn_srcptr v, for (i = h2_min; i < h1_len; i++) { - f[N - 1 - i] = nmod_mul(h[i], G->y[i], mod); + f[N - 1 - i] = nmod_mul(h[i], G->int_s2[i], mod); } 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); + _nmod_poly_mullow(h, G->int_f2->coeffs, G->int_f2->length, f, f2_len, N, mod); for (i = 0; i < len; i++) { - poly[i] = nmod_mul(h[N - 1 - i], G->z[i], mod); + poly[i] = nmod_mul(h[N - 1 - i], G->int_s3[i], mod); } _nmod_vec_clear(f); 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..c1e42f3e2f --- /dev/null +++ b/src/nmod_poly/profile/p-eval_interp_precomputations.c @@ -0,0 +1,245 @@ +/* + 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 long 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; + +/* 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 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); +} + +/* geometric: fast, counting only precomputations */ +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(G, r, npoints_precomp, mod); + prof_stop(); + + nmod_geometric_progression_clear(G); + } + + 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}; + /* 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_fast_onlyprecomp; + double time_geom_fast_onlyprecomp; + + info_t info; + flint_bitcnt_t i; + + flint_printf("unit: measurements in ms\n"); + + for (i = 62; i <= FLINT_BITS; i++) + { + info.bits = i; + + printf("==== nbits = %ld====\n", i); + + if (func_bench == 0) /* bench all */ + flint_printf("len\tpoints | general | geometric \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 | geometric\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_fast_onlyprecomp, &tmp, sample_nmod_vec_fast_onlyprecomp, (void *) &info); + + if (func_bench == 0 || func_bench == 2) + prof_repeat(&time_geom_fast_onlyprecomp, &tmp, sample_nmod_vec_geom_fast_onlyprecomp, (void *) &info); + + if (func_bench == 0) + { + flint_printf("%ld\t%7ld| %.1e | %.1e\n", + info.length, info.npoints_precomp, + time_fast_onlyprecomp/fac, + time_geom_fast_onlyprecomp/fac); + } + else if (func_bench == 1) + { + flint_printf("%ld\t%7ld| %.1e\n", + info.length, info.npoints_precomp, + time_fast_onlyprecomp/fac); + } + else if (func_bench == 2) + { + flint_printf("%ld\t%7ld| %.1e\n", + info.length, info.npoints_precomp, + time_geom_fast_onlyprecomp/fac); + } + } + + flint_printf("\n"); + } + + return 0; +} + +#undef __NB_ITER From 74e5e2940f24261cdcf47e3ec7aba2f5688bdccd Mon Sep 17 00:00:00 2001 From: Vincent Neiger Date: Fri, 1 May 2026 18:10:03 +0200 Subject: [PATCH 02/13] adding Shoup variants to accelerate precomputations --- src/nmod_poly.h | 1 - src/nmod_poly/geometric_progression.c | 156 ++++++++++++++++-- src/nmod_poly/impl.h | 17 ++ .../profile/p-eval_interp_precomputations.c | 95 +++++++++-- 4 files changed, 234 insertions(+), 35 deletions(-) diff --git a/src/nmod_poly.h b/src/nmod_poly.h index 978d914a2e..957ae13660 100644 --- a/src/nmod_poly.h +++ b/src/nmod_poly.h @@ -68,7 +68,6 @@ typedef struct slong len; // number of points } nmod_geometric_progression_struct; - typedef nmod_geometric_progression_struct nmod_geometric_progression_t[1]; /* Memory management ********************************************************/ diff --git a/src/nmod_poly/geometric_progression.c b/src/nmod_poly/geometric_progression.c index fc874f1443..b4eac8b00c 100644 --- a/src/nmod_poly/geometric_progression.c +++ b/src/nmod_poly/geometric_progression.c @@ -10,12 +10,49 @@ (at your option) any later version. See . */ +#include "nmod.h" #include "nmod_vec.h" #include "nmod_poly.h" +#include "impl.h" -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) +/* specialized for moduli that support n_mulmod_shoup */ +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) +{ + /* 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; + + /* precomputations for Shoup multiplication */ + ulong r_pow_i2 = r; + ulong r_pr_quo = n_mulmod_precomp_shoup(r_pow_i2, mod.n); + + G->ev_f->coeffs[0] = 1; + for (slong i = 1; i < 2*len - 1; i++) + { + 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); + } + + r_pr_quo = n_mulmod_precomp_shoup(inv_r, mod.n); + G->ev_s[0] = 1; + for (slong i = 1; i < len; i++) + { + 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); + } +} + +/* general variant */ +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) */ @@ -27,7 +64,7 @@ _nmod_geometric_progression_evaluate_init(nmod_geometric_progression_t G, ulong for (slong i = 1; i < 2*len - 1; i++) { G->ev_f->coeffs[i] = nmod_mul(G->ev_f->coeffs[i - 1], r, mod); - r = nmod_mul(r, q, mod); + r = nmod_mul(r, q, mod); /* r**(2*i+1) */ } G->ev_s[0] = 1; @@ -38,9 +75,10 @@ _nmod_geometric_progression_evaluate_init(nmod_geometric_progression_t G, ulong } } -void -_nmod_geometric_progression_interpolate_init(nmod_geometric_progression_t G, slong len, nmod_t mod, - ulong q, ulong inv_q) +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] */ @@ -107,23 +145,109 @@ _nmod_geometric_progression_interpolate_init(nmod_geometric_progression_t G, slo } } -void -nmod_geometric_progression_init(nmod_geometric_progression_t G, ulong r, slong len, nmod_t mod) +void _nmod_geometric_progression_interpolate_init(nmod_geometric_progression_t G, + slong len, nmod_t mod, + ulong q, ulong inv_q) { - 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); + /* 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_f1, i) = (-1)**i * q_i / u_i */ + /* coeff(G->int_f2, i) = q_i / u_i, i.e. G->int_f2 == G->int_f1(-x) */ + nmod_poly_init2_preinv(G->int_f1, mod.n, mod.ninv, len); + nmod_poly_init2_preinv(G->int_f2, mod.n, mod.ninv, len); + G->int_f1->length = len; + G->int_f2->length = len; + /* G->int_s2[i] = (-1)**i * u_i / q_i = inverse of coeff(G->int_f1, i) */ + /* G->int_s3[i] = (-1)**i / u_i */ + /* G->int_s1[i] = 1 / u_i */ + G->int_s1 = _nmod_vec_init(len); + G->int_s2 = _nmod_vec_init(len); + G->int_s3 = _nmod_vec_init(len); + + G->int_f1->coeffs[0] = 1; + G->int_f2->coeffs[0] = 1; + G->int_s1[0] = 1; + G->int_s2[0] = 1; + G->int_s3[0] = 1; + + ulong q_pow_i = 1; + ulong inv_q_pow_i = 1; + ulong inv_q_i = 1; + ulong prod_diff = 1; + + /* TODO improve with Shoup */ + 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_f2->coeffs[i] = nmod_mul(G->int_f2->coeffs[i-1], q_pow_i, mod); /* q_i */ + q_pow_i = nmod_mul(q_pow_i, q, mod); /* q**i */ + G->int_s3[i] = q_pow_i - 1; /* temporarily, G->int_s3[i] = q**i - 1 */ + prod_diff = nmod_mul(q_pow_i - 1, prod_diff, mod); /* u_i */ + if (i % 2) /* i is odd */ + G->int_s2[i] = nmod_mul(mod.n - prod_diff, inv_q_i, mod); /* (-1)**i * u_i / q_i */ + else /* i is even */ + G->int_s2[i] = nmod_mul(prod_diff, inv_q_i, mod); /* (-1)**i * 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 */ + ulong tmp = nmod_mul(G->int_f2->coeffs[i], w_i, mod); + G->int_f2->coeffs[i] = tmp; /* q_i / u_i */ + G->int_s1[i-1] = nmod_mul(G->int_s3[i], w_i, mod); /* 1 / u_{i-1} */ + if (i % 2) /* i is odd */ + { + G->int_s3[i] = mod.n - w_i; /* (-1)**i * G->int_s1[i] */ + G->int_f1->coeffs[i] = mod.n - tmp; /* (-1)**i * G->int_f2[i] */ + } + else /* i is even */ + { + G->int_s3[i] = w_i; /* (-1)**i * G->int_s1[i] */ + G->int_f1->coeffs[i] = tmp; /* (-1)**i * G->int_f2[i] */ + } + } +} + +void nmod_geometric_progression_init(nmod_geometric_progression_t G, + ulong r, slong len, nmod_t mod) +{ G->len = len; G->mod = mod; - _nmod_geometric_progression_evaluate_init(G, r, len, mod, q, inv_r, inv_q); + 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); - _nmod_geometric_progression_interpolate_init(G, len, mod, q, inv_q); + 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); + + _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); + _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); + } + 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); + + _nmod_geometric_progression_evaluate_init(G, r, len, mod, q, inv_r, inv_q); + _nmod_geometric_progression_interpolate_init(G, len, mod, q, inv_q); + } } -void -nmod_geometric_progression_clear(nmod_geometric_progression_t G) +void nmod_geometric_progression_clear(nmod_geometric_progression_t G) { _nmod_vec_clear(G->ev_s); nmod_poly_clear(G->ev_f); diff --git a/src/nmod_poly/impl.h b/src/nmod_poly/impl.h index 80be7214a7..97dc5c0a9f 100644 --- a/src/nmod_poly/impl.h +++ b/src/nmod_poly/impl.h @@ -13,6 +13,7 @@ #define NMOD_POLY_IMPL_H #include "nmod_types.h" +#include "nmod_poly.h" /* TODO is this ok? */ void _nmod_poly_inv_series_basecase_preinv1(nn_ptr Qinv, nn_srcptr Q, slong Qlen, slong n, ulong q, nmod_t mod); void _nmod_poly_div_series_basecase_preinv1(nn_ptr Qinv, nn_srcptr P, slong Plen, nn_srcptr Q, slong Qlen, slong n, ulong q, nmod_t mod); @@ -22,4 +23,20 @@ int nmod_poly_irreducible_tetranomial(nmod_poly_t res, ulong n); int nmod_poly_irreducible_pentanomial(nmod_poly_t res, ulong n); void _nmod_poly_divrem_q0_preinv1(nn_ptr Q, nn_ptr R, nn_srcptr A, nn_srcptr B, slong lenA, ulong invL, nmod_t mod); +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); +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); +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); +void _nmod_geometric_progression_interpolate_init(nmod_geometric_progression_t G, slong len, nmod_t mod, + ulong q, ulong inv_q); + + #endif diff --git a/src/nmod_poly/profile/p-eval_interp_precomputations.c b/src/nmod_poly/profile/p-eval_interp_precomputations.c index c1e42f3e2f..6cc32403ce 100644 --- a/src/nmod_poly/profile/p-eval_interp_precomputations.c +++ b/src/nmod_poly/profile/p-eval_interp_precomputations.c @@ -16,6 +16,7 @@ #include "profiler.h" #include "ulong_extras.h" #include "nmod_poly.h" +#include "nmod_poly/impl.h" #include "nmod_vec.h" #define __NB_ITER 10 @@ -54,8 +55,8 @@ typedef struct slong npoints_precomp; } info_t; -/* tree, counting only precomputations */ -void sample_nmod_vec_fast_onlyprecomp(void * arg, ulong count) +/* precomputation for general points (subproduct tree) */ +void sample_general_precomp(void * arg, ulong count) { ulong n; nmod_t mod; @@ -91,8 +92,8 @@ void sample_nmod_vec_fast_onlyprecomp(void * arg, ulong count) FLINT_TEST_CLEAR(state); } -/* geometric: fast, counting only precomputations */ -void sample_nmod_vec_geom_fast_onlyprecomp(void * arg, ulong count) +/* precomputations for geometric progression (eval+interp) */ +void sample_geometric_precomp(void * arg, ulong count) { ulong n; nmod_t mod; @@ -116,10 +117,61 @@ void sample_nmod_vec_geom_fast_onlyprecomp(void * arg, ulong count) 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(); - nmod_geometric_progression_clear(G); + } + + 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; + 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); + + prof_start(); + if (NMOD_CAN_USE_SHOUP(mod) && npoints_precomp >= 4) + for (ulong ii = 0; ii < __NB_ITER; ii++) + { + _nmod_geometric_progression_evaluate_init_nonfullword(G, r, npoints_precomp, mod, q, inv_r, inv_q); + nmod_poly_clear(G->ev_f); + flint_free(G->ev_s); + } + else + for (ulong ii = 0; ii < __NB_ITER; ii++) + { + _nmod_geometric_progression_evaluate_init(G, r, npoints_precomp, mod, q, inv_r, inv_q); + nmod_poly_clear(G->ev_f); + flint_free(G->ev_s); + } + prof_stop(); } FLINT_TEST_CLEAR(state); @@ -181,26 +233,28 @@ int main(int argc, char * argv[]) /* }; */ double tmp; - double time_fast_onlyprecomp; - double time_geom_fast_onlyprecomp; + 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 = 62; i <= FLINT_BITS; i++) + 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 | geometric \n"); + flint_printf("len\tpoints | general | geom | g-eval \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 | geometric\n"); + flint_printf("len\tpoints | geom | g-eval \n"); for (int len = 0; len < nb_lens; ++len) { @@ -210,29 +264,34 @@ int main(int argc, char * argv[]) info.length = lengths[len]; if (func_bench == 0 || func_bench == 1) - prof_repeat(&time_fast_onlyprecomp, &tmp, sample_nmod_vec_fast_onlyprecomp, (void *) &info); + prof_repeat(&time_general, &tmp, sample_general_precomp, (void *) &info); if (func_bench == 0 || func_bench == 2) - prof_repeat(&time_geom_fast_onlyprecomp, &tmp, sample_nmod_vec_geom_fast_onlyprecomp, (void *) &info); + { + prof_repeat(&time_geometric_all, &tmp, sample_geometric_precomp, (void *) &info); + prof_repeat(&time_geometric_eval, &tmp, sample_geometric_precomp_eval, (void *) &info); + } if (func_bench == 0) { - flint_printf("%ld\t%7ld| %.1e | %.1e\n", + flint_printf("%ld\t%7ld| %.1e | %.1e | %.1e\n", info.length, info.npoints_precomp, - time_fast_onlyprecomp/fac, - time_geom_fast_onlyprecomp/fac); + time_general/fac, + time_geometric_all/fac, + time_geometric_eval/fac); } else if (func_bench == 1) { flint_printf("%ld\t%7ld| %.1e\n", info.length, info.npoints_precomp, - time_fast_onlyprecomp/fac); + time_general/fac); } else if (func_bench == 2) { - flint_printf("%ld\t%7ld| %.1e\n", + flint_printf("%ld\t%7ld| %.1e | %.1e\n", info.length, info.npoints_precomp, - time_geom_fast_onlyprecomp/fac); + time_geometric_all/fac, + time_geometric_eval/fac); } } From f50912127e63ebf7278dfc57347e57f4e0278ba4 Mon Sep 17 00:00:00 2001 From: Vincent Neiger Date: Fri, 1 May 2026 18:43:28 +0200 Subject: [PATCH 03/13] adding mask for selection of functionality + starting interpolation with Shoup --- src/nmod_poly.h | 8 +- src/nmod_poly/geometric_progression.c | 60 +++++++++++---- .../profile/p-eval_interp_precomputations.c | 77 +++++++++++++------ 3 files changed, 106 insertions(+), 39 deletions(-) diff --git a/src/nmod_poly.h b/src/nmod_poly.h index 957ae13660..a9d1cf1287 100644 --- a/src/nmod_poly.h +++ b/src/nmod_poly.h @@ -593,8 +593,14 @@ 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_mask(nmod_geometric_progression_t G, + ulong r, slong len, nmod_t mod, + ulong mask); +void nmod_geometric_progression_init(nmod_geometric_progression_t G, + ulong r, slong len, nmod_t mod); + +void _nmod_geometric_progression_clear_mask(nmod_geometric_progression_t G, ulong mask); 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); diff --git a/src/nmod_poly/geometric_progression.c b/src/nmod_poly/geometric_progression.c index b4eac8b00c..4751ab0ae4 100644 --- a/src/nmod_poly/geometric_progression.c +++ b/src/nmod_poly/geometric_progression.c @@ -214,8 +214,13 @@ void _nmod_geometric_progression_interpolate_init(nmod_geometric_progression_t G } } -void nmod_geometric_progression_init(nmod_geometric_progression_t G, - ulong r, slong len, nmod_t mod) +/* initialize for selection of functionalities: */ +/* the lowest 3 bits of `mask` act as a mask for the */ +/* three functionalities, in this order: */ +/* evaluate(bit0)+interpolate(bit1)+extrapolate(bit2) */ +void _nmod_geometric_progression_init_mask(nmod_geometric_progression_t G, + ulong r, slong len, nmod_t mod, + ulong mask) { G->len = len; G->mod = mod; @@ -233,8 +238,12 @@ void nmod_geometric_progression_init(nmod_geometric_progression_t G, 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); - _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); - _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 (mask & 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 ((mask>>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 ((mask>>2) & UWORD(1)) */ + /* extrapolate */ } else { @@ -242,18 +251,43 @@ void nmod_geometric_progression_init(nmod_geometric_progression_t G, const ulong inv_r = nmod_inv(r, mod); const ulong inv_q = nmod_mul(inv_r, inv_r, mod); - _nmod_geometric_progression_evaluate_init(G, r, len, mod, q, inv_r, inv_q); - _nmod_geometric_progression_interpolate_init(G, len, mod, q, inv_q); + if (mask & UWORD(1)) /* evaluate */ + _nmod_geometric_progression_evaluate_init(G, r, len, mod, q, inv_r, inv_q); + if ((mask>>1) & UWORD(1)) /* interpolate */ + _nmod_geometric_progression_interpolate_init(G, len, mod, q, inv_q); + /* if ((mask>>2) & UWORD(1)) */ + /* extrapolate */ + } +} + +/* clear for selection of functionalities (see init) */ +void _nmod_geometric_progression_clear_mask(nmod_geometric_progression_t G, ulong mask) +{ + if (mask & UWORD(1)) /* evaluate */ + { + _nmod_vec_clear(G->ev_s); + nmod_poly_clear(G->ev_f); } + if ((mask>>1) & UWORD(1)) /* interpolate */ + { + _nmod_vec_clear(G->int_s1); + _nmod_vec_clear(G->int_s2); + _nmod_vec_clear(G->int_s3); + nmod_poly_clear(G->int_f1); + nmod_poly_clear(G->int_f2); + } + /* if ((mask>>2) & UWORD(1)) */ +} + +/* 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_mask(G, r, len, mod, UWORD(7)); } +/* clear for all: evaluate+interpolate+extrapolate */ void nmod_geometric_progression_clear(nmod_geometric_progression_t G) { - _nmod_vec_clear(G->ev_s); - nmod_poly_clear(G->ev_f); - _nmod_vec_clear(G->int_s1); - _nmod_vec_clear(G->int_s2); - _nmod_vec_clear(G->int_s3); - nmod_poly_clear(G->int_f1); - nmod_poly_clear(G->int_f2); + _nmod_geometric_progression_clear_mask(G, UWORD(7)); } diff --git a/src/nmod_poly/profile/p-eval_interp_precomputations.c b/src/nmod_poly/profile/p-eval_interp_precomputations.c index 6cc32403ce..80e4879726 100644 --- a/src/nmod_poly/profile/p-eval_interp_precomputations.c +++ b/src/nmod_poly/profile/p-eval_interp_precomputations.c @@ -16,7 +16,6 @@ #include "profiler.h" #include "ulong_extras.h" #include "nmod_poly.h" -#include "nmod_poly/impl.h" #include "nmod_vec.h" #define __NB_ITER 10 @@ -152,25 +151,50 @@ void sample_geometric_precomp_eval(void * arg, ulong count) nmod_geometric_progression_t G; G->len = npoints_precomp; G->mod = 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); prof_start(); - if (NMOD_CAN_USE_SHOUP(mod) && npoints_precomp >= 4) - for (ulong ii = 0; ii < __NB_ITER; ii++) - { - _nmod_geometric_progression_evaluate_init_nonfullword(G, r, npoints_precomp, mod, q, inv_r, inv_q); - nmod_poly_clear(G->ev_f); - flint_free(G->ev_s); - } - else - for (ulong ii = 0; ii < __NB_ITER; ii++) - { - _nmod_geometric_progression_evaluate_init(G, r, npoints_precomp, mod, q, inv_r, inv_q); - nmod_poly_clear(G->ev_f); - flint_free(G->ev_s); - } + for (ulong ii = 0; ii < __NB_ITER; ii++) + { + _nmod_geometric_progression_init_mask(G, r, npoints_precomp, mod, UWORD(1)); + _nmod_geometric_progression_clear_mask(G, UWORD(1)); + } + 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_mask(G, r, npoints_precomp, mod, UWORD(2)); + _nmod_geometric_progression_clear_mask(G, UWORD(2)); + } prof_stop(); } @@ -236,7 +260,7 @@ int main(int argc, char * argv[]) double time_general; double time_geometric_all; double time_geometric_eval; - /* double time_geometric_interp; */ + double time_geometric_interp; info_t info; flint_bitcnt_t i; @@ -250,11 +274,11 @@ int main(int argc, char * argv[]) printf("==== nbits = %ld====\n", i); if (func_bench == 0) /* bench all */ - flint_printf("len\tpoints | general | geom | g-eval \n"); + 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 \n"); + flint_printf("len\tpoints | geom | g-eval | g-interp \n"); for (int len = 0; len < nb_lens; ++len) { @@ -270,15 +294,17 @@ int main(int argc, char * argv[]) { 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\n", + 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_eval/fac, + time_geometric_interp/fac); } else if (func_bench == 1) { @@ -288,10 +314,11 @@ int main(int argc, char * argv[]) } else if (func_bench == 2) { - flint_printf("%ld\t%7ld| %.1e | %.1e\n", + flint_printf("%ld\t%7ld| %.1e | %.1e | %.1e\n", info.length, info.npoints_precomp, time_geometric_all/fac, - time_geometric_eval/fac); + time_geometric_eval/fac, + time_geometric_interp/fac); } } From fc0fda1e48313df94c39478af50fdc2fc9e757b7 Mon Sep 17 00:00:00 2001 From: Vincent Neiger Date: Fri, 1 May 2026 22:33:13 +0200 Subject: [PATCH 04/13] add version of geometric_progression_init which exploit Shoup's modular multiplication with precomputation when the modulus has <= 63 bits --- src/nmod_poly/geometric_progression.c | 20 +++++++------ .../profile/p-eval_interp_precomputations.c | 28 ------------------- 2 files changed, 12 insertions(+), 36 deletions(-) diff --git a/src/nmod_poly/geometric_progression.c b/src/nmod_poly/geometric_progression.c index 4751ab0ae4..3773f5958e 100644 --- a/src/nmod_poly/geometric_progression.c +++ b/src/nmod_poly/geometric_progression.c @@ -105,24 +105,29 @@ void _nmod_geometric_progression_interpolate_init_nonfullword(nmod_geometric_pro G->int_s2[0] = 1; G->int_s3[0] = 1; + 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; - /* TODO improve with Shoup */ 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_f2->coeffs[i] = nmod_mul(G->int_f2->coeffs[i-1], q_pow_i, mod); /* q_i */ - q_pow_i = nmod_mul(q_pow_i, q, mod); /* q**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 */ + /* inv_q_pow_i = n_mulmod_shoup(inv_q, inv_q_pow_i, inv_q_pr_quo, mod.n); */ + 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_f2->coeffs[i] = n_mulmod_shoup(q_pow_i, G->int_f2->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_s3[i] = q_pow_i - 1; /* temporarily, G->int_s3[i] = q**i - 1 */ prod_diff = nmod_mul(q_pow_i - 1, prod_diff, mod); /* u_i */ if (i % 2) /* i is odd */ - G->int_s2[i] = nmod_mul(mod.n - prod_diff, inv_q_i, mod); /* (-1)**i * u_i / q_i */ + G->int_s2[i] = n_mulmod_shoup(inv_q_i, mod.n - prod_diff, inv_q_i_pr, mod.n); /* (-1)**i * u_i / q_i */ else /* i is even */ - G->int_s2[i] = nmod_mul(prod_diff, inv_q_i, mod); /* (-1)**i * u_i / q_i */ + G->int_s2[i] = n_mulmod_shoup(inv_q_i, prod_diff, inv_q_i_pr, mod.n); /* (-1)**i * u_i / q_i */ } G->int_s1[len-1] = nmod_inv(prod_diff, mod); /* 1 / u_{len-1} */ @@ -179,7 +184,6 @@ void _nmod_geometric_progression_interpolate_init(nmod_geometric_progression_t G ulong inv_q_i = 1; ulong prod_diff = 1; - /* TODO improve with Shoup */ for (slong i = 1; i < len; i++) { inv_q_i = nmod_mul(inv_q_i, inv_q_pow_i, mod); /* 1 / q_i */ diff --git a/src/nmod_poly/profile/p-eval_interp_precomputations.c b/src/nmod_poly/profile/p-eval_interp_precomputations.c index 80e4879726..087a4ebe51 100644 --- a/src/nmod_poly/profile/p-eval_interp_precomputations.c +++ b/src/nmod_poly/profile/p-eval_interp_precomputations.c @@ -227,34 +227,6 @@ int main(int argc, char * argv[]) 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_general; From 1a5bb66ae3a2429018f32a92ab9fdae49d32502b Mon Sep 17 00:00:00 2001 From: Vincent Neiger Date: Fri, 1 May 2026 23:51:37 +0200 Subject: [PATCH 05/13] improve handling of precomputation selection --- src/nmod_poly.h | 27 ++++----- src/nmod_poly/evaluate_geometric_nmod_vec.c | 2 +- src/nmod_poly/geometric_progression.c | 39 ++++++------- .../profile/p-eval_interp_precomputations.c | 8 +-- src/nmod_poly/profile/p-evaluate_nmod_vec.c | 57 +++++++++++++++++-- 5 files changed, 90 insertions(+), 43 deletions(-) diff --git a/src/nmod_poly.h b/src/nmod_poly.h index a9d1cf1287..523f80b615 100644 --- a/src/nmod_poly.h +++ b/src/nmod_poly.h @@ -64,8 +64,9 @@ typedef struct nmod_poly_t ev_f; /* evaluate: polynomial */ nn_ptr int_s1, int_s2, int_s3; /* interpolate: scaling constants */ nmod_poly_t int_f1, int_f2; /* interpolate: polynomials */ - nmod_t mod; - slong len; // number of points + 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]; @@ -464,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, @@ -594,13 +595,13 @@ void _nmod_poly_tree_build(nn_ptr * tree, nn_srcptr roots, slong len, nmod_t mod /* Geometric evaluation / interpolation *************************************/ -void _nmod_geometric_progression_init_mask(nmod_geometric_progression_t G, - ulong r, slong len, nmod_t mod, - ulong mask); +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_mask(nmod_geometric_progression_t G, ulong mask); +void _nmod_geometric_progression_clear_function(nmod_geometric_progression_t G, ulong function); 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); @@ -702,7 +703,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 a1c42f0f41..7a9e5d3c3c 100644 --- a/src/nmod_poly/evaluate_geometric_nmod_vec.c +++ b/src/nmod_poly/evaluate_geometric_nmod_vec.c @@ -129,7 +129,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 3773f5958e..3908038285 100644 --- a/src/nmod_poly/geometric_progression.c +++ b/src/nmod_poly/geometric_progression.c @@ -218,16 +218,17 @@ void _nmod_geometric_progression_interpolate_init(nmod_geometric_progression_t G } } -/* initialize for selection of functionalities: */ -/* the lowest 3 bits of `mask` act as a mask for the */ -/* three functionalities, in this order: */ -/* evaluate(bit0)+interpolate(bit1)+extrapolate(bit2) */ -void _nmod_geometric_progression_init_mask(nmod_geometric_progression_t G, - ulong r, slong len, nmod_t mod, - ulong mask) +/* 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)) { @@ -242,11 +243,11 @@ void _nmod_geometric_progression_init_mask(nmod_geometric_progression_t G, 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 (mask & UWORD(1)) /* evaluate */ + 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 ((mask>>1) & UWORD(1)) /* interpolate */ + 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 ((mask>>2) & UWORD(1)) */ + /* if ((function>>2) & UWORD(1)) */ /* extrapolate */ } else @@ -255,24 +256,24 @@ void _nmod_geometric_progression_init_mask(nmod_geometric_progression_t G, const ulong inv_r = nmod_inv(r, mod); const ulong inv_q = nmod_mul(inv_r, inv_r, mod); - if (mask & UWORD(1)) /* evaluate */ + if (function & UWORD(1)) /* evaluate */ _nmod_geometric_progression_evaluate_init(G, r, len, mod, q, inv_r, inv_q); - if ((mask>>1) & UWORD(1)) /* interpolate */ + if ((function>>1) & UWORD(1)) /* interpolate */ _nmod_geometric_progression_interpolate_init(G, len, mod, q, inv_q); - /* if ((mask>>2) & UWORD(1)) */ + /* if ((function>>2) & UWORD(1)) */ /* extrapolate */ } } /* clear for selection of functionalities (see init) */ -void _nmod_geometric_progression_clear_mask(nmod_geometric_progression_t G, ulong mask) +void _nmod_geometric_progression_clear_function(nmod_geometric_progression_t G, ulong function) { - if (mask & UWORD(1)) /* evaluate */ + if (function & UWORD(1)) /* evaluate */ { _nmod_vec_clear(G->ev_s); nmod_poly_clear(G->ev_f); } - if ((mask>>1) & UWORD(1)) /* interpolate */ + if ((function>>1) & UWORD(1)) /* interpolate */ { _nmod_vec_clear(G->int_s1); _nmod_vec_clear(G->int_s2); @@ -280,18 +281,18 @@ void _nmod_geometric_progression_clear_mask(nmod_geometric_progression_t G, ulon nmod_poly_clear(G->int_f1); nmod_poly_clear(G->int_f2); } - /* if ((mask>>2) & UWORD(1)) */ + /* if ((function>>2) & UWORD(1)) */ } /* 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_mask(G, r, len, mod, UWORD(7)); + _nmod_geometric_progression_init_function(G, r, len, mod, UWORD(7)); } /* clear for all: evaluate+interpolate+extrapolate */ void nmod_geometric_progression_clear(nmod_geometric_progression_t G) { - _nmod_geometric_progression_clear_mask(G, UWORD(7)); + _nmod_geometric_progression_clear_function(G, G->function); } diff --git a/src/nmod_poly/profile/p-eval_interp_precomputations.c b/src/nmod_poly/profile/p-eval_interp_precomputations.c index 087a4ebe51..e2c388dd70 100644 --- a/src/nmod_poly/profile/p-eval_interp_precomputations.c +++ b/src/nmod_poly/profile/p-eval_interp_precomputations.c @@ -155,8 +155,8 @@ void sample_geometric_precomp_eval(void * arg, ulong count) prof_start(); for (ulong ii = 0; ii < __NB_ITER; ii++) { - _nmod_geometric_progression_init_mask(G, r, npoints_precomp, mod, UWORD(1)); - _nmod_geometric_progression_clear_mask(G, UWORD(1)); + _nmod_geometric_progression_init_function(G, r, npoints_precomp, mod, UWORD(1)); + _nmod_geometric_progression_clear_function(G, UWORD(1)); } prof_stop(); } @@ -192,8 +192,8 @@ void sample_geometric_precomp_interp(void * arg, ulong count) prof_start(); for (ulong ii = 0; ii < __NB_ITER; ii++) { - _nmod_geometric_progression_init_mask(G, r, npoints_precomp, mod, UWORD(2)); - _nmod_geometric_progression_clear_mask(G, UWORD(2)); + _nmod_geometric_progression_init_function(G, r, npoints_precomp, mod, UWORD(2)); + _nmod_geometric_progression_clear_function(G, UWORD(2)); } prof_stop(); } diff --git a/src/nmod_poly/profile/p-evaluate_nmod_vec.c b/src/nmod_poly/profile/p-evaluate_nmod_vec.c index 9440bf0807..75e84238cf 100644 --- a/src/nmod_poly/profile/p-evaluate_nmod_vec.c +++ b/src/nmod_poly/profile/p-evaluate_nmod_vec.c @@ -367,6 +367,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)) @@ -439,6 +479,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; @@ -456,7 +497,7 @@ 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("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 */ @@ -468,7 +509,7 @@ int main(int argc, char * argv[]) 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("len\tpoints | GEOMETRIC PROGRESSION | POLY_MUL\n"); flint_printf("len\tpoints |iter\tfast\tw/ prec\tprecomp\n"); } @@ -510,14 +551,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 +571,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); } } } From 6282e7049ff8dad3ef8ab3c47525107eca2a401c Mon Sep 17 00:00:00 2001 From: Vincent Neiger Date: Sat, 2 May 2026 11:28:50 +0200 Subject: [PATCH 06/13] adding profile for interpolate_nmod_vec (general points + geometric progression), and minor enhancements --- src/nmod_poly/evaluate_geometric_nmod_vec.c | 1 + .../interpolate_geometric_nmod_vec.c | 26 +- src/nmod_poly/profile/p-evaluate_nmod_vec.c | 22 +- .../profile/p-interpolate_nmod_vec.c | 592 ++++++++++++++++++ 4 files changed, 619 insertions(+), 22 deletions(-) create mode 100644 src/nmod_poly/profile/p-interpolate_nmod_vec.c diff --git a/src/nmod_poly/evaluate_geometric_nmod_vec.c b/src/nmod_poly/evaluate_geometric_nmod_vec.c index 7a9e5d3c3c..c4a85768fe 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) diff --git a/src/nmod_poly/interpolate_geometric_nmod_vec.c b/src/nmod_poly/interpolate_geometric_nmod_vec.c index 32368d6b5d..554ea5c278 100644 --- a/src/nmod_poly/interpolate_geometric_nmod_vec.c +++ b/src/nmod_poly/interpolate_geometric_nmod_vec.c @@ -14,17 +14,16 @@ #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) { + FLINT_ASSERT((G->function >> 1) & 1); + slong i, N, f1_len, f2_len, h1_len, h2_min; nn_ptr f, h; N = G->len; - FLINT_ASSERT(len == N); - if (len == 1) { poly[0] = v[0]; @@ -85,23 +84,22 @@ _nmod_poly_interpolate_geometric_nmod_vec_fast_precomp(nn_ptr poly, nn_srcptr v, _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 n) { nmod_geometric_progression_t G; - nmod_geometric_progression_init(G, r, n, poly->mod); - + _nmod_geometric_progression_init_function(G, r, n, poly->mod, UWORD(2)); nmod_poly_interpolate_geometric_nmod_vec_fast_precomp(poly, ys, G, n); - nmod_geometric_progression_clear(G); + _nmod_geometric_progression_clear_function(G, UWORD(2)); } diff --git a/src/nmod_poly/profile/p-evaluate_nmod_vec.c b/src/nmod_poly/profile/p-evaluate_nmod_vec.c index 75e84238cf..63c228d166 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 @@ -430,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 `pre_fact * points` points\n"); + 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; @@ -496,19 +499,22 @@ 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("==== 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("==== 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"); } 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..ad806dc2a3 --- /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 long 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 From 0c91e56867e7fc19616093dc2244e03f0d266741 Mon Sep 17 00:00:00 2001 From: Vincent Neiger Date: Sat, 2 May 2026 16:33:18 +0200 Subject: [PATCH 07/13] augment geometric interpolate test to cover cases where the number of points in precomputation differs from the number of values to interpolate --- .../interpolate_geometric_nmod_vec.c | 8 ++ src/nmod_poly/profile/p-evaluate_nmod_vec.c | 2 +- .../t-interpolate_geometric_nmod_vec_fast.c | 97 +++++++++++++++---- 3 files changed, 85 insertions(+), 22 deletions(-) diff --git a/src/nmod_poly/interpolate_geometric_nmod_vec.c b/src/nmod_poly/interpolate_geometric_nmod_vec.c index 554ea5c278..c6787cb80b 100644 --- a/src/nmod_poly/interpolate_geometric_nmod_vec.c +++ b/src/nmod_poly/interpolate_geometric_nmod_vec.c @@ -35,6 +35,9 @@ void _nmod_poly_interpolate_geometric_nmod_vec_fast_precomp(nn_ptr poly, f = _nmod_vec_init(N); h = _nmod_vec_init(N); + /* step1: Newton interpolation */ + /* [Bostan - Schost, J.Complexity 2005, Section 5.1] */ + for (i = 0; i < N; i++) { if (v[N - i - 1] != 0) @@ -52,6 +55,11 @@ void _nmod_poly_interpolate_geometric_nmod_vec_fast_precomp(nn_ptr poly, } _nmod_poly_mullow(h, G->int_f1->coeffs, G->int_f1->length, f, f1_len, N, mod); + /* FIXME division by q_i ?? */ + + /* step2: Newton basis -> monomial basis */ + /* [Bostan - Schost, J.Complexity 2005, Section 5.2] */ + while (h1_len > 0 && h[h1_len - 1] == 0) { h1_len--; diff --git a/src/nmod_poly/profile/p-evaluate_nmod_vec.c b/src/nmod_poly/profile/p-evaluate_nmod_vec.c index 63c228d166..a63d77db82 100644 --- a/src/nmod_poly/profile/p-evaluate_nmod_vec.c +++ b/src/nmod_poly/profile/p-evaluate_nmod_vec.c @@ -433,7 +433,7 @@ int main(int argc, char * argv[]) if (pre_fact != 1) { - flint_printf("[pre_fact provided] -> will build geometric progression with `pre_fact * points` points\n"); + 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"); } 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..dca0cabffd 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,96 @@ 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)); + 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 (`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); } - nmod_poly_clear(P); - nmod_poly_clear(Q); - _nmod_vec_clear(y); + flint_free(y); } TEST_FUNCTION_END(state); From 4a6928cd301964cbc8f5c46b86d2d7780ccabe89 Mon Sep 17 00:00:00 2001 From: Vincent Neiger Date: Sat, 2 May 2026 17:35:59 +0200 Subject: [PATCH 08/13] first step of geometric interpolate: supporter smaller lengths second step of geometric interpolate: explain computation --- .../interpolate_geometric_nmod_vec.c | 120 +++++++++++------- .../t-interpolate_geometric_nmod_vec_fast.c | 2 +- 2 files changed, 78 insertions(+), 44 deletions(-) diff --git a/src/nmod_poly/interpolate_geometric_nmod_vec.c b/src/nmod_poly/interpolate_geometric_nmod_vec.c index c6787cb80b..4e2c60f193 100644 --- a/src/nmod_poly/interpolate_geometric_nmod_vec.c +++ b/src/nmod_poly/interpolate_geometric_nmod_vec.c @@ -17,9 +17,10 @@ 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) { + FLINT_ASSERT(len <= G->len); FLINT_ASSERT((G->function >> 1) & 1); - slong i, N, f1_len, f2_len, h1_len, h2_min; + slong i, N, f_len, h_len; nn_ptr f, h; N = G->len; @@ -30,63 +31,96 @@ void _nmod_poly_interpolate_geometric_nmod_vec_fast_precomp(nn_ptr poly, return; } - /* TODO what is f? */ - /* TODO what is h? */ f = _nmod_vec_init(N); h = _nmod_vec_init(N); - /* step1: Newton interpolation */ - /* [Bostan - Schost, J.Complexity 2005, Section 5.1] */ - - 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_f1 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 = 0; val < len; val++) + if (v[val] != 0) break; - } - } - - f1_len = N - i; - h1_len = FLINT_MIN(G->int_f1->length + f1_len - 1, N); - - for (i = 0; i < f1_len; i++) + if (val == len) { - f[i] = nmod_mul(v[i], G->int_s1[i], mod); + _nmod_vec_zero(poly, len); + return; } - _nmod_poly_mullow(h, G->int_f1->coeffs, G->int_f1->length, f, f1_len, N, mod); - - /* FIXME division by q_i ?? */ - /* step2: Newton basis -> monomial basis */ - /* [Bostan - Schost, J.Complexity 2005, Section 5.2] */ + /* i = number of trailing zero coefficients */ + 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 (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_f1 mod x**len */ + /* == x**val (f * G->int_f1 mod x**(len-val)) */ + /* note: len - val is <= G->int_f1->length, since G->intf_1 has */ + /* length G->len >= len (all its coefficients are nonzero) */ + _nmod_vec_zero(h, val); + _nmod_poly_mullow(h+val, G->int_f1->coeffs, len - val, f, f_len, len - val, mod); + + for (h_len = len; h_len > val; h_len--) + if (h[h_len-1] != 0) + break; - while (h1_len > 0 && h[h1_len - 1] == 0) - { - h1_len--; - } + /* FIXME division by q_i ?? */ - for (i = 0; i < h1_len; i++) - { - if (h[i] != 0) - { + /* TODO remove debug */ + /* if (f_len != len) */ + /* flint_printf("len = %ld, G->len = %ld, flen = %ld\n", len, G->len, f_len); */ + + /** step2: Newton basis -> monomial basis + * [Bostan - Schost, J.Complexity 2005, Section 5.2] + * in the Newton basis, poly has coefficients h[i] / q_i, i=0..len-1 + * -> Convert it to monomial basis, through the transposed product + * mul_t(len-1, sum_{i=0}^{len-1} uu_i x**i, + * 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 - q**k) + * == (-1)**i * q_i / u_i, for u_i as above + * This gives `len` coefficients and it just remains to scale the i-th one + * by (-1)**i * uu_i / q_i == 1 / u_i + * -> so, with the precomputed data, the two polynomials have coefficients + * uu_i == G->int_f1[i] and (-1)**i * h[i]/uu_i == (-1)**i * h[i] * G->int_s2[i] + * and the final scaling factors are 1/u_i == G->int_s1[i] + * meaning that we are looking to compute + * mul_t(len-1, G->int_f1, f) + * where f = sum_{i=0}^{len-1} h[i] * G->int_s2[i] x**i + * -> from that, it remains to scale coefficients by 1/u_i == G->int_s1[i] + */ + + for (val = 0; val < h_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->int_s2[i], mod); - } + for (i = val; i < h_len; i++) + f[N - 1 - i] = nmod_mul(h[i], G->int_s2[i], mod); - f2_len = N - h2_min; - _nmod_vec_zero(f, N - h1_len); - _nmod_poly_mullow(h, G->int_f2->coeffs, G->int_f2->length, f, f2_len, N, mod); + f_len = N - val; + _nmod_vec_zero(f, N - h_len); + _nmod_poly_mullow(h, G->int_f2->coeffs, G->int_f2->length, f, f_len, N, mod); for (i = 0; i < len; i++) - { poly[i] = nmod_mul(h[N - 1 - i], G->int_s3[i], mod); - } _nmod_vec_clear(f); _nmod_vec_clear(h); 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 dca0cabffd..de470ae7c0 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 @@ -100,7 +100,7 @@ TEST_FUNCTION_START(nmod_poly_interpolate_geometric_nmod_vec_fast, state) result = nmod_poly_equal(P, Q); if (!result) { - flint_printf("FAIL (`len` points):\n"); + 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"); From b207261888f7e94d046817b524f2f82d5710e72f Mon Sep 17 00:00:00 2001 From: Vincent Neiger Date: Sat, 2 May 2026 22:19:20 +0200 Subject: [PATCH 09/13] support smaller lengths in second step of geometric interpolate; simplify precomputations by removing data that is not used anymore --- src/nmod_poly.h | 14 +-- src/nmod_poly/geometric_progression.c | 104 ++++++------------ .../interpolate_geometric_nmod_vec.c | 77 +++++++------ .../t-interpolate_geometric_nmod_vec_fast.c | 3 +- 4 files changed, 83 insertions(+), 115 deletions(-) diff --git a/src/nmod_poly.h b/src/nmod_poly.h index 523f80b615..a90be1f599 100644 --- a/src/nmod_poly.h +++ b/src/nmod_poly.h @@ -60,13 +60,13 @@ nmod_poly_compose_mod_precomp_preinv_arg_t; typedef struct { - nn_ptr ev_s; /* evaluate: scaling constants */ - nmod_poly_t ev_f; /* evaluate: polynomial */ - nn_ptr int_s1, int_s2, int_s3; /* interpolate: scaling constants */ - nmod_poly_t int_f1, int_f2; /* interpolate: polynomials */ - nmod_t mod; /* modulus */ - slong len; /* number of points */ - ulong function; /* choice of precomputations */ + 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]; diff --git a/src/nmod_poly/geometric_progression.c b/src/nmod_poly/geometric_progression.c index 3908038285..82ce6934cd 100644 --- a/src/nmod_poly/geometric_progression.c +++ b/src/nmod_poly/geometric_progression.c @@ -85,25 +85,17 @@ void _nmod_geometric_progression_interpolate_init_nonfullword(nmod_geometric_pro /* write u_i for prod_{1 <= k <= i} (q**k - 1), */ /* and q_i for q**(i*(i-1)/2) */ - /* coeff(G->int_f1, i) = (-1)**i * q_i / u_i */ - /* coeff(G->int_f2, i) = q_i / u_i, i.e. G->int_f2 == G->int_f1(-x) */ - nmod_poly_init2_preinv(G->int_f1, mod.n, mod.ninv, len); - nmod_poly_init2_preinv(G->int_f2, mod.n, mod.ninv, len); - G->int_f1->length = len; - G->int_f2->length = len; - - /* G->int_s2[i] = (-1)**i * u_i / q_i = inverse of coeff(G->int_f1, i) */ - /* G->int_s3[i] = (-1)**i / u_i */ + /* 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_s3 = _nmod_vec_init(len); - G->int_f1->coeffs[0] = 1; - G->int_f2->coeffs[0] = 1; - G->int_s1[0] = 1; + G->int_f->coeffs[0] = 1; G->int_s2[0] = 1; - G->int_s3[0] = 1; const ulong one_precomp = n_mulmod_precomp_shoup(UWORD(1), mod.n); ulong q_pow_i = 1; @@ -118,35 +110,23 @@ void _nmod_geometric_progression_interpolate_init_nonfullword(nmod_geometric_pro { 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 */ - /* inv_q_pow_i = n_mulmod_shoup(inv_q, inv_q_pow_i, inv_q_pr_quo, mod.n); */ 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_f2->coeffs[i] = n_mulmod_shoup(q_pow_i, G->int_f2->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_s3[i] = q_pow_i - 1; /* temporarily, G->int_s3[i] = q**i - 1 */ + 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 */ - if (i % 2) /* i is odd */ - G->int_s2[i] = n_mulmod_shoup(inv_q_i, mod.n - prod_diff, inv_q_i_pr, mod.n); /* (-1)**i * u_i / q_i */ - else /* i is even */ - G->int_s2[i] = n_mulmod_shoup(inv_q_i, prod_diff, inv_q_i_pr, mod.n); /* (-1)**i * u_i / q_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--) { - ulong w_i = G->int_s1[i]; /* 1 / u_i */ - ulong tmp = nmod_mul(G->int_f2->coeffs[i], w_i, mod); - G->int_f2->coeffs[i] = tmp; /* q_i / u_i */ - G->int_s1[i-1] = nmod_mul(G->int_s3[i], w_i, mod); /* 1 / u_{i-1} */ - if (i % 2) /* i is odd */ - { - G->int_s3[i] = mod.n - w_i; /* (-1)**i * G->int_s1[i] */ - G->int_f1->coeffs[i] = mod.n - tmp; /* (-1)**i * G->int_f2[i] */ - } - else /* i is even */ - { - G->int_s3[i] = w_i; /* (-1)**i * G->int_s1[i] */ - G->int_f1->coeffs[i] = tmp; /* (-1)**i * G->int_f2[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} */ } } @@ -159,25 +139,17 @@ void _nmod_geometric_progression_interpolate_init(nmod_geometric_progression_t G /* write u_i for prod_{1 <= k <= i} (q**k - 1), */ /* and q_i for q**(i*(i-1)/2) */ - /* coeff(G->int_f1, i) = (-1)**i * q_i / u_i */ - /* coeff(G->int_f2, i) = q_i / u_i, i.e. G->int_f2 == G->int_f1(-x) */ - nmod_poly_init2_preinv(G->int_f1, mod.n, mod.ninv, len); - nmod_poly_init2_preinv(G->int_f2, mod.n, mod.ninv, len); - G->int_f1->length = len; - G->int_f2->length = len; + /* 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_s2[i] = (-1)**i * u_i / q_i = inverse of coeff(G->int_f1, i) */ - /* G->int_s3[i] = (-1)**i / u_i */ /* 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_s3 = _nmod_vec_init(len); - G->int_f1->coeffs[0] = 1; - G->int_f2->coeffs[0] = 1; - G->int_s1[0] = 1; + G->int_f->coeffs[0] = 1; G->int_s2[0] = 1; - G->int_s3[0] = 1; ulong q_pow_i = 1; ulong inv_q_pow_i = 1; @@ -188,33 +160,22 @@ void _nmod_geometric_progression_interpolate_init(nmod_geometric_progression_t G { 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_f2->coeffs[i] = nmod_mul(G->int_f2->coeffs[i-1], q_pow_i, mod); /* 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_s3[i] = q_pow_i - 1; /* temporarily, G->int_s3[i] = q**i - 1 */ + 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 */ - if (i % 2) /* i is odd */ - G->int_s2[i] = nmod_mul(mod.n - prod_diff, inv_q_i, mod); /* (-1)**i * u_i / q_i */ - else /* i is even */ - G->int_s2[i] = nmod_mul(prod_diff, inv_q_i, mod); /* (-1)**i * u_i / q_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 */ - ulong tmp = nmod_mul(G->int_f2->coeffs[i], w_i, mod); - G->int_f2->coeffs[i] = tmp; /* q_i / u_i */ - G->int_s1[i-1] = nmod_mul(G->int_s3[i], w_i, mod); /* 1 / u_{i-1} */ - if (i % 2) /* i is odd */ - { - G->int_s3[i] = mod.n - w_i; /* (-1)**i * G->int_s1[i] */ - G->int_f1->coeffs[i] = mod.n - tmp; /* (-1)**i * G->int_f2[i] */ - } - else /* i is even */ - { - G->int_s3[i] = w_i; /* (-1)**i * G->int_s1[i] */ - G->int_f1->coeffs[i] = tmp; /* (-1)**i * G->int_f2[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} */ } } @@ -277,11 +238,10 @@ void _nmod_geometric_progression_clear_function(nmod_geometric_progression_t G, { _nmod_vec_clear(G->int_s1); _nmod_vec_clear(G->int_s2); - _nmod_vec_clear(G->int_s3); - nmod_poly_clear(G->int_f1); - nmod_poly_clear(G->int_f2); + nmod_poly_clear(G->int_f); } /* if ((function>>2) & UWORD(1)) */ + /* extrapolate */ } /* initialize for all: evaluate+interpolate+extrapolate */ diff --git a/src/nmod_poly/interpolate_geometric_nmod_vec.c b/src/nmod_poly/interpolate_geometric_nmod_vec.c index 4e2c60f193..181632ccda 100644 --- a/src/nmod_poly/interpolate_geometric_nmod_vec.c +++ b/src/nmod_poly/interpolate_geometric_nmod_vec.c @@ -45,7 +45,7 @@ void _nmod_poly_interpolate_geometric_nmod_vec_fast_precomp(nn_ptr poly, * 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_f1 mod x**len + * f * G->int_f mod x**len * where f = sum_{i=0}^{len-1} v[i] * G->int_s1[i] x**i */ @@ -60,7 +60,7 @@ void _nmod_poly_interpolate_geometric_nmod_vec_fast_precomp(nn_ptr poly, return; } - /* i = number of trailing zero coefficients */ + /* actual length of f */ for (f_len = len; f_len > val; f_len--) if (v[f_len-1] != 0) break; @@ -71,56 +71,63 @@ void _nmod_poly_interpolate_geometric_nmod_vec_fast_precomp(nn_ptr poly, for (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_f1 mod x**len */ - /* == x**val (f * G->int_f1 mod x**(len-val)) */ - /* note: len - val is <= G->int_f1->length, since G->intf_1 has */ + /* 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_vec_zero(h, val); - _nmod_poly_mullow(h+val, G->int_f1->coeffs, len - val, f, f_len, len - val, mod); + _nmod_poly_mullow(h+val, G->int_f->coeffs, len - val, f, f_len, len - val, mod); - for (h_len = len; h_len > val; h_len--) - if (h[h_len-1] != 0) - break; - - /* FIXME division by q_i ?? */ - - /* TODO remove debug */ - /* if (f_len != len) */ - /* flint_printf("len = %ld, G->len = %ld, flen = %ld\n", len, G->len, f_len); */ + /* for Newton interpolation, here we should compute h[i] = h[i]/q_i */ + /* yet this will simplify just below, so we just leave h as it is */ - /** step2: Newton basis -> monomial basis + /** step2: Newton basis -> monomial basis * [Bostan - Schost, J.Complexity 2005, Section 5.2] - * in the Newton basis, poly has coefficients h[i] / q_i, i=0..len-1 - * -> Convert it to monomial basis, through the transposed product - * mul_t(len-1, sum_{i=0}^{len-1} uu_i x**i, - * sum_{i=0}^{len-1} (-1)**i * (h[i]/q_i)*q_i/uu_i x**i) + * -> 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 - q**k) + * and uu_i = prod_{1 <= k <= i} q**(k-1) / (1 - q**k) * == (-1)**i * q_i / u_i, for u_i as above - * This gives `len` coefficients and it just remains to scale the i-th one + * (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 - * -> so, with the precomputed data, the two polynomials have coefficients - * uu_i == G->int_f1[i] and (-1)**i * h[i]/uu_i == (-1)**i * h[i] * G->int_s2[i] - * and the final scaling factors are 1/u_i == G->int_s1[i] - * meaning that we are looking to compute - * mul_t(len-1, G->int_f1, f) - * where f = sum_{i=0}^{len-1} h[i] * G->int_s2[i] x**i - * -> from that, it remains to scale coefficients by 1/u_i == G->int_s1[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] */ - for (val = 0; val < h_len; val++) + /* valuation of h is at least val, see if it is higher */ + for (; val < len; val++) if (h[val] != 0) break; + /* actual length of h */ + for (h_len = len; h_len > val; h_len--) + if (h[h_len-1] != 0) + break; + + /* compute reversed and scaled f */ + _nmod_vec_zero(f, len - h_len); /* TODO necessary? */ for (i = val; i < h_len; i++) - f[N - 1 - i] = nmod_mul(h[i], G->int_s2[i], mod); + f[len-1-i] = nmod_mul(h[i], G->int_s2[i], mod); + _nmod_vec_zero(f+len-val, val); /* TODO necessary? */ - f_len = N - val; - _nmod_vec_zero(f, N - h_len); - _nmod_poly_mullow(h, G->int_f2->coeffs, G->int_f2->length, f, f_len, N, mod); + /* transposed short product */ + _nmod_poly_mullow(h, f, len, G->int_f->coeffs, len, len, mod); + /* final scaling */ for (i = 0; i < len; i++) - poly[i] = nmod_mul(h[N - 1 - i], G->int_s3[i], mod); + poly[i] = nmod_mul(h[len-1-i], G->int_s1[i], mod); _nmod_vec_clear(f); _nmod_vec_clear(h); 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 de470ae7c0..561c82568f 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 @@ -25,7 +25,8 @@ TEST_FUNCTION_START(nmod_poly_interpolate_geometric_nmod_vec_fast, state) ulong mod, r; slong len, npoints; - npoints = 1 + (i == 0 ? 1 : n_randint(state, 200)); + /* npoints = 1 + (i == 0 ? 1 : n_randint(state, 200)); */ + npoints = 1 + n_randint(state, 200); len = 1 + n_randint(state, npoints); do { From cdc1d3ae9a9809f3d46940cb9639fd8592666dc5 Mon Sep 17 00:00:00 2001 From: Vincent Neiger Date: Sat, 2 May 2026 23:15:48 +0200 Subject: [PATCH 10/13] improve mullow in interpolate when valuation improve mulmid in evaluate --- .../interpolate_geometric_nmod_vec.c | 28 ++++++++----------- 1 file changed, 11 insertions(+), 17 deletions(-) diff --git a/src/nmod_poly/interpolate_geometric_nmod_vec.c b/src/nmod_poly/interpolate_geometric_nmod_vec.c index 181632ccda..f4e791bfc6 100644 --- a/src/nmod_poly/interpolate_geometric_nmod_vec.c +++ b/src/nmod_poly/interpolate_geometric_nmod_vec.c @@ -20,19 +20,15 @@ void _nmod_poly_interpolate_geometric_nmod_vec_fast_precomp(nn_ptr poly, FLINT_ASSERT(len <= G->len); FLINT_ASSERT((G->function >> 1) & 1); - slong i, N, f_len, h_len; - nn_ptr f, h; - - N = G->len; - if (len == 1) { poly[0] = v[0]; return; } - f = _nmod_vec_init(N); - h = _nmod_vec_init(N); + slong f_len, h_len, val; + nn_ptr f = _nmod_vec_init(len); + nn_ptr h = _nmod_vec_init(len); /** step1: Newton interpolation * [Bostan - Schost, J.Complexity 2005, Section 5.1] @@ -50,7 +46,6 @@ void _nmod_poly_interpolate_geometric_nmod_vec_fast_precomp(nn_ptr poly, */ /* val = valuation of output poly in Newton basis */ - slong val = 0; for (val = 0; val < len; val++) if (v[val] != 0) break; @@ -68,18 +63,18 @@ void _nmod_poly_interpolate_geometric_nmod_vec_fast_precomp(nn_ptr poly, /* 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 (i = 0; i < f_len; 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_vec_zero(h, val); _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 will simplify just below, so we just leave h as it is */ + /* 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] @@ -117,16 +112,15 @@ void _nmod_poly_interpolate_geometric_nmod_vec_fast_precomp(nn_ptr poly, break; /* compute reversed and scaled f */ - _nmod_vec_zero(f, len - h_len); /* TODO necessary? */ - for (i = val; i < h_len; i++) - f[len-1-i] = nmod_mul(h[i], G->int_s2[i], mod); - _nmod_vec_zero(f+len-val, val); /* TODO necessary? */ + 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); /* transposed short product */ - _nmod_poly_mullow(h, f, len, G->int_f->coeffs, len, len, mod); + _nmod_poly_mullow(h+len-h_len, f, h_len-val, G->int_f->coeffs, h_len, h_len, mod); /* final scaling */ - for (i = 0; i < len; i++) + _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); From d7ba9d1214dc5b33a5bf366ffc004e8caa0ea6bb Mon Sep 17 00:00:00 2001 From: Vincent Neiger Date: Sun, 3 May 2026 00:40:56 +0200 Subject: [PATCH 11/13] simplify geometric evaluate following improvements of nmod_poly_mulmid --- src/nmod_poly/evaluate_geometric_nmod_vec.c | 62 ++++++------------- src/nmod_poly/profile/p-evaluate_nmod_vec.c | 36 +++++------ .../t-interpolate_geometric_nmod_vec_fast.c | 1 + 3 files changed, 38 insertions(+), 61 deletions(-) diff --git a/src/nmod_poly/evaluate_geometric_nmod_vec.c b/src/nmod_poly/evaluate_geometric_nmod_vec.c index c4a85768fe..b6738c74a0 100644 --- a/src/nmod_poly/evaluate_geometric_nmod_vec.c +++ b/src/nmod_poly/evaluate_geometric_nmod_vec.c @@ -52,47 +52,19 @@ void _nmod_poly_evaluate_geometric_nmod_vec_fast_precomp(nn_ptr vs, nn_srcptr po * [rev(a) * (G->ev_f >> x**val)]_{alen - 1}^{len} (that is, coeffs [alen - 1, alen - 1 + len)) */ - /* below are 3 different versions (one requires fft_small) */ - /* TODO once some optimized middle product is written, the "version 1" should probably be discarded */ + /* 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->ev_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 - **/ -#if FLINT_HAVE_FFT_SMALL - if (2 * (plen - val) - 2 + len <= 192) -#else - if (1) -#endif - { - slong blen = alen - 1 + len; - nn_ptr a = _nmod_vec_init(alen); - nn_ptr b = _nmod_vec_init(blen); - - for (slong i = val; i < plen; i++) - a[i - val] = nmod_mul(G->ev_s[i], poly[i], mod); - - nn_ptr Frev = _nmod_vec_init(blen); - _nmod_poly_reverse(Frev, G->ev_f->coeffs + val, blen, blen); - _nmod_poly_mullow(b, Frev, blen, a, alen, blen, mod); - - for (slong i = 0; i < len; i++) - vs[i] = nmod_mul(G->ev_s[i], b[alen - 1 + len - 1 - i], mod); - - _nmod_vec_clear(Frev); - _nmod_vec_clear(a); - _nmod_vec_clear(b); - } - else - { - /* version 2 */ /* 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 - /* version 2.a uses fft_small directly (2025-12-04: fastest in medium and large lengths, like 100 and more) */ + if (2 * (plen - val) - 2 + len > 192) + { + /* 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++) @@ -104,23 +76,27 @@ void _nmod_poly_evaluate_geometric_nmod_vec_fast_precomp(nn_ptr vs, nn_srcptr po vs[i] = nmod_mul(G->ev_s[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) */ + } + else +#endif + { + /* 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->ev_s[i], poly[i], mod); - _nmod_poly_mulhigh(b, G->ev_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->ev_s[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 } } diff --git a/src/nmod_poly/profile/p-evaluate_nmod_vec.c b/src/nmod_poly/profile/p-evaluate_nmod_vec.c index a63d77db82..ae63472f5b 100644 --- a/src/nmod_poly/profile/p-evaluate_nmod_vec.c +++ b/src/nmod_poly/profile/p-evaluate_nmod_vec.c @@ -447,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; 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 561c82568f..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 @@ -111,6 +111,7 @@ TEST_FUNCTION_START(nmod_poly_interpolate_geometric_nmod_vec_fast, state) nmod_poly_clear(P); nmod_poly_clear(Q); + nmod_geometric_progression_clear(G); } flint_free(y); From 135090416c7bc1dd9a82555f88eadc6effd574a8 Mon Sep 17 00:00:00 2001 From: Vincent Neiger Date: Sun, 3 May 2026 13:55:18 +0200 Subject: [PATCH 12/13] - update documentation after modifications of geometric progression and interpolation - documentation contained sometimes `len`, sometimes `n` for the number of points, with `n` also used for the modulus value in close-by functions: unify to `len` for the number of points; or `ilen` (input length) and `olen` (output length) when necessary - some cleaning in impl.h, making these functions static - minor fix in the profile and interpolation function to avoid some memory leak --- doc/source/nmod_poly.rst | 123 +++++++++--------- src/nmod_poly.h | 40 +++--- src/nmod_poly/geometric_progression.c | 8 +- src/nmod_poly/impl.h | 17 --- .../interpolate_geometric_nmod_vec.c | 19 +-- .../profile/p-eval_interp_precomputations.c | 4 +- 6 files changed, 102 insertions(+), 109 deletions(-) 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 a90be1f599..29b0e74801 100644 --- a/src/nmod_poly.h +++ b/src/nmod_poly.h @@ -554,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); @@ -594,40 +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_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_function(nmod_geometric_progression_t G, ulong function); 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); diff --git a/src/nmod_poly/geometric_progression.c b/src/nmod_poly/geometric_progression.c index 82ce6934cd..139fcde853 100644 --- a/src/nmod_poly/geometric_progression.c +++ b/src/nmod_poly/geometric_progression.c @@ -13,9 +13,9 @@ #include "nmod.h" #include "nmod_vec.h" #include "nmod_poly.h" -#include "impl.h" /* 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, @@ -23,7 +23,7 @@ void _nmod_geometric_progression_evaluate_init_nonfullword(nmod_geometric_progre ulong inv_q, ulong inv_q_pr_quo, ulong inv_q_pr_rem) { /* G->ev_f = sum_{0 <= i < 2*len - 1} q**(i*i/2) * x**i */ - /* G->ev_s[i] = 1 / q**(i*i/2) */ + /* 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); @@ -50,6 +50,7 @@ void _nmod_geometric_progression_evaluate_init_nonfullword(nmod_geometric_progre } /* 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) @@ -75,6 +76,7 @@ void _nmod_geometric_progression_evaluate_init(nmod_geometric_progression_t G, } } +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, @@ -130,6 +132,7 @@ void _nmod_geometric_progression_interpolate_init_nonfullword(nmod_geometric_pro } } +static void _nmod_geometric_progression_interpolate_init(nmod_geometric_progression_t G, slong len, nmod_t mod, ulong q, ulong inv_q) @@ -227,6 +230,7 @@ void _nmod_geometric_progression_init_function(nmod_geometric_progression_t G, } /* 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 */ diff --git a/src/nmod_poly/impl.h b/src/nmod_poly/impl.h index 97dc5c0a9f..80be7214a7 100644 --- a/src/nmod_poly/impl.h +++ b/src/nmod_poly/impl.h @@ -13,7 +13,6 @@ #define NMOD_POLY_IMPL_H #include "nmod_types.h" -#include "nmod_poly.h" /* TODO is this ok? */ void _nmod_poly_inv_series_basecase_preinv1(nn_ptr Qinv, nn_srcptr Q, slong Qlen, slong n, ulong q, nmod_t mod); void _nmod_poly_div_series_basecase_preinv1(nn_ptr Qinv, nn_srcptr P, slong Plen, nn_srcptr Q, slong Qlen, slong n, ulong q, nmod_t mod); @@ -23,20 +22,4 @@ int nmod_poly_irreducible_tetranomial(nmod_poly_t res, ulong n); int nmod_poly_irreducible_pentanomial(nmod_poly_t res, ulong n); void _nmod_poly_divrem_q0_preinv1(nn_ptr Q, nn_ptr R, nn_srcptr A, nn_srcptr B, slong lenA, ulong invL, nmod_t mod); -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); -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); -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); -void _nmod_geometric_progression_interpolate_init(nmod_geometric_progression_t G, slong len, nmod_t mod, - ulong q, ulong inv_q); - - #endif diff --git a/src/nmod_poly/interpolate_geometric_nmod_vec.c b/src/nmod_poly/interpolate_geometric_nmod_vec.c index f4e791bfc6..9f80619894 100644 --- a/src/nmod_poly/interpolate_geometric_nmod_vec.c +++ b/src/nmod_poly/interpolate_geometric_nmod_vec.c @@ -26,10 +26,6 @@ void _nmod_poly_interpolate_geometric_nmod_vec_fast_precomp(nn_ptr poly, return; } - slong f_len, h_len, val; - nn_ptr f = _nmod_vec_init(len); - nn_ptr h = _nmod_vec_init(len); - /** step1: Newton interpolation * [Bostan - Schost, J.Complexity 2005, Section 5.1] * -> The coefficients of the interpolant, in the Newton basis associated @@ -46,7 +42,8 @@ void _nmod_poly_interpolate_geometric_nmod_vec_fast_precomp(nn_ptr poly, */ /* val = valuation of output poly in Newton basis */ - for (val = 0; val < len; val++) + slong val = 0; + for (; val < len; val++) if (v[val] != 0) break; if (val == len) @@ -55,6 +52,10 @@ void _nmod_poly_interpolate_geometric_nmod_vec_fast_precomp(nn_ptr poly, return; } + slong f_len, h_len; + nn_ptr f = _nmod_vec_init(len); + nn_ptr h = _nmod_vec_init(len); + /* actual length of f */ for (f_len = len; f_len > val; f_len--) if (v[f_len-1] != 0) @@ -139,10 +140,10 @@ void nmod_poly_interpolate_geometric_nmod_vec_fast_precomp(nmod_poly_t poly, } void nmod_poly_interpolate_geometric_nmod_vec_fast(nmod_poly_t poly, - ulong r, nn_srcptr ys, slong n) + ulong r, nn_srcptr ys, slong len) { nmod_geometric_progression_t G; - _nmod_geometric_progression_init_function(G, r, n, poly->mod, UWORD(2)); - nmod_poly_interpolate_geometric_nmod_vec_fast_precomp(poly, ys, G, n); - _nmod_geometric_progression_clear_function(G, UWORD(2)); + _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 index e2c388dd70..d0ec69839b 100644 --- a/src/nmod_poly/profile/p-eval_interp_precomputations.c +++ b/src/nmod_poly/profile/p-eval_interp_precomputations.c @@ -156,7 +156,7 @@ void sample_geometric_precomp_eval(void * arg, ulong count) for (ulong ii = 0; ii < __NB_ITER; ii++) { _nmod_geometric_progression_init_function(G, r, npoints_precomp, mod, UWORD(1)); - _nmod_geometric_progression_clear_function(G, UWORD(1)); + nmod_geometric_progression_clear(G); } prof_stop(); } @@ -193,7 +193,7 @@ void sample_geometric_precomp_interp(void * arg, ulong count) for (ulong ii = 0; ii < __NB_ITER; ii++) { _nmod_geometric_progression_init_function(G, r, npoints_precomp, mod, UWORD(2)); - _nmod_geometric_progression_clear_function(G, UWORD(2)); + nmod_geometric_progression_clear(G); } prof_stop(); } From 6df0657d9d3f750c1b2a64a117c9d13766f44197 Mon Sep 17 00:00:00 2001 From: Vincent Neiger Date: Mon, 4 May 2026 16:56:34 +0200 Subject: [PATCH 13/13] fix long -> ulong in profile code ; nmod_find_root --- src/nmod_poly/profile/p-eval_interp_precomputations.c | 2 +- src/nmod_poly/profile/p-evaluate_nmod_vec.c | 2 +- src/nmod_poly/profile/p-interpolate_nmod_vec.c | 2 +- 3 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/nmod_poly/profile/p-eval_interp_precomputations.c b/src/nmod_poly/profile/p-eval_interp_precomputations.c index d0ec69839b..d5571f3479 100644 --- a/src/nmod_poly/profile/p-eval_interp_precomputations.c +++ b/src/nmod_poly/profile/p-eval_interp_precomputations.c @@ -24,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++) diff --git a/src/nmod_poly/profile/p-evaluate_nmod_vec.c b/src/nmod_poly/profile/p-evaluate_nmod_vec.c index ae63472f5b..9fbeb63804 100644 --- a/src/nmod_poly/profile/p-evaluate_nmod_vec.c +++ b/src/nmod_poly/profile/p-evaluate_nmod_vec.c @@ -24,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++) diff --git a/src/nmod_poly/profile/p-interpolate_nmod_vec.c b/src/nmod_poly/profile/p-interpolate_nmod_vec.c index ad806dc2a3..7ff726b26b 100644 --- a/src/nmod_poly/profile/p-interpolate_nmod_vec.c +++ b/src/nmod_poly/profile/p-interpolate_nmod_vec.c @@ -24,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++)