diff --git a/doc/source/nmod_poly.rst b/doc/source/nmod_poly.rst
index e7c7230ba9..89ed733495 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.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.c b/src/nmod_poly/extrapolate_geometric.c
new file mode 100644
index 0000000000..5aca851eab
--- /dev/null
+++ b/src/nmod_poly/extrapolate_geometric.c
@@ -0,0 +1,126 @@
+/*
+ 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 & 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));
+
+ if (olen == 0)
+ return;
+
+ 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)
+ {
+ TMP_INIT;
+ TMP_START;
+ nn_ptr tmp = TMP_ALLOC(ilen * sizeof(ulong));
+
+ /* first scaling */
+ for (slong i = 0; i < ilen; i++)
+ {
+ 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 */
+ _nmod_poly_mulmid(oval, G->ext_ff->coeffs + offset-ilen, ilen+olen-1, tmp, ilen, ilen-1, ilen+olen-1, G->mod);
+
+ /* second scaling */
+ 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);
+ }
+
+ TMP_END;
+ }
+
+ /* backward extrapolation */
+ else
+ {
+ TMP_INIT;
+ TMP_START;
+ nn_ptr tmp = TMP_ALLOC(FLINT_MAX(ilen, olen) * sizeof(ulong));
+
+ /* first scaling */
+ for (slong i = 0; i < ilen; i++)
+ {
+ 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 */
+ _nmod_poly_mulmid(oval, G->ext_fb->coeffs - (offset+olen), ilen+olen-1, tmp, ilen, ilen-1, ilen+olen-1, G->mod);
+
+ /* second scaling */
+ 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);
+ }
+
+ TMP_END;
+ }
+}
+
+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++)
+ 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..468a02964a 100644
--- a/src/nmod_poly/geometric_progression.c
+++ b/src/nmod_poly/geometric_progression.c
@@ -182,10 +182,131 @@ void _nmod_geometric_progression_interpolate_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_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);
+
+ 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; /* 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);
+ 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);
+ }
+}
+
+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)
@@ -209,10 +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);
- /* 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);
}
else
{
@@ -222,10 +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 */
+
+ 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);
}
}
@@ -238,14 +365,33 @@ 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);
+ _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/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);
diff --git a/src/nmod_poly/profile/p-eval_interp_precomputations.c b/src/nmod_poly/profile/p-eval_interp_precomputations.c
index d5571f3479..e8fc5d994c 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,80 @@ 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);
+}
+
+/* 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))
@@ -233,6 +307,8 @@ int main(int argc, char * argv[])
double time_geometric_all;
double time_geometric_eval;
double time_geometric_interp;
+ double time_geometric_extrap;
+ double time_geometric_extra_inter;
info_t info;
flint_bitcnt_t i;
@@ -246,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 \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 \n");
+ flint_printf("len\tpoints | geom | g-eval | g-interp| g-extrap| g-int-ext \n");
for (int len = 0; len < nb_lens; ++len)
{
@@ -267,16 +343,20 @@ 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);
+ 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\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_interp/fac,
+ time_geometric_extrap/fac,
+ time_geometric_extra_inter/fac);
}
else if (func_bench == 1)
{
@@ -286,11 +366,13 @@ 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 | %.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,
+ time_geometric_extra_inter/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
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..8bb22de34d
--- /dev/null
+++ b/src/nmod_poly/test/t-extrapolate_geometric.c
@@ -0,0 +1,171 @@
+/*
+ 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);
+
+ ulong c = 1 + n_randint(state, modn-1);
+
+ /* 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 */
+ 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);
+
+ 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);
+
+ ulong c = 1 + n_randint(state, modn-1);
+
+ /* 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 */
+ 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);
+
+ 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);
+}