Extrapolation at points geometric progression, over nmod

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
 --------------------------------------------------------------------------------
 

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);

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 <https://www.gnu.org/licenses/>.
+*/
+
+#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);
+}