Port PML's geometric evaluation / interpolation integration for univariate polynomial

src/nmod_poly.h

@@ -569,6 +569,24 @@ 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  *************************************/
+
+void nmod_geometric_progression_init(nmod_geometric_progression_t G, ulong r, slong n, nmod_t mod);
+
+void nmod_geometric_progression_clear(nmod_geometric_progression_t G);
+
+void _nmod_poly_evaluate_geometric_nmod_vec_iter(nn_ptr ys, nn_srcptr coeffs, slong len, ulong r, slong n, nmod_t mod);
+void nmod_poly_evaluate_geometric_nmod_vec_iter(nn_ptr ys, const nmod_poly_t poly, ulong r, slong n);
+
+void _nmod_poly_evaluate_geometric_nmod_vec_fast_precomp(nn_ptr vs, nn_srcptr poly, slong plen, const nmod_geometric_progression_t G, slong len);
+
+void _nmod_poly_evaluate_geometric_nmod_vec_fast(nn_ptr ys, nn_srcptr coeffs, slong len, ulong r, slong n, nmod_t mod);
+void nmod_poly_evaluate_geometric_nmod_vec_fast(nn_ptr ys, const nmod_poly_t poly, ulong r, slong n);
+
+void nmod_poly_interpolate_geometric_nmod_vec_fast(nmod_poly_t poly, ulong r, nn_srcptr ys, slong n);
+void nmod_poly_interpolate_geometric_nmod_vec_fast_precomp(nmod_poly_t poly, nn_srcptr v,
+    const nmod_geometric_progression_t G, slong len);
+
 /* Interpolation  ************************************************************/
 
 void _nmod_poly_interpolate_nmod_vec_newton(nn_ptr poly, nn_srcptr xs, nn_srcptr ys, slong n, nmod_t mod);

src/nmod_poly/evaluate_geometric_nmod_vec.c

@@ -0,0 +1,97 @@
+/*
+    Copyright (C) 2025, Vincent Neiger, Éric Schost
+    Copyright (C) 2025, Mael Hostettler
+
+    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"
+#include "nmod_vec.h"
+
+void
+_nmod_poly_evaluate_geometric_nmod_vec_fast_precomp(nn_ptr vs, nn_srcptr poly, 
+    slong plen, const nmod_geometric_progression_t G, slong len)
+{
+    nmod_poly_t a, b;
+    slong i, d;
+
+    d = G->d;
+    FLINT_ASSERT(len <= d);
+
+    if (plen == 0)
+    {
+        for (i = 0; i < d; i++)
+        {
+            vs[i] = 0;
+        }
+        return;
+    }
+
+    nmod_poly_init2(a, G->mod.n, d);
+    nmod_poly_init(b, G->mod.n);
+
+    for (i = 0; i < plen; i++)
+    {
+        nmod_poly_set_coeff_ui(a, d - 1 - i, nmod_mul(G->x[i], poly[i], G->mod));
+    }
+
+    _nmod_vec_zero(a->coeffs, d - plen + 1); // TODO: better approch mentionned in https://github.com/flintlib/flint/pull/2449#discussion_r2476772717
+    nmod_poly_mul(b, a, G->f);
+
+    for (i = 0; i < len; i++)
+    {
+        vs[i] = nmod_mul(G->x[i], nmod_poly_get_coeff_ui(b, i+d-1), G->mod);
+    }
+
+    nmod_poly_clear(b);
+    nmod_poly_clear(a);
+}
+
+void 
+_nmod_poly_evaluate_geometric_nmod_vec_fast(nn_ptr ys, nn_srcptr poly, 
+    slong plen, ulong r, slong n, nmod_t mod)
+{
+    nmod_geometric_progression_t G;
+    nmod_geometric_progression_init(G, r, n, mod);
+
+    _nmod_poly_evaluate_geometric_nmod_vec_fast_precomp(ys, poly, plen, G, n);
+    nmod_geometric_progression_clear(G);
+}
+
+void
+nmod_poly_evaluate_geometric_nmod_vec_fast(nn_ptr ys,
+        const nmod_poly_t poly, ulong r, slong n)
+{
+    _nmod_poly_evaluate_geometric_nmod_vec_fast(ys, poly->coeffs,
+                                                poly->length, r, n, poly->mod);
+}
+
+void
+_nmod_poly_evaluate_geometric_nmod_vec_iter(nn_ptr ys, nn_srcptr coeffs, slong len,
+    ulong r, slong n, nmod_t mod)
+{
+    slong i;
+    ulong rpow = 1;
+    ulong r2 = nmod_mul(r, r, mod);
+
+    for (i = 0; i < n; i++)
+    {
+        ys[i] = _nmod_poly_evaluate_nmod(coeffs, len, rpow, mod);
+        rpow = nmod_mul(rpow, r2, mod);
+    }
+}
+
+void
+nmod_poly_evaluate_geometric_nmod_vec_iter(nn_ptr ys,
+    const nmod_poly_t poly, ulong r, slong n)
+{
+    _nmod_poly_evaluate_geometric_nmod_vec_iter(ys, poly->coeffs,
+                                        poly->length, r, n, poly->mod);
+}
+

src/nmod_poly/geometric_progression.c

@@ -0,0 +1,133 @@
+/*
+    Copyright (C) 2025, Vincent Neiger, Éric Schost
+    Copyright (C) 2025, Mael Hostettler
+
+    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_vec.h"
+#include "nmod_poly.h"
+
+void 
+nmod_geometric_progression_init(nmod_geometric_progression_t G, ulong r, slong d, nmod_t mod)
+{
+    ulong q, inv_r, inv_q, tmp, qk, inv_qk, qq, s;
+    nn_ptr diff, inv_diff, prod_diff;
+    slong i;
+
+    G->d = d;
+    G->mod = mod;
+
+    nmod_poly_init2(G->f, mod.n, 2*d-1);
+    nmod_poly_init(G->g1, mod.n);
+    nmod_poly_init(G->g2, mod.n);
+    nmod_poly_set_coeff_ui(G->g1, 0, 1);
+    nmod_poly_set_coeff_ui(G->g2, 0, 1);
+
+    G->x = _nmod_vec_init(d);
+    G->w = _nmod_vec_init(d);
+    G->z = _nmod_vec_init(d);
+    G->y = _nmod_vec_init(d);
+
+    G->x[0] = 1;
+    G->y[0] = 1;
+    G->w[0] = 1;
+    G->z[0] = 1;
+
+    q = nmod_mul(r, r, mod);
+    inv_r = nmod_inv(r, mod);
+    inv_q = nmod_mul(inv_r, inv_r, mod);
+
+    nmod_poly_set_coeff_ui(G->f, 0, 1);
+    tmp = r;
+    for (i = 1; i < 2*d-1; i++)
+    {
+        nmod_poly_set_coeff_ui(G->f, i, nmod_mul(nmod_poly_get_coeff_ui(G->f, i-1), tmp, mod));
+        tmp = nmod_mul(tmp, q, mod);
+    }
+
+    tmp = inv_r;
+    for (i = 1; i < d; i++)
+    {
+        G->x[i] = nmod_mul(G->x[i-1], tmp, mod);
+        tmp = nmod_mul(tmp, inv_q, mod);
+    }
+
+    inv_diff = _nmod_vec_init(d);
+    diff = _nmod_vec_init(d);
+    prod_diff = _nmod_vec_init(d);
+    inv_diff[0] = 1;
+    diff[0] = 1;
+    prod_diff[0] = 1;
+
+    qk = q;  // montgomery inversion
+    for (i = 1; i < d; i++)
+    {
+        diff[i] = qk - 1;
+        inv_diff[i] = diff[i];
+        qk = nmod_mul(qk, q, mod);
+        prod_diff[i] = nmod_mul(diff[i], prod_diff[i-1], mod);
+    }
+
+    tmp = n_invmod(prod_diff[d-1], mod.n);
+    for (i = d-1; i > 0; i--)
+    {
+        inv_diff[i] = nmod_mul(prod_diff[i-1], tmp, mod);
+        tmp = nmod_mul(tmp, diff[i], mod);
+    }
+    inv_diff[0] = tmp;
+    // end montgomery inversion
+
+    // sets sequences w, y, z and polynomials g1, g2
+    qk = 1;
+    inv_qk = 1;
+    qq = 1;
+    s = 1;
+
+    for (i = 1; i < d; i++)
+    {
+        qq = nmod_mul(qq, qk, mod);   // prod q^i
+        s = nmod_mul(s, inv_qk, mod); // prod 1/q^i
+        G->w[i] = nmod_mul(G->w[i-1], inv_diff[i], mod); // prod 1/(q^i-1)
+        tmp = nmod_mul(qq, G->w[i], mod); // prod q^i/(q^i-1)
+        nmod_poly_set_coeff_ui(G->g2, i, tmp);
+
+        if ((i & 1) == 1)
+        {
+            nmod_poly_set_coeff_ui(G->g1, i, nmod_neg(tmp, mod));
+            G->y[i] = nmod_neg(prod_diff[i], mod);
+            G->z[i] = nmod_neg(G->w[i], mod);
+        }
+        else
+        {
+            nmod_poly_set_coeff_ui(G->g1, i, tmp);
+            G->y[i] = prod_diff[i];
+            G->z[i] = G->w[i];
+        }
+        G->y[i] = nmod_mul(G->y[i], s, mod);
+
+        qk = nmod_mul(qk, q, mod);
+        inv_qk = nmod_mul(inv_qk, inv_q, mod);
+    }
+
+    _nmod_vec_clear(prod_diff);
+    _nmod_vec_clear(inv_diff);
+    _nmod_vec_clear(diff);
+}
+
+void 
+nmod_geometric_progression_clear(nmod_geometric_progression_t G)
+{
+    nmod_poly_clear(G->f);
+    nmod_poly_clear(G->g2);
+    nmod_poly_clear(G->g1);
+    _nmod_vec_clear(G->x);
+    _nmod_vec_clear(G->z);
+    _nmod_vec_clear(G->y);
+    _nmod_vec_clear(G->w);
+}

src/nmod_poly/interpolate_geometric_nmod_vec.c

@@ -0,0 +1,71 @@
+/*
+    Copyright (C) 2025, Vincent Neiger
+    Copyright (C) 2025, Mael Hostettler
+
+    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_interpolate_geometric_nmod_vec_fast_precomp(nmod_poly_t poly, nn_srcptr v,
+    const nmod_geometric_progression_t G, slong len)
+{
+    slong i, N;
+    nmod_poly_t f, h;
+    nmod_t mod;
+
+    N = G->d;
+    FLINT_ASSERT(len == N);
+    nmod_poly_fit_length(poly, N);
+    nmod_poly_zero(poly);
+
+    if (N == 1)
+    {
+        nmod_poly_set_coeff_ui(poly, 0, v[0]);
+        return;
+    }
+
+    mod = G->mod;
+    nmod_poly_init2(f, mod.n, N);
+    nmod_poly_init2(h, mod.n, N);
+
+    for(i = 0; i < N; i++)
+    {
+        nmod_poly_set_coeff_ui(f, i, nmod_mul(v[i], G->w[i], mod));
+    }
+
+    nmod_poly_mullow(h, f, G->g1, N);
+
+    for (i = 0; i < N; i++)
+    {
+        nmod_poly_set_coeff_ui(f, N - 1 - i, nmod_mul(nmod_poly_get_coeff_ui(h, i), G->y[i], mod));
+    }
+
+    nmod_poly_mullow(h, f, G->g2, N);
+
+    for (i = 0; i < N; i++)
+    {
+        nmod_poly_set_coeff_ui(poly, i, nmod_mul(nmod_poly_get_coeff_ui(h, N - 1 - i), G->z[i], mod));
+    }
+
+    nmod_poly_clear(f);
+    nmod_poly_clear(h);
+}
+
+void
+nmod_poly_interpolate_geometric_nmod_vec_fast(nmod_poly_t poly, ulong r, 
+    nn_srcptr ys, slong n)
+{
+    nmod_geometric_progression_t G;
+    nmod_geometric_progression_init(G, r, n, poly->mod);
+
+    nmod_poly_interpolate_geometric_nmod_vec_fast_precomp(poly, ys, G, n);
+    nmod_geometric_progression_clear(G);
+}

src/nmod_poly/test/main.c

@@ -52,6 +52,7 @@
 #include "t-evaluate_mat_paterson_stockmeyer.c"
 #include "t-evaluate_nmod.c"
 #include "t-evaluate_nmod_vec_fast.c"
+#include "t-evaluate_geometric_nmod_vec_fast.c"
 #include "t-exp_series.c"
 #include "t-find_distinct_nonzero_roots.c"
 #include "t-fread_print.c"
@@ -178,6 +179,7 @@ test_struct tests[] =
     TEST_FUNCTION(nmod_poly_evaluate_mat_paterson_stockmeyer),
     TEST_FUNCTION(nmod_poly_evaluate_nmod),
     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_find_distinct_nonzero_roots),
     TEST_FUNCTION(nmod_poly_fread_print),