From d56dd96e992e8d4c7b57ba2299976130b1c904a3 Mon Sep 17 00:00:00 2001 From: Vincent Neiger Date: Sat, 9 May 2026 10:22:25 +0200 Subject: [PATCH 1/6] - add function for polynomial extrapolation at points in geometric sequence - add test file for geometric extrapolation Co-authored-by: ktran11 --- src/nmod_poly.h | 10 +- .../extrapolate_geometric_nmod_vec.c | 128 ++++++++++++++ src/nmod_poly/geometric_progression.c | 72 +++++++- src/nmod_poly/test/main.c | 2 + src/nmod_poly/test/t-extrapolate_geometric.c | 157 ++++++++++++++++++ 5 files changed, 362 insertions(+), 7 deletions(-) create mode 100644 src/nmod_poly/extrapolate_geometric_nmod_vec.c create mode 100644 src/nmod_poly/test/t-extrapolate_geometric.c diff --git a/src/nmod_poly.h b/src/nmod_poly.h index 251eab32be..3e24c1a317 100644 --- a/src/nmod_poly.h +++ b/src/nmod_poly.h @@ -64,6 +64,11 @@ typedef struct nmod_poly_t ev_f; /* evaluate: polynomial */ nn_ptr int_s1, int_s2; /* interpolate: scaling constants */ nmod_poly_t int_f; /* interpolate: polynomial */ + nmod_poly_t ext_ff; /* extrapolate forward: polynomial */ + nmod_poly_t ext_fb; /* extrapolate backward: polynomial */ + nn_ptr ext_s1f; /* extrapolate: forward scaling constant */ + nn_ptr ext_s1b; /* extrapolate: backward scaling constant */ + nn_ptr ext_s2, ext_s3; /* extrapolate: shared scaling constants */ nmod_t mod; /* modulus */ slong len; /* number of points */ ulong function; /* choice of precomputations */ @@ -592,7 +597,7 @@ void _nmod_poly_tree_free(nn_ptr * tree, slong len); void _nmod_poly_tree_build(nn_ptr * tree, nn_srcptr roots, slong len, nmod_t mod); -/* Geometric evaluation / interpolation *************************************/ +/* Geometric evaluation / interpolation / extrapolation *********************/ void _nmod_geometric_progression_init_function(nmod_geometric_progression_t G, ulong r, slong len, nmod_t mod, @@ -613,6 +618,9 @@ void _nmod_poly_interpolate_geometric_nmod_vec_fast_precomp(nn_ptr poly, nn_srcp 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 len); +void nmod_poly_extrapolate_geometric(nn_ptr oval, slong olen, nn_srcptr ival, slong ilen, slong offset, ulong r, nmod_t mod); +void nmod_poly_extrapolate_geometric_precomp(nn_ptr oval, slong olen, nn_srcptr ival, slong ilen, slong offset, const nmod_geometric_progression_t G); + /* Interpolation ************************************************************/ void _nmod_poly_interpolate_nmod_vec_newton(nn_ptr poly, nn_srcptr xs, nn_srcptr ys, slong len, nmod_t mod); diff --git a/src/nmod_poly/extrapolate_geometric_nmod_vec.c b/src/nmod_poly/extrapolate_geometric_nmod_vec.c new file mode 100644 index 0000000000..b11a244084 --- /dev/null +++ b/src/nmod_poly/extrapolate_geometric_nmod_vec.c @@ -0,0 +1,128 @@ +/* + Copyright (C) 2026 Vincent Neiger, Kevin Tran + + 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 "nmod.h" +#include "nmod_poly.h" + +void nmod_poly_extrapolate_geometric_precomp(nn_ptr oval, slong olen, + nn_srcptr ival, slong ilen, + slong offset, + const nmod_geometric_progression_t G) +{ + /* precomputation has been done */ + FLINT_ASSERT((G->function >> 2) & 1); + /* input/output points are disjoint, and stay within precomputed data length */ + FLINT_ASSERT((offset >= ilen && G->len >= offset+olen) + || (offset <= -olen && G->len >= ilen-offset)); + + if (ilen == 0) + { + for (slong i = 0; i < olen; i++) + oval[i] = 0; + return; + } + + if (ilen == 1) + { + for (slong i = 0; i < olen; i++) + oval[i] = ival[0]; + return; + } + + /* forward extrapolation */ + if (offset > 0) + { + /* TODO TMP ALLOC? */ + nn_ptr tmp = FLINT_ARRAY_ALLOC(ilen, ulong); + /* first scaling */ + // svals = [ext_s2[i] * ext_s3[m-1-i] * vals[i] for i in range(m)] + for (slong i = 0; i < ilen; i++) + { + /* TODO use some NMOD_RED3?? */ + tmp[i] = nmod_mul(G->ext_s3[ilen-1-i], ival[i], G->mod); + tmp[i] = nmod_mul(G->ext_s2[i], tmp[i], G->mod); + } + + /* middle product */ + // g = xring(svals) + // f = ext_ff.shift(-(lmk - m)).truncate(m+n-1) + // mp = (f * g).shift(-(m-1)).truncate(n) + _nmod_poly_mulmid(oval, G->ext_ff->coeffs + offset-ilen, ilen+olen-1, tmp, ilen, ilen-1, ilen+olen-1, G->mod); + + /* second scaling */ + // w = [ext_s1f[lmk+j] * ext_s2[lmk-m+j] * mp[j] for j in range(n)] + for (slong j = 0; j < olen; j++) + { + oval[j] = nmod_mul(G->ext_s2[offset-ilen+j], oval[j], G->mod); + oval[j] = nmod_mul(G->ext_s1f[offset+j], oval[j], G->mod); + } + flint_free(tmp); + } + + /* backward extrapolation */ + else + { + /* TODO TMP ALLOC? */ + nn_ptr tmp = FLINT_ARRAY_ALLOC(FLINT_MAX(ilen, olen), ulong); + /* first scaling */ + // svals = [ext_s2[m-1-i] * ext_s3[i] * vals[m-1-i] for i in range(m)] + for (slong i = 0; i < ilen; i++) + { + /* TODO use some NMOD_RED3?? */ + tmp[i] = nmod_mul(G->ext_s2[ilen-1-i], ival[ilen-1-i], G->mod); + tmp[i] = nmod_mul(G->ext_s3[i], tmp[i], G->mod); + } + + /* middle product */ + // g = xring(svals) + // f = ext_fb.shift(lmk + n).truncate(m+n-1) + // mp = (f * g).shift(-(m-1)).truncate(n) + _nmod_poly_mulmid(oval, G->ext_fb->coeffs - (offset+olen), ilen+olen-1, tmp, ilen, ilen-1, ilen+olen-1, G->mod); + + /* second scaling */ + // w = [ext_s1b[m-1-lmk-j] * ext_s3[-lmk-j-1] * mp[n-1-j] for j in range(n)] + for (slong j = 0; j < olen; j++) + tmp[j] = oval[olen - 1 - j]; + for (slong j = 0; j < olen; j++) + { + oval[j] = nmod_mul(G->ext_s3[-offset-1-j], tmp[j], G->mod); + oval[j] = nmod_mul(G->ext_s1b[ilen-1-offset-j], oval[j], G->mod); + } + + flint_free(tmp); + } +} + +void nmod_poly_extrapolate_geometric(nn_ptr oval, slong olen, + nn_srcptr ival, slong ilen, + slong offset, ulong r, nmod_t mod) +{ + /* TODO handle overlap */ + if (ilen == 0) + { + for (slong i = 0; i < olen; i++) + oval[i] = 0; + return; + } + + if (ilen == 1) + { + for (slong i = 0; i < olen; i++) + oval[i] = ival[0]; + return; + } + + nmod_geometric_progression_t G; + slong len = (offset > 0) ? offset+olen : ilen-offset; + _nmod_geometric_progression_init_function(G, r, len, mod, UWORD(4)); + nmod_poly_extrapolate_geometric_precomp(oval, olen, ival, ilen, offset, G); + nmod_geometric_progression_clear(G); +} diff --git a/src/nmod_poly/geometric_progression.c b/src/nmod_poly/geometric_progression.c index 139fcde853..6a2b9fb952 100644 --- a/src/nmod_poly/geometric_progression.c +++ b/src/nmod_poly/geometric_progression.c @@ -182,6 +182,55 @@ void _nmod_geometric_progression_interpolate_init(nmod_geometric_progression_t G } } +static +void _nmod_geometric_progression_extrapolate_init(nmod_geometric_progression_t G, + slong len, nmod_t mod, + ulong q, ulong inv_q) +{ + /* ext_ff: sum_{i=0}^{len-1} x**i / (q**(i+1) - 1) [for forward] */ + /* ext_fb: sum_{i=0}^{len-1} x**i / (q**(-i-1) - 1) [for backward] */ + nmod_poly_init2_preinv(G->ext_ff, mod.n, mod.ninv, len-1); + nmod_poly_init2_preinv(G->ext_fb, mod.n, mod.ninv, len-1); + G->ext_ff->length = len-1; + G->ext_fb->length = len-1; + + /* ext_s1f[i] = prod_{k=1}^{i} (q**k - 1) [for forward] */ /* TODO s0f */ + /* ext_s1b[i] = prod_{k=1}^{i} (q**(-k) - 1) [for backward] */ /* TODO s0b */ + /* ext_s2[i] = 1 / ext_s1f[i] [for both] */ /* TODO s1p */ + /* ext_s3[i] = 1 / ext_s1b[i] [for both] */ /* TODO s1m */ + G->ext_s1f = _nmod_vec_init(len); + G->ext_s1b = _nmod_vec_init(len); + G->ext_s2 = _nmod_vec_init(len); + G->ext_s3 = _nmod_vec_init(len); + + G->ext_s1f[0] = 1; + G->ext_s1b[0] = 1; + + ulong q_pow_i = q; + ulong inv_q_pow_i = inv_q; + + for (slong i = 1; i < len; i++) + { + G->ext_ff->coeffs[i-1] = q_pow_i - 1; + G->ext_fb->coeffs[i-1] = inv_q_pow_i - 1; + G->ext_s1f[i] = nmod_mul(G->ext_s1f[i-1], q_pow_i - 1, mod); + G->ext_s1b[i] = nmod_mul(G->ext_s1b[i-1], inv_q_pow_i - 1, mod); + q_pow_i = nmod_mul(q_pow_i, q, mod); + inv_q_pow_i = nmod_mul(inv_q_pow_i, inv_q, mod); + } + + G->ext_s2[len-1] = nmod_inv(G->ext_s1f[len-1], mod); + G->ext_s3[len-1] = nmod_inv(G->ext_s1b[len-1], mod); + + for (slong i = len-1; i > 0; i--) + { + G->ext_s2[i-1] = nmod_mul(G->ext_s2[i], G->ext_ff->coeffs[i-1], mod); + G->ext_s3[i-1] = nmod_mul(G->ext_s3[i], G->ext_fb->coeffs[i-1], mod); + G->ext_ff->coeffs[i-1] = nmod_mul(G->ext_s2[i], G->ext_s1f[i-1], mod); + G->ext_fb->coeffs[i-1] = nmod_mul(G->ext_s3[i], G->ext_s1b[i-1], mod); + } +} + /* initialize for selection of functionalities: */ /* the lowest 3 bits of `function` act as a mask for the */ /* three functionalities, in this order: */ @@ -211,8 +260,11 @@ void _nmod_geometric_progression_init_function(nmod_geometric_progression_t G, _nmod_geometric_progression_evaluate_init_nonfullword(G, r, len, mod, q, q_pr_quo, q_pr_rem, inv_r, inv_q, inv_q_pr_quo, inv_q_pr_rem); if ((function>>1) & UWORD(1)) /* interpolate */ _nmod_geometric_progression_interpolate_init_nonfullword(G, len, mod, q, q_pr_quo, q_pr_rem, inv_q, inv_q_pr_quo, inv_q_pr_rem); - /* if ((function>>2) & UWORD(1)) */ - /* extrapolate */ + /* TODO improve Shoup interpolate */ + if ((function>>2) & UWORD(1)) /* extrapolate */ + _nmod_geometric_progression_extrapolate_init(G, len, mod, q, inv_q); + /* TODO Shoup extrapolate */ + /* TODO Shoup version of combine extrapolate&interpolate */ } else { @@ -224,8 +276,9 @@ void _nmod_geometric_progression_init_function(nmod_geometric_progression_t G, _nmod_geometric_progression_evaluate_init(G, r, len, mod, q, inv_r, inv_q); if ((function>>1) & UWORD(1)) /* interpolate */ _nmod_geometric_progression_interpolate_init(G, len, mod, q, inv_q); - /* if ((function>>2) & UWORD(1)) */ - /* extrapolate */ + if ((function>>2) & UWORD(1)) /* extrapolate */ + _nmod_geometric_progression_extrapolate_init(G, len, mod, q, inv_q); + /* TODO combine extrapolate&interpolate (note: ext_s1f is the same as G->int_s1) */ } } @@ -244,8 +297,15 @@ void _nmod_geometric_progression_clear_function(nmod_geometric_progression_t G, _nmod_vec_clear(G->int_s2); nmod_poly_clear(G->int_f); } - /* if ((function>>2) & UWORD(1)) */ - /* extrapolate */ + if ((function>>2) & UWORD(1)) /* extrapolate */ + { + _nmod_vec_clear(G->ext_s1f); + _nmod_vec_clear(G->ext_s1b); + _nmod_vec_clear(G->ext_s2); + _nmod_vec_clear(G->ext_s3); + nmod_poly_clear(G->ext_ff); + nmod_poly_clear(G->ext_fb); + } } /* initialize for all: evaluate+interpolate+extrapolate */ diff --git a/src/nmod_poly/test/main.c b/src/nmod_poly/test/main.c index bd8f9e874e..8ad8dcf9de 100644 --- a/src/nmod_poly/test/main.c +++ b/src/nmod_poly/test/main.c @@ -54,6 +54,7 @@ #include "t-evaluate_nmod_vec_fast.c" #include "t-evaluate_geometric_nmod_vec_fast.c" #include "t-exp_series.c" +#include "t-extrapolate_geometric.c" #include "t-find_distinct_nonzero_roots.c" #include "t-fread_print.c" #include "t-gcd.c" @@ -187,6 +188,7 @@ test_struct tests[] = TEST_FUNCTION(nmod_poly_evaluate_nmod_vec_fast), TEST_FUNCTION(nmod_poly_evaluate_geometric_nmod_vec_fast), TEST_FUNCTION(nmod_poly_exp_series), + TEST_FUNCTION(nmod_poly_extrapolate_geometric), TEST_FUNCTION(nmod_poly_find_distinct_nonzero_roots), TEST_FUNCTION(nmod_poly_fread_print), TEST_FUNCTION(nmod_poly_gcd), diff --git a/src/nmod_poly/test/t-extrapolate_geometric.c b/src/nmod_poly/test/t-extrapolate_geometric.c new file mode 100644 index 0000000000..37e59e86a2 --- /dev/null +++ b/src/nmod_poly/test/t-extrapolate_geometric.c @@ -0,0 +1,157 @@ +/* + Copyright (C) 2026 Vincent Neiger, Kevin Tran + + 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 "nmod.h" +#include "test_helpers.h" +#include "nmod_vec.h" +#include "nmod_poly.h" + +TEST_FUNCTION_START(nmod_poly_extrapolate_geometric, state) +{ + int i, result = 1; + + for (i = 0; i < 200 * flint_test_multiplier(); i++) + { + nmod_poly_t P; + nn_ptr ival, oval, val; + nn_ptr ipts, opts; + ulong modn, r, q; + nmod_t mod; + slong ilen, olen, istart, offset; + + ilen = (i < 10) ? i : n_randint(state, 500); + olen = (i < 10) ? i : n_randint(state, 500); + + ival = FLINT_ARRAY_ALLOC(ilen, ulong); + oval = FLINT_ARRAY_ALLOC(olen, ulong); + val = FLINT_ARRAY_ALLOC(olen, ulong); + ipts = FLINT_ARRAY_ALLOC(ilen, ulong); + opts = FLINT_ARRAY_ALLOC(olen, ulong); + + /* forward: offset >= ilen */ + { + istart = n_randint(state, 50); + offset = ilen + n_randint(state, 10); + + /* choose modn large enough for r to exist: */ + do + { + modn = n_randtest_prime(state, 1); + } + while (modn <= 2*((ulong)istart+offset+olen) + 1); + + nmod_init(&mod, modn); + nmod_poly_init(P, modn); + r = n_primitive_root_prime(modn); + q = nmod_mul(r, r, mod); + + _nmod_vec_randtest(ival, state, ilen, mod); + + /* input points are c * q**k, k = istart...istart+ilen-1 */ + ulong c = 1 + n_randint(state, modn-1); + ipts[0] = nmod_pow_ui(q, istart, mod); + ipts[0] = nmod_mul(ipts[0], c, mod); + for (slong k = 1; k < ilen; k++) + ipts[k] = nmod_mul(ipts[k-1], q, mod); + + /* output points are c * q**k, k = istart+offset...istart+offset+olen-1 */ + opts[0] = nmod_pow_ui(q, istart+offset, mod); + opts[0] = nmod_mul(opts[0], c, mod); + for (slong k = 1; k < olen; k++) + opts[k] = nmod_mul(opts[k-1], q, mod); + + nmod_poly_interpolate_nmod_vec(P, ipts, ival, ilen); + nmod_poly_evaluate_nmod_vec(val, P, opts, olen); + + nmod_poly_extrapolate_geometric(oval, olen, ival, ilen, offset, r, mod); + + result = _nmod_vec_equal(oval, val, olen); + + if (!result) + { + flint_printf("FAIL (forward):\n"); + flint_printf("mod=%wu, ilen=%wd, olen=%wd, offset=%wd\n\n", modn, ilen, olen, offset); + flint_printf("ival: "); _nmod_vec_print_pretty(ival, ilen, mod); flint_printf("\n\n"); + flint_printf("oval: "); _nmod_vec_print_pretty(oval, olen, mod); flint_printf("\n\n"); + flint_printf("correct: "); _nmod_vec_print_pretty(val, olen, mod); flint_printf("\n\n"); + fflush(stdout); + flint_abort(); + } + + nmod_poly_clear(P); + } + + /* backward: offset + olen <= 0 */ + { + /* here we build the points explicitly, so we also need istart+offset >= 0 */ + istart = olen + 10 + n_randint(state, 50); + offset = - (olen + (slong)n_randint(state, 10)); + + /* choose modn large enough for r to exist */ + do + { + modn = n_randtest_prime(state, 1); + } + while (modn <= 2*((ulong)istart+ilen) + 1); + + nmod_init(&mod, modn); + nmod_poly_init(P, modn); + r = n_primitive_root_prime(modn); + q = nmod_mul(r, r, mod); + + _nmod_vec_randtest(ival, state, ilen, mod); + + /* input points are c * q**k, k = istart...istart+ilen-1 */ + ulong c = 1 + n_randint(state, modn-1); + ipts[0] = nmod_pow_ui(q, istart, mod); + ipts[0] = nmod_mul(ipts[0], c, mod); + for (slong k = 1; k < ilen; k++) + ipts[k] = nmod_mul(ipts[k-1], q, mod); + + /* output points are c * q**k, k = istart+offset...istart+offset+olen-1 */ + opts[0] = nmod_pow_ui(q, istart+offset, mod); + opts[0] = nmod_mul(opts[0], c, mod); + for (slong k = 1; k < olen; k++) + opts[k] = nmod_mul(opts[k-1], q, mod); + + nmod_poly_interpolate_nmod_vec(P, ipts, ival, ilen); + nmod_poly_evaluate_nmod_vec(val, P, opts, olen); + + nmod_poly_extrapolate_geometric(oval, olen, ival, ilen, offset, r, mod); + + result = _nmod_vec_equal(oval, val, olen); + + if (!result) + { + flint_printf("FAIL (backward):\n"); + flint_printf("mod=%wu, ilen=%wd, olen=%wd, ", modn, ilen, olen); + flint_printf("istart=%wd, offset=%wd, q=%wu, c=%wu\n\n", istart, offset, q, c); + if (ilen < 20 && olen < 20) + { + flint_printf("ival: "); _nmod_vec_print_pretty(ival, ilen, mod); flint_printf("\n\n"); + flint_printf("oval: "); _nmod_vec_print_pretty(oval, olen, mod); flint_printf("\n\n"); + flint_printf("correct: "); _nmod_vec_print_pretty(val, olen, mod); flint_printf("\n\n"); + } + fflush(stdout); + flint_abort(); + } + nmod_poly_clear(P); + } + + flint_free(ival); + flint_free(oval); + flint_free(val); + flint_free(ipts); + flint_free(opts); + } + + TEST_FUNCTION_END(state); +} From 410a54146b01d2a7c01e2c7964d15b6241fda5d1 Mon Sep 17 00:00:00 2001 From: Vincent Neiger Date: Sat, 9 May 2026 15:31:14 +0200 Subject: [PATCH 2/6] Add profile for geometric extrapolate + minor enhancements Co-authored-by: ktran11 --- .../extrapolate_geometric_nmod_vec.c | 29 +- .../profile/p-eval_interp_precomputations.c | 55 +- .../profile/p-extrapolate_geometric.c | 683 ++++++++++++++++++ 3 files changed, 742 insertions(+), 25 deletions(-) create mode 100644 src/nmod_poly/profile/p-extrapolate_geometric.c diff --git a/src/nmod_poly/extrapolate_geometric_nmod_vec.c b/src/nmod_poly/extrapolate_geometric_nmod_vec.c index b11a244084..0a21f4cfa8 100644 --- a/src/nmod_poly/extrapolate_geometric_nmod_vec.c +++ b/src/nmod_poly/extrapolate_geometric_nmod_vec.c @@ -40,55 +40,48 @@ void nmod_poly_extrapolate_geometric_precomp(nn_ptr oval, slong olen, /* forward extrapolation */ if (offset > 0) { - /* TODO TMP ALLOC? */ - nn_ptr tmp = FLINT_ARRAY_ALLOC(ilen, ulong); + TMP_INIT; + TMP_START; + nn_ptr tmp = TMP_ALLOC(ilen * sizeof(ulong)); + /* first scaling */ - // svals = [ext_s2[i] * ext_s3[m-1-i] * vals[i] for i in range(m)] for (slong i = 0; i < ilen; i++) { - /* TODO use some NMOD_RED3?? */ tmp[i] = nmod_mul(G->ext_s3[ilen-1-i], ival[i], G->mod); tmp[i] = nmod_mul(G->ext_s2[i], tmp[i], G->mod); } /* middle product */ - // g = xring(svals) - // f = ext_ff.shift(-(lmk - m)).truncate(m+n-1) - // mp = (f * g).shift(-(m-1)).truncate(n) _nmod_poly_mulmid(oval, G->ext_ff->coeffs + offset-ilen, ilen+olen-1, tmp, ilen, ilen-1, ilen+olen-1, G->mod); /* second scaling */ - // w = [ext_s1f[lmk+j] * ext_s2[lmk-m+j] * mp[j] for j in range(n)] for (slong j = 0; j < olen; j++) { oval[j] = nmod_mul(G->ext_s2[offset-ilen+j], oval[j], G->mod); oval[j] = nmod_mul(G->ext_s1f[offset+j], oval[j], G->mod); } - flint_free(tmp); + + TMP_END; } /* backward extrapolation */ else { - /* TODO TMP ALLOC? */ - nn_ptr tmp = FLINT_ARRAY_ALLOC(FLINT_MAX(ilen, olen), ulong); + TMP_INIT; + TMP_START; + nn_ptr tmp = TMP_ALLOC(FLINT_MAX(ilen, olen) * sizeof(ulong)); + /* first scaling */ - // svals = [ext_s2[m-1-i] * ext_s3[i] * vals[m-1-i] for i in range(m)] for (slong i = 0; i < ilen; i++) { - /* TODO use some NMOD_RED3?? */ tmp[i] = nmod_mul(G->ext_s2[ilen-1-i], ival[ilen-1-i], G->mod); tmp[i] = nmod_mul(G->ext_s3[i], tmp[i], G->mod); } /* middle product */ - // g = xring(svals) - // f = ext_fb.shift(lmk + n).truncate(m+n-1) - // mp = (f * g).shift(-(m-1)).truncate(n) _nmod_poly_mulmid(oval, G->ext_fb->coeffs - (offset+olen), ilen+olen-1, tmp, ilen, ilen-1, ilen+olen-1, G->mod); /* second scaling */ - // w = [ext_s1b[m-1-lmk-j] * ext_s3[-lmk-j-1] * mp[n-1-j] for j in range(n)] for (slong j = 0; j < olen; j++) tmp[j] = oval[olen - 1 - j]; for (slong j = 0; j < olen; j++) @@ -97,7 +90,7 @@ void nmod_poly_extrapolate_geometric_precomp(nn_ptr oval, slong olen, oval[j] = nmod_mul(G->ext_s1b[ilen-1-offset-j], oval[j], G->mod); } - flint_free(tmp); + TMP_END; } } diff --git a/src/nmod_poly/profile/p-eval_interp_precomputations.c b/src/nmod_poly/profile/p-eval_interp_precomputations.c index d5571f3479..d278d1723b 100644 --- a/src/nmod_poly/profile/p-eval_interp_precomputations.c +++ b/src/nmod_poly/profile/p-eval_interp_precomputations.c @@ -91,7 +91,7 @@ void sample_general_precomp(void * arg, ulong count) FLINT_TEST_CLEAR(state); } -/* precomputations for geometric progression (eval+interp) */ +/* precomputations for geometric progression (eval+interp+extrap) */ void sample_geometric_precomp(void * arg, ulong count) { ulong n; @@ -201,6 +201,43 @@ void sample_geometric_precomp_interp(void * arg, ulong count) FLINT_TEST_CLEAR(state); } +/* precomputations for geometric progression (extrap only) */ +void sample_geometric_precomp_extrap(void * arg, ulong count) +{ + ulong n; + nmod_t mod; + ulong i; + + info_t * info = (info_t *) arg; + flint_bitcnt_t bits = info->bits; + slong npoints_precomp = info->npoints_precomp; + + FLINT_TEST_INIT(state); + + for (i = 0; i < count; i++) + { + n = n_randprime(state, bits, 1); + nmod_init(&mod, n); + ulong r = nmod_find_root(2*npoints_precomp, mod); + if (r == 0) + flint_printf("\n...could not find element of suitable order for geometric progression...\n"); + + nmod_geometric_progression_t G; + G->len = npoints_precomp; + G->mod = mod; + + prof_start(); + for (ulong ii = 0; ii < __NB_ITER; ii++) + { + _nmod_geometric_progression_init_function(G, r, npoints_precomp, mod, UWORD(4)); + nmod_geometric_progression_clear(G); + } + prof_stop(); + } + + FLINT_TEST_CLEAR(state); +} + int main(int argc, char * argv[]) { if (argc > 1 && (strcmp(argv[1], "-h") == 0 || strcmp(argv[1], "--help") == 0)) @@ -233,6 +270,7 @@ int main(int argc, char * argv[]) double time_geometric_all; double time_geometric_eval; double time_geometric_interp; + double time_geometric_extrap; info_t info; flint_bitcnt_t i; @@ -246,11 +284,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 | g-interp \n"); + flint_printf("len\tpoints | general | geom | g-eval | g-interp| g-extrap \n"); else if (func_bench == 1) /* general only */ flint_printf("len\tpoints | general\n"); else if (func_bench == 2) /* geometric only */ - flint_printf("len\tpoints | geom | g-eval | g-interp \n"); + flint_printf("len\tpoints | geom | g-eval | g-interp| g-extrap \n"); for (int len = 0; len < nb_lens; ++len) { @@ -267,16 +305,18 @@ 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); + prof_repeat(&time_geometric_extrap, &tmp, sample_geometric_precomp_extrap, (void *) &info); } if (func_bench == 0) { - flint_printf("%ld\t%7ld| %.1e | %.1e | %.1e | %.1e\n", + flint_printf("%ld\t%7ld| %.1e | %.1e | %.1e | %.1e | %.1e\n", info.length, info.npoints_precomp, time_general/fac, time_geometric_all/fac, time_geometric_eval/fac, - time_geometric_interp/fac); + time_geometric_interp/fac, + time_geometric_extrap/fac); } else if (func_bench == 1) { @@ -286,11 +326,12 @@ int main(int argc, char * argv[]) } else if (func_bench == 2) { - 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_geometric_all/fac, time_geometric_eval/fac, - time_geometric_interp/fac); + time_geometric_interp/fac, + time_geometric_extrap/fac); } } diff --git a/src/nmod_poly/profile/p-extrapolate_geometric.c b/src/nmod_poly/profile/p-extrapolate_geometric.c new file mode 100644 index 0000000000..47c7a73370 --- /dev/null +++ b/src/nmod_poly/profile/p-extrapolate_geometric.c @@ -0,0 +1,683 @@ +/* + Copyright (C) 2026 Vincent Neiger, Kevin Tran + + This file is part of FLINT. + + FLINT is free software: you can redistribute it and/or modify it under + the terms of the GNU Lesser General Public License (LGPL) as published + by the Free Software Foundation; either version 3 of the License, or + (at your option) any later version. See . +*/ + +#include /* for atoi */ +#include + +#include "nmod.h" +#include "profiler.h" +#include "ulong_extras.h" +#include "nmod_poly.h" +#include "nmod_vec.h" + +#define __NB_ITER 10 + +/*------------------------------------------------------------*/ +/* finds an element of order at least n */ +/* returns 0 if not found */ +/*------------------------------------------------------------*/ +static ulong nmod_find_root(slong n, nmod_t mod) +{ + ulong attempts = 0; + for (ulong q = 2; q < mod.n; q++) + { + slong k = 1; + slong qk = q; + while (qk != 1 && k < n) + { + qk = nmod_mul(qk, q, mod); + k++; + } + if (qk != 1) + { + return q; + } + attempts += 1; + if (attempts >= 10) + return 0; + } + return 0; +} + +typedef struct +{ + flint_bitcnt_t bits; + slong ilen; + slong olen; +} info_t; + +/* general points */ +void sample_general(void * arg, ulong count) +{ + ulong n; + nmod_t mod; + ulong i; + + info_t * info = (info_t *) arg; + flint_bitcnt_t bits = info->bits; + slong ilen = info->ilen; + slong olen = info->olen; + + FLINT_TEST_INIT(state); + + nn_ptr ipts = FLINT_ARRAY_ALLOC(ilen, ulong); + nn_ptr opts = FLINT_ARRAY_ALLOC(olen, ulong); + nn_ptr ival = FLINT_ARRAY_ALLOC(ilen, ulong); + nn_ptr oval = FLINT_ARRAY_ALLOC(olen, ulong); + + for (i = 0; i < count; i++) + { + n = n_randprime(state, bits, 0); + if (n == UWORD(0)) n++; + nmod_init(&mod, n); + + /* need pairwise distinct ipts for interpolation, let's take ipts[k] = k */ + for (slong k = 0; k < ilen; k++) + ipts[k] = k; + _nmod_vec_rand(opts, state, olen, mod); + _nmod_vec_rand(ival, state, ilen, mod); + + prof_start(); + nn_ptr poly = FLINT_ARRAY_ALLOC(ilen, ulong); + for (ulong ii = 0; ii < __NB_ITER; ii++) + { + _nmod_poly_interpolate_nmod_vec(poly, ipts, ival, ilen, mod); + _nmod_poly_evaluate_nmod_vec(oval, poly, ilen, opts, olen, mod); + } + flint_free(poly); + prof_stop(); + } + + flint_free(ipts); + flint_free(opts); + flint_free(ival); + flint_free(oval); + + FLINT_TEST_CLEAR(state); +} + +/* general points, not counting tree precomputations */ +void sample_general_noprecomp(void * arg, ulong count) +{ + ulong n; + nmod_t mod; + ulong i; + + info_t * info = (info_t *) arg; + flint_bitcnt_t bits = info->bits; + slong ilen = info->ilen; + slong olen = info->olen; + + FLINT_TEST_INIT(state); + + nn_ptr ipts = FLINT_ARRAY_ALLOC(ilen, ulong); + nn_ptr opts = FLINT_ARRAY_ALLOC(olen, ulong); + nn_ptr ival = FLINT_ARRAY_ALLOC(ilen, ulong); + nn_ptr oval = FLINT_ARRAY_ALLOC(olen, ulong); + + for (i = 0; i < count; i++) + { + n = n_randprime(state, bits, 0); + if (n == UWORD(0)) n++; + nmod_init(&mod, n); + + /* need pairwise distinct ipts for interpolation, let's take ipts[k] = k */ + for (slong k = 0; k < ilen; k++) + ipts[k] = k; + _nmod_vec_rand(opts, state, olen, mod); + _nmod_vec_rand(ival, state, ilen, mod); + + nn_ptr * itree = _nmod_poly_tree_alloc(ilen); + _nmod_poly_tree_build(itree, ipts, ilen, mod); + nn_ptr iw = _nmod_vec_init(ilen); + _nmod_poly_interpolation_weights(iw, itree, ilen, mod); + + nn_ptr * otree = _nmod_poly_tree_alloc(olen); + _nmod_poly_tree_build(otree, opts, olen, mod); + nn_ptr ow = _nmod_vec_init(olen); + _nmod_poly_interpolation_weights(ow, otree, olen, mod); + + prof_start(); + nn_ptr poly = FLINT_ARRAY_ALLOC(ilen, ulong); + for (ulong ii = 0; ii < __NB_ITER; ii++) + { + _nmod_poly_interpolate_nmod_vec_fast_precomp(poly, ival, itree, iw, ilen, mod); + _nmod_poly_evaluate_nmod_vec_fast_precomp(oval, poly, ilen, otree, olen, mod); + } + flint_free(poly); + prof_stop(); + + _nmod_poly_tree_free(itree, ilen); + flint_free(iw); + _nmod_poly_tree_free(otree, olen); + flint_free(ow); + } + + flint_free(ipts); + flint_free(opts); + flint_free(ival); + flint_free(oval); + + FLINT_TEST_CLEAR(state); +} + +/* tree, counting only precomputations */ +void sample_general_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 ilen = info->ilen; + slong olen = info->olen; + + FLINT_TEST_INIT(state); + + nn_ptr ipts = FLINT_ARRAY_ALLOC(ilen, ulong); + nn_ptr opts = FLINT_ARRAY_ALLOC(olen, ulong); + + for (i = 0; i < count; i++) + { + n = n_randprime(state, bits, 0); + if (n == UWORD(0)) n++; + nmod_init(&mod, n); + + /* need pairwise distinct ipts for interpolation, let's take ipts[k] = k */ + for (slong k = 0; k < ilen; k++) + ipts[k] = k; + _nmod_vec_rand(opts, state, olen, mod); + + nn_ptr * itree = _nmod_poly_tree_alloc(ilen); + nn_ptr iw = _nmod_vec_init(ilen); + nn_ptr * otree = _nmod_poly_tree_alloc(olen); + nn_ptr ow = _nmod_vec_init(olen); + + prof_start(); + for (ulong ii = 0; ii < __NB_ITER; ii++) + { + _nmod_poly_tree_build(itree, ipts, ilen, mod); + _nmod_poly_interpolation_weights(iw, itree, ilen, mod); + + _nmod_poly_tree_build(otree, opts, olen, mod); + _nmod_poly_interpolation_weights(ow, otree, olen, mod); + } + prof_stop(); + _nmod_poly_tree_free(itree, ilen); + flint_free(iw); + _nmod_poly_tree_free(otree, olen); + flint_free(ow); + } + + flint_free(ipts); + flint_free(opts); + + FLINT_TEST_CLEAR(state); +} + +/* geometric, including precomp */ +void sample_geometric(void * arg, ulong count) +{ + ulong n; + nmod_t mod; + ulong i; + + info_t * info = (info_t *) arg; + flint_bitcnt_t bits = info->bits; + slong ilen = info->ilen; + slong olen = info->olen; + + FLINT_TEST_INIT(state); + + nn_ptr ival = FLINT_ARRAY_ALLOC(ilen, ulong); + nn_ptr oval = FLINT_ARRAY_ALLOC(olen, ulong); + + for (i = 0; i < count; i++) + { + n = n_randprime(state, bits, 0); + nmod_init(&mod, n); + + _nmod_vec_rand(ival, state, ilen, mod); + ulong r = nmod_find_root(2*(ilen + olen), 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_extrapolate_geometric(oval, olen, ival, ilen, ilen, r, mod); + prof_stop(); + } + + flint_free(ival); + flint_free(oval); + + FLINT_TEST_CLEAR(state); +} + +/* geometric, not counting precomputations */ +void sample_geometric_noprecomp(void * arg, ulong count) +{ + ulong n; + nmod_t mod; + ulong i; + + info_t * info = (info_t *) arg; + flint_bitcnt_t bits = info->bits; + slong ilen = info->ilen; + slong olen = info->olen; + + FLINT_TEST_INIT(state); + + nn_ptr ival = FLINT_ARRAY_ALLOC(ilen, ulong); + nn_ptr oval = FLINT_ARRAY_ALLOC(olen, ulong); + + for (i = 0; i < count; i++) + { + nmod_geometric_progression_t G; + n = n_randprime(state, bits, 0); + nmod_init(&mod, n); + + _nmod_vec_rand(ival, state, ilen, mod); + ulong r = nmod_find_root(2*(ilen + olen), mod); + if (r == 0) + flint_printf("\n...could not find element of suitable order for geometric progression...\n"); + _nmod_geometric_progression_init_function(G, r, ilen+olen, mod, UWORD(4)); + + prof_start(); + for (ulong ii = 0; ii < __NB_ITER; ii++) + nmod_poly_extrapolate_geometric_precomp(oval, olen, ival, ilen, ilen, G); + prof_stop(); + + nmod_geometric_progression_clear(G); + } + + flint_free(ival); + flint_free(oval); + + FLINT_TEST_CLEAR(state); +} + +/* geometric (backward), including precomp */ +void sample_geometric_backward(void * arg, ulong count) +{ + ulong n; + nmod_t mod; + ulong i; + + info_t * info = (info_t *) arg; + flint_bitcnt_t bits = info->bits; + slong ilen = info->ilen; + slong olen = info->olen; + + FLINT_TEST_INIT(state); + + nn_ptr ival = FLINT_ARRAY_ALLOC(ilen, ulong); + nn_ptr oval = FLINT_ARRAY_ALLOC(olen, ulong); + + for (i = 0; i < count; i++) + { + n = n_randprime(state, bits, 0); + nmod_init(&mod, n); + + _nmod_vec_rand(ival, state, ilen, mod); + ulong r = nmod_find_root(2*(ilen + olen), 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_extrapolate_geometric(oval, olen, ival, ilen, -olen, r, mod); + prof_stop(); + } + + flint_free(ival); + flint_free(oval); + + FLINT_TEST_CLEAR(state); +} + +/* geometric, not counting precomputations */ +void sample_geometric_backward_noprecomp(void * arg, ulong count) +{ + ulong n; + nmod_t mod; + ulong i; + + info_t * info = (info_t *) arg; + flint_bitcnt_t bits = info->bits; + slong ilen = info->ilen; + slong olen = info->olen; + + FLINT_TEST_INIT(state); + + nn_ptr ival = FLINT_ARRAY_ALLOC(ilen, ulong); + nn_ptr oval = FLINT_ARRAY_ALLOC(olen, ulong); + + for (i = 0; i < count; i++) + { + nmod_geometric_progression_t G; + n = n_randprime(state, bits, 0); + nmod_init(&mod, n); + + _nmod_vec_rand(ival, state, ilen, mod); + ulong r = nmod_find_root(2*(ilen + olen), mod); + if (r == 0) + flint_printf("\n...could not find element of suitable order for geometric progression...\n"); + _nmod_geometric_progression_init_function(G, r, ilen+olen, mod, UWORD(4)); + + prof_start(); + for (ulong ii = 0; ii < __NB_ITER; ii++) + nmod_poly_extrapolate_geometric_precomp(oval, olen, ival, ilen, -olen, G); + prof_stop(); + + nmod_geometric_progression_clear(G); + } + + flint_free(ival); + flint_free(oval); + + FLINT_TEST_CLEAR(state); +} + +/* geometric, not counting precomputations, via interpolation+evaluation */ +void sample_geometric_noprecomp_intev(void * arg, ulong count) +{ + ulong n; + nmod_t mod; + ulong i; + + info_t * info = (info_t *) arg; + flint_bitcnt_t bits = info->bits; + slong ilen = info->ilen; + slong olen = info->olen; + + FLINT_TEST_INIT(state); + + nn_ptr ival = FLINT_ARRAY_ALLOC(ilen, ulong); + nn_ptr oval = FLINT_ARRAY_ALLOC(olen, ulong); + + for (i = 0; i < count; i++) + { + nmod_geometric_progression_t G; + n = n_randprime(state, bits, 0); + nmod_init(&mod, n); + + _nmod_vec_rand(ival, state, ilen, mod); + ulong r = nmod_find_root(2*(ilen + olen), mod); + if (r == 0) + flint_printf("\n...could not find element of suitable order for geometric progression...\n"); + _nmod_geometric_progression_init_function(G, r, ilen+olen, mod, UWORD(3)); + + prof_start(); + nn_ptr poly = FLINT_ARRAY_ALLOC(ilen, ulong); + for (ulong ii = 0; ii < __NB_ITER; ii++) + { + _nmod_poly_interpolate_geometric_nmod_vec_fast_precomp(poly, ival, G, ilen, mod); + /* scale so that next evaluations are at q**j for j=ilen...ilen+olen-1 */ + ulong q_pow_ilen = nmod_pow_ui(r, 2*ilen, mod); + ulong qk = q_pow_ilen; + for (slong k = 1; k < ilen; k++) + { + poly[k] = nmod_mul(poly[k], qk, mod); + qk = nmod_mul(qk, q_pow_ilen, mod); + } + _nmod_poly_evaluate_geometric_nmod_vec_fast_precomp(oval, poly, ilen, G, olen, mod); + } + flint_free(poly); + prof_stop(); + + nmod_geometric_progression_clear(G); + } + + flint_free(ival); + flint_free(oval); + + FLINT_TEST_CLEAR(state); +} + +/* geometric, counting only precomputations (on npoints_precomp points) */ +void sample_geometric_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 ilen = info->ilen; + slong olen = info->olen; + + FLINT_TEST_INIT(state); + + for (i = 0; i < count; i++) + { + nmod_geometric_progression_t G; + n = n_randprime(state, bits, 0); + nmod_init(&mod, n); + + ulong r = nmod_find_root(2*(ilen + olen), 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_geometric_progression_init_function(G, r, ilen+olen, mod, UWORD(4)); + nmod_geometric_progression_clear(G); + } + prof_stop(); + } + + FLINT_TEST_CLEAR(state); +} + +/* comparison with relevant middle product */ +void sample_nmod_poly_mulmid(void * arg, ulong count) +{ + ulong n; + nmod_t mod; + ulong i; + + info_t * info = (info_t *) arg; + flint_bitcnt_t bits = info->bits; + slong ilen = info->ilen; + slong olen = info->olen; + + FLINT_TEST_INIT(state); + + for (i = 0; i < count; i++) + { + n = n_randprime(state, bits, 1); + nmod_init(&mod, n); + nn_ptr poly1 = FLINT_ARRAY_ALLOC(ilen+olen-1, ulong); + nn_ptr poly2 = FLINT_ARRAY_ALLOC(ilen, ulong); + nn_ptr poly3 = FLINT_ARRAY_ALLOC(olen, ulong); + + _nmod_vec_rand(poly1, state, ilen+olen-1, mod); + _nmod_vec_rand(poly2, state, ilen, mod); + + prof_start(); + for (ulong ii = 0; ii < __NB_ITER; ii++) + _nmod_poly_mulmid(poly3, poly1, ilen+olen-1, poly2, ilen, ilen-1, ilen+olen-1, mod); + prof_stop(); + + flint_free(poly1); + flint_free(poly2); + flint_free(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] [growth] [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" + " [growth] is optional (default 1.0)\n" + " positive float: extrapolate from `ilen` points to `growth * ilen` points\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 double growth = (argc >= 3) ? atof(argv[2]) : 1.; + + if (growth != 1.) + flint_printf("[growth provided] -> will extrapolate from `ilen` points to `%lf * ilen` points\n", growth); + + /* 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_general; + double time_general_noprecomp; + double time_general_onlyprecomp; + double time_geometric; + double time_geometric_noprecomp; + double time_geometric_backward; + double time_geometric_backward_noprecomp; + double time_geometric_intev; + double time_geometric_onlyprecomp; + double time_poly_mulmid; + + 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("ilen\t olen |general\tw/ tree\t tree | geom\tprec\tgeom-b\tprec-b\tint+ev |precomp\t|mulmid\n"); + else if (func_bench == 1) /* general only */ + flint_printf("ilen\t olen |general\tw/ tree\t tree |\n"); + else if (func_bench == 2) /* geometric only */ + flint_printf("ilen\t olen |geom\tprec\tgeom-b\tprec-b\tint+ev |precomp\t|mulmid\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.ilen = lengths[len]; + info.olen = growth * lengths[len]; + + if (func_bench == 0 || func_bench == 1) + { + prof_repeat(&time_general, &tmp, sample_general, (void *) &info); + + prof_repeat(&time_general_noprecomp, &tmp, sample_general_noprecomp, (void *) &info); + + prof_repeat(&time_general_onlyprecomp, &tmp, sample_general_onlyprecomp, (void *) &info); + } + + if (func_bench == 0 || func_bench == 2) + { + prof_repeat(&time_geometric, &tmp, sample_geometric, (void *) &info); + + prof_repeat(&time_geometric_noprecomp, &tmp, sample_geometric_noprecomp, (void *) &info); + + prof_repeat(&time_geometric_backward, &tmp, sample_geometric_backward, (void *) &info); + + prof_repeat(&time_geometric_backward_noprecomp, &tmp, sample_geometric_backward_noprecomp, (void *) &info); + + prof_repeat(&time_geometric_intev, &tmp, sample_geometric_noprecomp_intev, (void *) &info); + + prof_repeat(&time_geometric_onlyprecomp, &tmp, sample_geometric_onlyprecomp, (void *) &info); + + prof_repeat(&time_poly_mulmid, &tmp, sample_nmod_poly_mulmid, (void *) &info); + } + + if (func_bench == 0) + { + flint_printf("%ld\t%7ld|%.1e %.1e %.1e|%.1e %.1e %.1e %.1e %.1e|%.1e |%.1e\n", + info.ilen, info.olen, + time_general/fac, time_general_noprecomp/fac, time_general_onlyprecomp/fac, + time_geometric/fac, time_geometric_noprecomp/fac, + time_geometric_backward/fac, time_geometric_backward_noprecomp/fac, + time_geometric_intev/fac, time_geometric_onlyprecomp/fac, + time_poly_mulmid/fac); + } + else if (func_bench == 1) + { + flint_printf("%ld\t%7ld|%.1e %.1e %.1e\n", + info.ilen, info.olen, + time_general/fac, time_general_noprecomp/fac, time_general_onlyprecomp/fac); + } + else if (func_bench == 2) + { + flint_printf("%ld\t%7ld|%.1e %.1e %.1e %.1e %.1e|%.1e |%.1e\n", + info.ilen, info.olen, + time_geometric/fac, time_geometric_noprecomp/fac, + time_geometric_backward/fac, time_geometric_backward_noprecomp/fac, + time_geometric_intev/fac, time_geometric_onlyprecomp/fac, + time_poly_mulmid/fac); + } + } + + flint_printf("\n"); + } + + return 0; +} + +#undef __NB_ITER From eadb226a30f999b868b5685415a49025975cfb96 Mon Sep 17 00:00:00 2001 From: Vincent Neiger Date: Sat, 9 May 2026 19:43:06 +0200 Subject: [PATCH 3/6] Clean init functions and add function for init with both extrapolation+interpolation Co-authored-by: ktran11 --- ...ric_nmod_vec.c => extrapolate_geometric.c} | 1 - src/nmod_poly/geometric_progression.c | 130 +++++++++++++++--- .../profile/p-eval_interp_precomputations.c | 53 ++++++- 3 files changed, 155 insertions(+), 29 deletions(-) rename src/nmod_poly/{extrapolate_geometric_nmod_vec.c => extrapolate_geometric.c} (99%) diff --git a/src/nmod_poly/extrapolate_geometric_nmod_vec.c b/src/nmod_poly/extrapolate_geometric.c similarity index 99% rename from src/nmod_poly/extrapolate_geometric_nmod_vec.c rename to src/nmod_poly/extrapolate_geometric.c index 0a21f4cfa8..6fc27593b1 100644 --- a/src/nmod_poly/extrapolate_geometric_nmod_vec.c +++ b/src/nmod_poly/extrapolate_geometric.c @@ -98,7 +98,6 @@ void nmod_poly_extrapolate_geometric(nn_ptr oval, slong olen, nn_srcptr ival, slong ilen, slong offset, ulong r, nmod_t mod) { - /* TODO handle overlap */ if (ilen == 0) { for (slong i = 0; i < olen; i++) diff --git a/src/nmod_poly/geometric_progression.c b/src/nmod_poly/geometric_progression.c index 6a2b9fb952..468a02964a 100644 --- a/src/nmod_poly/geometric_progression.c +++ b/src/nmod_poly/geometric_progression.c @@ -194,10 +194,10 @@ void _nmod_geometric_progression_extrapolate_init(nmod_geometric_progression_t G G->ext_ff->length = len-1; G->ext_fb->length = len-1; - /* ext_s1f[i] = prod_{k=1}^{i} (q**k - 1) [for forward] */ /* TODO s0f */ - /* ext_s1b[i] = prod_{k=1}^{i} (q**(-k) - 1) [for backward] */ /* TODO s0b */ - /* ext_s2[i] = 1 / ext_s1f[i] [for both] */ /* TODO s1p */ - /* ext_s3[i] = 1 / ext_s1b[i] [for both] */ /* TODO s1m */ + /* ext_s1f[i] = prod_{k=1}^{i} (q**k - 1) [for forward] */ + /* ext_s1b[i] = prod_{k=1}^{i} (q**(-k) - 1) [for backward] */ + /* ext_s2[i] = 1 / ext_s1f[i] [for both] */ + /* ext_s3[i] = 1 / ext_s1b[i] [for both] */ G->ext_s1f = _nmod_vec_init(len); G->ext_s1b = _nmod_vec_init(len); G->ext_s2 = _nmod_vec_init(len); @@ -211,8 +211,8 @@ void _nmod_geometric_progression_extrapolate_init(nmod_geometric_progression_t G for (slong i = 1; i < len; i++) { - G->ext_ff->coeffs[i-1] = q_pow_i - 1; - G->ext_fb->coeffs[i-1] = inv_q_pow_i - 1; + G->ext_ff->coeffs[i-1] = q_pow_i - 1; /* temporarily, q**i - 1 */ + G->ext_fb->coeffs[i-1] = inv_q_pow_i - 1; /* temporarily, q**(-i) - 1 */ G->ext_s1f[i] = nmod_mul(G->ext_s1f[i-1], q_pow_i - 1, mod); G->ext_s1b[i] = nmod_mul(G->ext_s1b[i-1], inv_q_pow_i - 1, mod); q_pow_i = nmod_mul(q_pow_i, q, mod); @@ -231,10 +231,82 @@ void _nmod_geometric_progression_extrapolate_init(nmod_geometric_progression_t G } } -/* 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) */ +static +void _nmod_geometric_progression_interpolate_extrapolate_init(nmod_geometric_progression_t G, + slong len, nmod_t mod, + ulong q, ulong inv_q) +{ + /* ext_ff: sum_{i=0}^{len-1} x**i / (q**(i+1) - 1) [for forward] */ + /* ext_fb: sum_{i=0}^{len-1} x**i / (q**(-i-1) - 1) [for backward] */ + nmod_poly_init2_preinv(G->ext_ff, mod.n, mod.ninv, len-1); + nmod_poly_init2_preinv(G->ext_fb, mod.n, mod.ninv, len-1); + G->ext_ff->length = len-1; + G->ext_fb->length = len-1; + + /* ext_s1f[i] = prod_{k=1}^{i} (q**k - 1) [for forward] */ + /* ext_s1b[i] = prod_{k=1}^{i} (q**(-k) - 1) [for backward] */ + /* ext_s2[i] = 1 / ext_s1f[i] [for both] */ + /* ext_s3[i] = 1 / ext_s1b[i] [for both] */ + G->ext_s1f = _nmod_vec_init(len); + G->ext_s1b = _nmod_vec_init(len); + G->ext_s2 = _nmod_vec_init(len); + G->ext_s3 = _nmod_vec_init(len); + + /* coeff(G->int_f, i) = (-1)**i * q**(i*(i-1)/2) / (prod_{k=1}^{i} (q**k - 1)) */ + /* == (-1)**i * q**(i*(i-1)/2) * ext_s2[i] */ + nmod_poly_init2_preinv(G->int_f, mod.n, mod.ninv, len); + G->int_f->length = len; + + /* G->int_s1[i] = ext_s2[i] */ + /* G->int_s2[i] = ext_s1f[i] / q**(i*(i-1)/2) */ + G->int_s1 = G->ext_s2; + G->int_s2 = _nmod_vec_init(len); + + G->ext_s1f[0] = 1; + G->ext_s1b[0] = 1; + G->int_f->coeffs[0] = 1; + G->int_s2[0] = 1; + + ulong q_pow_i = 1; + ulong inv_q_pow_i = 1; + + for (slong i = 1; i < len; i++) + { + G->int_f->coeffs[i] = nmod_mul(G->int_f->coeffs[i-1], + mod.n - q_pow_i, mod); /* temporarily, (-1)**i * q**(i*(i-1)/2) */ + G->int_s2[i] = nmod_mul(G->int_s2[i-1], inv_q_pow_i, mod); /* temporarily, 1 / q**(i*(i-1)/2) */ + G->int_s2[i-1] = nmod_mul(G->int_s2[i-1], G->ext_s1f[i-1], mod); /* int_s2[i-1] is now ok */ + q_pow_i = nmod_mul(q_pow_i, q, mod); + inv_q_pow_i = nmod_mul(inv_q_pow_i, inv_q, mod); + G->ext_ff->coeffs[i-1] = q_pow_i - 1; /* temporarily, q**i - 1 */ + G->ext_fb->coeffs[i-1] = inv_q_pow_i - 1; /* temporarily, q**(-i) - 1 */ + G->ext_s1f[i] = nmod_mul(G->ext_s1f[i-1], q_pow_i - 1, mod); + G->ext_s1b[i] = nmod_mul(G->ext_s1b[i-1], inv_q_pow_i - 1, mod); + } + + G->int_s2[len-1] = nmod_mul(G->int_s2[len-1], G->ext_s1f[len-1], mod); /* int_s2[len-1] is now ok */ + G->ext_s2[len-1] = nmod_inv(G->ext_s1f[len-1], mod); + G->ext_s3[len-1] = nmod_inv(G->ext_s1b[len-1], mod); + + for (slong i = len-1; i > 0; i--) + { + G->ext_s2[i-1] = nmod_mul(G->ext_s2[i], G->ext_ff->coeffs[i-1], mod); + G->int_f->coeffs[i] = nmod_mul(G->ext_s2[i], G->int_f->coeffs[i], mod); + G->ext_s3[i-1] = nmod_mul(G->ext_s3[i], G->ext_fb->coeffs[i-1], mod); + G->ext_ff->coeffs[i-1] = nmod_mul(G->ext_s2[i], G->ext_s1f[i-1], mod); + G->ext_fb->coeffs[i-1] = nmod_mul(G->ext_s3[i], G->ext_s1b[i-1], mod); + } +} + +/* 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) */ +/* */ +/* TODO possible improvements, if precomputation efficiency matters: */ +/* - extrapolate may benefit from n_mulmod_shoup */ +/* - interpolate version with mulmod_shoup is not fully "Shoupified" */ +/* - unrolling may help */ void _nmod_geometric_progression_init_function(nmod_geometric_progression_t G, ulong r, slong len, nmod_t mod, ulong function) @@ -258,13 +330,13 @@ void _nmod_geometric_progression_init_function(nmod_geometric_progression_t G, if (function & UWORD(1)) /* evaluate */ _nmod_geometric_progression_evaluate_init_nonfullword(G, r, len, mod, q, q_pr_quo, q_pr_rem, inv_r, inv_q, inv_q_pr_quo, inv_q_pr_rem); - if ((function>>1) & UWORD(1)) /* interpolate */ + + if ((function & UWORD(6)) == UWORD(6)) /* interpolate+extrapolate */ + _nmod_geometric_progression_interpolate_extrapolate_init(G, len, mod, q, inv_q); + else if ((function & UWORD(2)) == UWORD(2)) /* 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); - /* TODO improve Shoup interpolate */ - if ((function>>2) & UWORD(1)) /* extrapolate */ + else if ((function & UWORD(4)) == UWORD(4)) /* extrapolate */ _nmod_geometric_progression_extrapolate_init(G, len, mod, q, inv_q); - /* TODO Shoup extrapolate */ - /* TODO Shoup version of combine extrapolate&interpolate */ } else { @@ -274,11 +346,13 @@ void _nmod_geometric_progression_init_function(nmod_geometric_progression_t G, if (function & UWORD(1)) /* evaluate */ _nmod_geometric_progression_evaluate_init(G, r, len, mod, q, inv_r, inv_q); - if ((function>>1) & UWORD(1)) /* interpolate */ - _nmod_geometric_progression_interpolate_init(G, len, mod, q, inv_q); - if ((function>>2) & UWORD(1)) /* extrapolate */ - _nmod_geometric_progression_extrapolate_init(G, len, mod, q, inv_q); - /* TODO combine extrapolate&interpolate (note: ext_s1f is the same as G->int_s1) */ + + if ((function & UWORD(6)) == UWORD(6)) /* interpolate+extrapolate */ + _nmod_geometric_progression_interpolate_extrapolate_init(G, len, mod, q, inv_q); + else if ((function & UWORD(2)) == UWORD(2)) /* interpolate */ + _nmod_geometric_progression_interpolate_init(G, len, mod, q, inv_q); + else if ((function & UWORD(4)) == UWORD(4)) /* extrapolate */ + _nmod_geometric_progression_extrapolate_init(G, len, mod, q, inv_q); } } @@ -291,13 +365,25 @@ void _nmod_geometric_progression_clear_function(nmod_geometric_progression_t G, _nmod_vec_clear(G->ev_s); nmod_poly_clear(G->ev_f); } - if ((function>>1) & UWORD(1)) /* interpolate */ + if ((function & UWORD(6)) == UWORD(6)) /* interpolate+extrapolate */ + { + /* G->int_s1 is ext_s2: no alloc -> not cleared */ + _nmod_vec_clear(G->int_s2); + nmod_poly_clear(G->int_f); + _nmod_vec_clear(G->ext_s1f); + _nmod_vec_clear(G->ext_s1b); + _nmod_vec_clear(G->ext_s2); + _nmod_vec_clear(G->ext_s3); + nmod_poly_clear(G->ext_ff); + nmod_poly_clear(G->ext_fb); + } + else if ((function & UWORD(2)) == UWORD(2)) /* interpolate */ { _nmod_vec_clear(G->int_s1); _nmod_vec_clear(G->int_s2); nmod_poly_clear(G->int_f); } - if ((function>>2) & UWORD(1)) /* extrapolate */ + else if ((function & UWORD(4)) == UWORD(4)) /* extrapolate */ { _nmod_vec_clear(G->ext_s1f); _nmod_vec_clear(G->ext_s1b); diff --git a/src/nmod_poly/profile/p-eval_interp_precomputations.c b/src/nmod_poly/profile/p-eval_interp_precomputations.c index d278d1723b..e8fc5d994c 100644 --- a/src/nmod_poly/profile/p-eval_interp_precomputations.c +++ b/src/nmod_poly/profile/p-eval_interp_precomputations.c @@ -238,6 +238,43 @@ void sample_geometric_precomp_extrap(void * arg, ulong count) FLINT_TEST_CLEAR(state); } +/* precomputations for geometric progression (extra+inter only) */ +void sample_geometric_precomp_extra_inter(void * arg, ulong count) +{ + ulong n; + nmod_t mod; + ulong i; + + info_t * info = (info_t *) arg; + flint_bitcnt_t bits = info->bits; + slong npoints_precomp = info->npoints_precomp; + + FLINT_TEST_INIT(state); + + for (i = 0; i < count; i++) + { + n = n_randprime(state, bits, 1); + nmod_init(&mod, n); + ulong r = nmod_find_root(2*npoints_precomp, mod); + if (r == 0) + flint_printf("\n...could not find element of suitable order for geometric progression...\n"); + + nmod_geometric_progression_t G; + G->len = npoints_precomp; + G->mod = mod; + + prof_start(); + for (ulong ii = 0; ii < __NB_ITER; ii++) + { + _nmod_geometric_progression_init_function(G, r, npoints_precomp, mod, UWORD(6)); + nmod_geometric_progression_clear(G); + } + prof_stop(); + } + + FLINT_TEST_CLEAR(state); +} + int main(int argc, char * argv[]) { if (argc > 1 && (strcmp(argv[1], "-h") == 0 || strcmp(argv[1], "--help") == 0)) @@ -271,6 +308,7 @@ int main(int argc, char * argv[]) double time_geometric_eval; double time_geometric_interp; double time_geometric_extrap; + double time_geometric_extra_inter; info_t info; flint_bitcnt_t i; @@ -284,11 +322,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 | g-interp| g-extrap \n"); + flint_printf("len\tpoints | general | geom | g-eval | g-interp| g-extrap| g-int-ext \n"); else if (func_bench == 1) /* general only */ flint_printf("len\tpoints | general\n"); else if (func_bench == 2) /* geometric only */ - flint_printf("len\tpoints | geom | g-eval | g-interp| g-extrap \n"); + flint_printf("len\tpoints | geom | g-eval | g-interp| g-extrap| g-int-ext \n"); for (int len = 0; len < nb_lens; ++len) { @@ -306,17 +344,19 @@ int main(int argc, char * argv[]) prof_repeat(&time_geometric_eval, &tmp, sample_geometric_precomp_eval, (void *) &info); prof_repeat(&time_geometric_interp, &tmp, sample_geometric_precomp_interp, (void *) &info); prof_repeat(&time_geometric_extrap, &tmp, sample_geometric_precomp_extrap, (void *) &info); + prof_repeat(&time_geometric_extra_inter, &tmp, sample_geometric_precomp_extra_inter, (void *) &info); } if (func_bench == 0) { - flint_printf("%ld\t%7ld| %.1e | %.1e | %.1e | %.1e | %.1e\n", + flint_printf("%ld\t%7ld| %.1e | %.1e | %.1e | %.1e | %.1e | %.1e\n", info.length, info.npoints_precomp, time_general/fac, time_geometric_all/fac, time_geometric_eval/fac, time_geometric_interp/fac, - time_geometric_extrap/fac); + time_geometric_extrap/fac, + time_geometric_extra_inter/fac); } else if (func_bench == 1) { @@ -326,12 +366,13 @@ int main(int argc, char * argv[]) } else if (func_bench == 2) { - flint_printf("%ld\t%7ld| %.1e | %.1e | %.1e | %.1e\n", + flint_printf("%ld\t%7ld| %.1e | %.1e | %.1e | %.1e | %.1e\n", info.length, info.npoints_precomp, time_geometric_all/fac, time_geometric_eval/fac, time_geometric_interp/fac, - time_geometric_extrap/fac); + time_geometric_extrap/fac, + time_geometric_extra_inter/fac); } } From 555c4d10c07edf8cd3bdd3285048c69a9eff434c Mon Sep 17 00:00:00 2001 From: Vincent Neiger Date: Sat, 9 May 2026 20:51:59 +0200 Subject: [PATCH 4/6] Minor fix + documentation Co-authored-by: ktran11 --- doc/source/nmod_poly.rst | 43 +++++++++++++++++++ src/nmod_poly/extrapolate_geometric.c | 2 +- .../interpolate_geometric_nmod_vec.c | 4 +- 3 files changed, 46 insertions(+), 3 deletions(-) diff --git a/doc/source/nmod_poly.rst b/doc/source/nmod_poly.rst index e7c7230ba9..f938bf6ed7 100644 --- a/doc/source/nmod_poly.rst +++ b/doc/source/nmod_poly.rst @@ -1596,6 +1596,49 @@ Interpolation Uses fast geometric multipoint interpolation, building a temporary geometric progression precomputation. +Extrapolation +-------------------------------------------------------------------------------- + +.. function:: void nmod_poly_extrapolate_geometric_precomp(nn_ptr oval, slong olen, nn_srcptr ival, slong ilen, slong offset, const nmod_geometric_progression_t G) + void nmod_poly_extrapolate_geometric(nn_ptr oval, slong olen, nn_srcptr ival, slong ilen, slong offset, ulong r, nmod_t mod) + + This extrapolates the ``ilen`` input values ``ival`` to compute ``olen`` + output values ``oval``, based on points that are powers of `r^2` in + geometric progression; ``offset`` specifies the output points relative to + the input ones. + + More precisely, let `m` stand for ``ilen``, let `n` stand for ``olen``, and + let `c` be any integer that is invertible modulo ``mod.n``. In this + description we assume `r` has "large enough" order; more details are given + below. The input ``ival`` is interpreted as the list of values `(f(c \cdot + r^{2(k+i)}))_{0 \le i < m}` of some polynomial `f` of degree less than `m`, + evaluated on the subsequence `(c \cdot r^{2(k+i)})_{0 \le i < m}` of the + geometric progression, for some given `k \ge 0`. Then the output ``oval`` + is the list of values `f(c \cdot r^{2(\ell+j)})_{0 \le j < n}`, for the + starting index `\ell = k + \texttt{offset}`. + + The input constraints are as follows. The algorithm does not need to know + about `c`, nor about `k` and `\ell` except for their difference + `\texttt{offset} = \ell - k`. The value `r` should be reduced modulo + ``mod.n``. There are two situations: forward extrapolation (when `\ell > + k`, i.e., ``offset`` is positive) and backward extrapolation (when `\ell < + k`, i.e., ``offset`` is negative). + + - In the forward case, `r` should have sufficient multiplicative order so + that none of the first ``offset + olen`` powers of `r^2` is one, and + one should have ``offset >= ilen``. The latter means that the input and + output lists of points are disjoint (since `\ell \ge k+m`). + + - In the backward case, `r` should have sufficient multiplicative order so + that none of the first ``ilen - offset`` powers of `r^2` is one, and one + should have ``offset + olen <= 0`` (recall that here ``offset`` is + negative). The latter means that the input and output lists of points are + disjoint (since `\ell+n \le k`). + + The function without ``_precomp`` builds a temporary geometric progression + precomputation relative to ``r`` and ``mod``, and calls the version with + ``_precomp`` with this additional data. + Composition -------------------------------------------------------------------------------- diff --git a/src/nmod_poly/extrapolate_geometric.c b/src/nmod_poly/extrapolate_geometric.c index 6fc27593b1..6f7c324bbb 100644 --- a/src/nmod_poly/extrapolate_geometric.c +++ b/src/nmod_poly/extrapolate_geometric.c @@ -18,7 +18,7 @@ void nmod_poly_extrapolate_geometric_precomp(nn_ptr oval, slong olen, const nmod_geometric_progression_t G) { /* precomputation has been done */ - FLINT_ASSERT((G->function >> 2) & 1); + FLINT_ASSERT((G->function & UWORD(4)) == UWORD(4)); /* input/output points are disjoint, and stay within precomputed data length */ FLINT_ASSERT((offset >= ilen && G->len >= offset+olen) || (offset <= -olen && G->len >= ilen-offset)); diff --git a/src/nmod_poly/interpolate_geometric_nmod_vec.c b/src/nmod_poly/interpolate_geometric_nmod_vec.c index 9f80619894..6fe60fc4de 100644 --- a/src/nmod_poly/interpolate_geometric_nmod_vec.c +++ b/src/nmod_poly/interpolate_geometric_nmod_vec.c @@ -18,7 +18,7 @@ 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); + FLINT_ASSERT((G->function & UWORD(2)) == UWORD(2)); if (len == 1) { @@ -131,7 +131,7 @@ void _nmod_poly_interpolate_geometric_nmod_vec_fast_precomp(nn_ptr poly, 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); + FLINT_ASSERT((G->function & UWORD(2)) == UWORD(2)); nmod_poly_fit_length(poly, len); _nmod_poly_set_length(poly, len); From 520aab7921fc5252c718b0700d35d54f96a32449 Mon Sep 17 00:00:00 2001 From: Vincent Neiger Date: Sat, 9 May 2026 21:21:22 +0200 Subject: [PATCH 5/6] small fix in documentation --- doc/source/nmod_poly.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/doc/source/nmod_poly.rst b/doc/source/nmod_poly.rst index f938bf6ed7..89ed733495 100644 --- a/doc/source/nmod_poly.rst +++ b/doc/source/nmod_poly.rst @@ -1614,7 +1614,7 @@ Extrapolation r^{2(k+i)}))_{0 \le i < m}` of some polynomial `f` of degree less than `m`, evaluated on the subsequence `(c \cdot r^{2(k+i)})_{0 \le i < m}` of the geometric progression, for some given `k \ge 0`. Then the output ``oval`` - is the list of values `f(c \cdot r^{2(\ell+j)})_{0 \le j < n}`, for the + is the list of values `(f(c \cdot r^{2(\ell+j)}))_{0 \le j < n}`, for the starting index `\ell = k + \texttt{offset}`. The input constraints are as follows. The algorithm does not need to know From 3d5f466a597cfcee29889d7d358cce8e328931ff Mon Sep 17 00:00:00 2001 From: Vincent Neiger Date: Sat, 9 May 2026 22:24:56 +0200 Subject: [PATCH 6/6] fix corner cases with 0 points in input or output --- src/nmod_poly/extrapolate_geometric.c | 6 +++ src/nmod_poly/test/t-extrapolate_geometric.c | 50 +++++++++++++------- 2 files changed, 38 insertions(+), 18 deletions(-) diff --git a/src/nmod_poly/extrapolate_geometric.c b/src/nmod_poly/extrapolate_geometric.c index 6f7c324bbb..5aca851eab 100644 --- a/src/nmod_poly/extrapolate_geometric.c +++ b/src/nmod_poly/extrapolate_geometric.c @@ -23,6 +23,9 @@ void nmod_poly_extrapolate_geometric_precomp(nn_ptr oval, slong olen, FLINT_ASSERT((offset >= ilen && G->len >= offset+olen) || (offset <= -olen && G->len >= ilen-offset)); + if (olen == 0) + return; + if (ilen == 0) { for (slong i = 0; i < olen; i++) @@ -98,6 +101,9 @@ void nmod_poly_extrapolate_geometric(nn_ptr oval, slong olen, nn_srcptr ival, slong ilen, slong offset, ulong r, nmod_t mod) { + if (olen == 0) + return; + if (ilen == 0) { for (slong i = 0; i < olen; i++) diff --git a/src/nmod_poly/test/t-extrapolate_geometric.c b/src/nmod_poly/test/t-extrapolate_geometric.c index 37e59e86a2..8bb22de34d 100644 --- a/src/nmod_poly/test/t-extrapolate_geometric.c +++ b/src/nmod_poly/test/t-extrapolate_geometric.c @@ -55,18 +55,25 @@ TEST_FUNCTION_START(nmod_poly_extrapolate_geometric, state) _nmod_vec_randtest(ival, state, ilen, mod); - /* input points are c * q**k, k = istart...istart+ilen-1 */ ulong c = 1 + n_randint(state, modn-1); - ipts[0] = nmod_pow_ui(q, istart, mod); - ipts[0] = nmod_mul(ipts[0], c, mod); - for (slong k = 1; k < ilen; k++) - ipts[k] = nmod_mul(ipts[k-1], q, mod); + + /* input points are c * q**k, k = istart...istart+ilen-1 */ + if (ilen > 0) + { + ipts[0] = nmod_pow_ui(q, istart, mod); + ipts[0] = nmod_mul(ipts[0], c, mod); + for (slong k = 1; k < ilen; k++) + ipts[k] = nmod_mul(ipts[k-1], q, mod); + } /* output points are c * q**k, k = istart+offset...istart+offset+olen-1 */ - opts[0] = nmod_pow_ui(q, istart+offset, mod); - opts[0] = nmod_mul(opts[0], c, mod); - for (slong k = 1; k < olen; k++) - opts[k] = nmod_mul(opts[k-1], q, mod); + if (olen > 0) + { + opts[0] = nmod_pow_ui(q, istart+offset, mod); + opts[0] = nmod_mul(opts[0], c, mod); + for (slong k = 1; k < olen; k++) + opts[k] = nmod_mul(opts[k-1], q, mod); + } nmod_poly_interpolate_nmod_vec(P, ipts, ival, ilen); nmod_poly_evaluate_nmod_vec(val, P, opts, olen); @@ -109,18 +116,25 @@ TEST_FUNCTION_START(nmod_poly_extrapolate_geometric, state) _nmod_vec_randtest(ival, state, ilen, mod); - /* input points are c * q**k, k = istart...istart+ilen-1 */ ulong c = 1 + n_randint(state, modn-1); - ipts[0] = nmod_pow_ui(q, istart, mod); - ipts[0] = nmod_mul(ipts[0], c, mod); - for (slong k = 1; k < ilen; k++) - ipts[k] = nmod_mul(ipts[k-1], q, mod); + + /* input points are c * q**k, k = istart...istart+ilen-1 */ + if (ilen > 0) + { + ipts[0] = nmod_pow_ui(q, istart, mod); + ipts[0] = nmod_mul(ipts[0], c, mod); + for (slong k = 1; k < ilen; k++) + ipts[k] = nmod_mul(ipts[k-1], q, mod); + } /* output points are c * q**k, k = istart+offset...istart+offset+olen-1 */ - opts[0] = nmod_pow_ui(q, istart+offset, mod); - opts[0] = nmod_mul(opts[0], c, mod); - for (slong k = 1; k < olen; k++) - opts[k] = nmod_mul(opts[k-1], q, mod); + if (olen > 0) + { + opts[0] = nmod_pow_ui(q, istart+offset, mod); + opts[0] = nmod_mul(opts[0], c, mod); + for (slong k = 1; k < olen; k++) + opts[k] = nmod_mul(opts[k-1], q, mod); + } nmod_poly_interpolate_nmod_vec(P, ipts, ival, ilen); nmod_poly_evaluate_nmod_vec(val, P, opts, olen);