Skip to content

Speed up nmod_poly_evaluate_nmod asymptotically #2752

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
fredrik-johansson merged 1 commit into from
Jun 24, 2026
Merged

Speed up nmod_poly_evaluate_nmod asymptotically #2752

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
21 changes: 11 additions & 10 deletions doc/source/nmod_poly.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1369,22 +1369,23 @@ Evaluation
multiplications done as in :func:`n_mulmod_shoup` using the precomputed
``c_precomp`` obtained via :func:`n_mulmod_precomp_shoup`.

.. function:: ulong _nmod_poly_evaluate_nmod(nn_srcptr poly, slong len, ulong c, nmod_t mod)
.. function:: ulong _nmod_poly_evaluate_nmod_horner(nn_srcptr poly, slong len, ulong c, nmod_t mod)
ulong _nmod_poly_evaluate_nmod_rectangular(nn_srcptr poly, slong len, ulong c, nmod_t mod)
ulong _nmod_poly_evaluate_nmod(nn_srcptr poly, slong len, ulong c, nmod_t mod)

Evaluates ``poly`` at the value ``c`` and reduces modulo the given modulus
of ``poly``. The value ``c`` should be reduced modulo the modulus. The
algorithm used is Horner's method, with multiplications done via
:func:`nmod_mul`.
of ``poly``. The value ``c`` should be reduced modulo the modulus.
The implemented algorithms are Horner's method and rectangular splitting.
The default function :func:`_nmod_poly_evaluate_nmod` selects
automatically between :func:`_nmod_poly_evaluate_nmod_horner`,
:func:`_nmod_poly_evaluate_nmod_rectangular`
:func:`_nmod_poly_evaluate_nmod_precomp`, and
:func:`_nmod_poly_evaluate_nmod_precomp_lazy`.

.. function:: ulong nmod_poly_evaluate_nmod(const nmod_poly_t poly, ulong c)

Evaluates ``poly`` at the value ``c`` and reduces modulo the modulus of
``poly``. The value ``c`` should be reduced modulo the modulus. The
algorithm used is Horner's method, with multiplications and additions done
differently depending on the modulus ``poly->mod`` and on the degree (calls
one of :func:`_nmod_poly_evaluate_nmod`,
:func:`_nmod_poly_evaluate_nmod_precomp`,
:func:`_nmod_poly_evaluate_nmod_precomp_lazy`).
``poly``. The value ``c`` should be reduced modulo the modulus.

.. function:: void nmod_poly_evaluate_mat_horner(nmod_mat_t dest, const nmod_poly_t poly, const nmod_mat_t c)

Expand Down
2 changes: 2 additions & 0 deletions src/nmod_poly.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -554,6 +554,8 @@ void nmod_poly_integral(nmod_poly_t x_int, const nmod_poly_t x);

/* Evaluation ***************************************************************/

ulong _nmod_poly_evaluate_nmod_horner(nn_srcptr poly, slong len, ulong c, nmod_t mod);
ulong _nmod_poly_evaluate_nmod_rectangular(nn_srcptr poly, slong len, ulong c, nmod_t mod);
ulong _nmod_poly_evaluate_nmod(nn_srcptr poly, slong len, ulong c, nmod_t mod);
ulong _nmod_poly_evaluate_nmod_precomp(nn_srcptr poly, slong len, ulong c, ulong c_precomp, ulong modn);
ulong _nmod_poly_evaluate_nmod_precomp_lazy(nn_srcptr poly, slong len, ulong c, ulong c_precomp, ulong modn);
Expand Down
170 changes: 131 additions & 39 deletions src/nmod_poly/evaluate_nmod.c
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,6 +1,7 @@
/*
Copyright (C) 2010 William Hart
Copyright (C) 2024 Vincent Neiger
Copyright (C) 2026 Fredrik Johansson

This file is part of FLINT.

Expand All @@ -12,11 +13,96 @@

#include "flint-mparam.h"
#include "ulong_extras.h"
#include "nmod_vec.h"
#include "nmod_poly.h"
#include "nmod.h"

static void
_nmod_vec_set_powers(nn_ptr res, ulong x, slong len, nmod_t mod)
{
slong i;
res[0] = 1;
res[1] = x;

for (i = 2; i < len; i++)
res[i] = nmod_mul(res[i - 1], x, mod);
}

static void
_nmod_vec_set_powers_precomp(nn_ptr res, ulong x, slong len, nmod_t mod)
{
slong i;
res[0] = 1;
res[1] = x;
ulong xpre = n_mulmod_precomp_shoup(x, mod.n);

for (i = 2; i < len; i++)
res[i] = n_mulmod_shoup(x, res[i - 1], xpre, mod.n);
}

ulong
_nmod_poly_evaluate_nmod(nn_srcptr poly, slong len, ulong c, nmod_t mod)
_nmod_poly_evaluate_nmod_rectangular(nn_srcptr poly,
slong len, ulong x, nmod_t mod)
{
slong i, m, r;
nn_ptr xs;
ulong s, y, xmpre;
TMP_INIT;

if (len == 0)
return 0;
if (len == 1 || x == 0)
return poly[0];
if (len == 2)
return nmod_addmul(poly[0], poly[1], x, mod);

m = n_sqrt(len) + 1;
m = FLINT_MIN(m, 1024);
r = (len + m - 1) / m;

dot_params_t params = _nmod_vec_dot_params(m - 1, mod);

TMP_START;
xs = TMP_ALLOC((m + 1) * sizeof(ulong));

if (NMOD_CAN_USE_SHOUP(mod))
{
_nmod_vec_set_powers_precomp(xs, x, m + 1, mod);
xmpre = n_mulmod_precomp_shoup(xs[m], mod.n);

y = _nmod_vec_dot(poly + (r - 1) * m + 1, xs + 1, len - (r - 1) * m - 1, mod, params);
y = nmod_add(y, poly[(r - 1) * m], mod);

for (i = r - 2; i >= 0; i--)
{
s = _nmod_vec_dot(poly + i * m + 1, xs + 1, m - 1, mod, params);
s = nmod_add(s, poly[i * m], mod);
y = n_mulmod_shoup(xs[m], y, xmpre, mod.n);
y = nmod_add(y, s, mod);
}
}
else
{
_nmod_vec_set_powers(xs, x, m + 1, mod);
y = _nmod_vec_dot(poly + (r - 1) * m + 1, xs + 1, len - (r - 1) * m - 1, mod, params);
y = nmod_add(y, poly[(r - 1) * m], mod);

for (i = r - 2; i >= 0; i--)
{
s = _nmod_vec_dot(poly + i * m + 1, xs + 1, m - 1, mod, params);
s = nmod_add(s, poly[i * m], mod);
y = nmod_mul(xs[m], y, mod);
y = nmod_add(y, s, mod);
}
}

TMP_END;

return y;
}

ulong
_nmod_poly_evaluate_nmod_horner(nn_srcptr poly, slong len, ulong c, nmod_t mod)
{
slong m;
ulong val;
Expand Down Expand Up @@ -98,53 +184,59 @@ _nmod_poly_evaluate_nmod_precomp_lazy(nn_srcptr poly, slong len, ulong c, ulong
return val;
}

/* if 3*mod.n - 1 <= 2**FLINT_BITS, can use the lazy variant */
#if FLINT_BITS == 64
#define LAZY_MAX UWORD(6148914691236517205)
#else
#define LAZY_MAX UWORD(1431655765)
#endif

ulong
nmod_poly_evaluate_nmod(const nmod_poly_t poly, ulong c)
_nmod_poly_evaluate_nmod(nn_srcptr poly, slong len, ulong c, nmod_t mod)
{
if (poly->length == 0)
if (len == 0)
return 0;

if (poly->length == 1 || c == 0)
return poly->coeffs[0];
if (len == 1 || c == 0)
return poly[0];

// if degree below the n_mulmod_shoup threshold
// or modulus forbids n_mulmod_shoup usage, use nmod_mul
#if FLINT_MULMOD_SHOUP_THRESHOLD <= 2
if (poly->mod.norm == 0) // here poly->length >= threshold
#else
if ((poly->length < FLINT_MULMOD_SHOUP_THRESHOLD)
|| (poly->mod.norm == 0))
#endif
if (!NMOD_CAN_USE_SHOUP(mod))
{
return _nmod_poly_evaluate_nmod(poly->coeffs, poly->length, c, poly->mod);
if (len <= 13)
return _nmod_poly_evaluate_nmod_horner(poly, len, c, mod);
else
return _nmod_poly_evaluate_nmod_rectangular(poly, len, c, mod);
}

// if 3*mod.n - 1 <= 2**FLINT_BITS, use n_mulmod_shoup, lazy variant
else
{
const ulong modn = poly->mod.n;
if (len < FLINT_MULMOD_SHOUP_THRESHOLD)
return _nmod_poly_evaluate_nmod_horner(poly, len, c, mod);

#if FLINT_BITS == 64
if (modn <= UWORD(6148914691236517205))
#else // FLINT_BITS == 32
if (modn <= UWORD(1431655765))
#endif
{
const ulong c_precomp = n_mulmod_precomp_shoup(c, modn);
ulong val = _nmod_poly_evaluate_nmod_precomp_lazy(poly->coeffs, poly->length, c, c_precomp, modn);
// correct excess
if (val >= 2*modn)
val -= 2*modn;
else if (val >= modn)
val -= modn;
return val;
}
if (len > 50)
return _nmod_poly_evaluate_nmod_rectangular(poly, len, c, mod);

// use n_mulmod_shoup, non-lazy variant
else
{
const ulong c_precomp = n_mulmod_precomp_shoup(c, modn);
return _nmod_poly_evaluate_nmod_precomp(poly->coeffs, poly->length, c, c_precomp, modn);
}
const ulong modn = mod.n;

if (modn <= LAZY_MAX)
{
const ulong c_precomp = n_mulmod_precomp_shoup(c, modn);
ulong val = _nmod_poly_evaluate_nmod_precomp_lazy(poly, len, c, c_precomp, modn);
/* Correct excess. */
if (val >= 2*modn)
val -= 2*modn;
else if (val >= modn)
val -= modn;
return val;
}
else
{
const ulong c_precomp = n_mulmod_precomp_shoup(c, modn);
return _nmod_poly_evaluate_nmod_precomp(poly, len, c, c_precomp, modn);
}
}

ulong
nmod_poly_evaluate_nmod(const nmod_poly_t poly, ulong c)
{
return _nmod_poly_evaluate_nmod(poly->coeffs, poly->length, c, poly->mod);
}

Loading
Loading