diff --git a/doc/source/nmod_poly.rst b/doc/source/nmod_poly.rst
index a5aed2c419..d9798379bf 100644
--- a/doc/source/nmod_poly.rst
+++ b/doc/source/nmod_poly.rst
@@ -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)
 
diff --git a/src/nmod_poly.h b/src/nmod_poly.h
index 21f0d6e653..49e8e6d261 100644
--- a/src/nmod_poly.h
+++ b/src/nmod_poly.h
@@ -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);
diff --git a/src/nmod_poly/evaluate_nmod.c b/src/nmod_poly/evaluate_nmod.c
index e357c5b844..e6c65afee0 100644
--- a/src/nmod_poly/evaluate_nmod.c
+++ b/src/nmod_poly/evaluate_nmod.c
@@ -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.
 
@@ -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;
@@ -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);
+}
+
diff --git a/src/nmod_poly/profile/p-evaluate_nmod.c b/src/nmod_poly/profile/p-evaluate_nmod.c
index b9854c32c0..6a01b3b356 100644
--- a/src/nmod_poly/profile/p-evaluate_nmod.c
+++ b/src/nmod_poly/profile/p-evaluate_nmod.c
@@ -57,7 +57,7 @@ void sample_interface(void * arg, ulong count)
     FLINT_TEST_CLEAR(state);
 }
 
-void sample_generic(void * arg, ulong count)
+void sample_horner(void * arg, ulong count)
 {
     ulong n;
     nmod_t mod;
@@ -87,7 +87,7 @@ void sample_generic(void * arg, ulong count)
 
         prof_start();
         for (j = 0; j < 100; j++)
-            _nmod_poly_evaluate_nmod(poly, length, pt, mod);
+            _nmod_poly_evaluate_nmod_horner(poly, length, pt, mod);
         prof_stop();
     }
 
@@ -184,76 +184,116 @@ void sample_precomp_lazy(void * arg, ulong count)
     FLINT_TEST_CLEAR(state);
 }
 
+void sample_rectangular(void * arg, ulong count)
+{
+    ulong n;
+    nmod_t mod;
+    info_t * info = (info_t *) arg;
+    flint_bitcnt_t bits = info->bits;
+    slong length = info->length;
+    slong i, j;
+
+    ulong * poly;
+    ulong pt;
+
+    FLINT_TEST_INIT(state);
+
+    poly = _nmod_vec_init(length);
+
+    for (i = 0; i < count; i++)
+    {
+        n = n_randbits(state, bits);
+        if (n == UWORD(0)) n++;
+
+        nmod_init(&mod, n);
+
+        for (j = 0; j < length; j++)
+            poly[j] = n_randint(state, n);
+
+        pt = n_randint(state, n);
+
+        prof_start();
+        for (j = 0; j < 100; j++)
+        {
+            _nmod_poly_evaluate_nmod_rectangular(poly, length, pt, mod);
+        }
+        prof_stop();
+    }
+
+    _nmod_vec_clear(poly);
+
+    FLINT_TEST_CLEAR(state);
+}
+
+#define NUM_LENGTHS 19
+#define NUM_BITS 5
+
 int main(void)
 {
-    slong lengths[18] = {1, 2, 3, 4, 6, 8,
-                         10, 12, 16, 20, 32, 45,
+    slong lengths[NUM_LENGTHS] = {1, 2, 3, 4, 6, 8,
+                         10, 12, 16, 20, 32, 45, 56,
                          64, 128, 256, 1024, 8192, 65536};
+    int bits[NUM_BITS] = { 16, 60, 62, 63, 64 };
+    int ibits;
 
     double min, max;
-    double mins[18]; // note: max seems to be consistently identical or extremely close to min
-    double mins_generic[18];
-    double mins_precomp[18];
-    double mins_precomp_lazy[18];
+    double mins[NUM_LENGTHS]; // note: max seems to be consistently identical or extremely close to min
+    double mins_horner[NUM_LENGTHS];
+    double mins_precomp[NUM_LENGTHS];
+    double mins_precomp_lazy[NUM_LENGTHS];
+    double mins_rectangular[NUM_LENGTHS];
     info_t info;
     flint_bitcnt_t i;
+    
 
     flint_printf("unit: all measurements in c/l\n");
-    flint_printf("profiled: interface | generic | precomp | precomp_lazy\n");
+    flint_printf("profiled: interface | horner | precomp | precomp_lazy | rectangular\n");
 
-    for (i = 62; i <= FLINT_BITS; i++)
+    for (ibits = 0; ibits < NUM_BITS; ibits++)
     {
+        i = bits[ibits];
         info.bits = i;
 
         printf("nbits = %ld\n", i);
-        for (int len = 0; len < 18; ++len)
+        for (int len = 0; len < NUM_LENGTHS; ++len)
         {
             info.length = lengths[len];
 
             prof_repeat(&min, &max, sample_interface, (void *) &info);
-            mins[len-1] = min;
-            prof_repeat(&min, &max, sample_generic, (void *) &info);
-            mins_generic[len-1] = min;
+            mins[len] = min;
+            prof_repeat(&min, &max, sample_horner, (void *) &info);
+            mins_horner[len] = min;
             prof_repeat(&min, &max, sample_precomp, (void *) &info);
-            mins_precomp[len-1] = min;
+            mins_precomp[len] = min;
             prof_repeat(&min, &max, sample_precomp_lazy, (void *) &info);
-            mins_precomp_lazy[len-1] = min;
+            mins_precomp_lazy[len] = min;
+            prof_repeat(&min, &max, sample_rectangular, (void *) &info);
+            mins_rectangular[len] = min;
         }
 
-        if (i < FLINT_BITS-1)
-        {
-            for (int len = 0; len < 18; ++len)
-            {
-                flint_printf("   len %ld\t%.2lf\t%.2lf\t%.2lf\t%.2lf",
-                        lengths[len],
-                        (mins[len-1]/(double)FLINT_CLOCK_SCALE_FACTOR)/(lengths[len]*100),
-                        (mins_generic[len-1]/(double)FLINT_CLOCK_SCALE_FACTOR)/(lengths[len]*100),
-                        (mins_precomp[len-1]/(double)FLINT_CLOCK_SCALE_FACTOR)/(lengths[len]*100),
-                        (mins_precomp_lazy[len-1]/(double)FLINT_CLOCK_SCALE_FACTOR)/(lengths[len]*100));
-                flint_printf("\n");
-            }
-        }
-        else if (i < FLINT_BITS)
+        if (i < FLINT_BITS)
         {
-            for (int len = 0; len < 18; ++len)
+            for (int len = 0; len < NUM_LENGTHS; ++len)
             {
-                flint_printf("   len %ld\t%.2lf\t%.2lf\t%.2lf\t na",
+                flint_printf("   len %ld\t%.2lf\t%.2lf\t%.2lf\t%.2lf\t%.2lf",
                         lengths[len],
-                        (mins[len-1]/(double)FLINT_CLOCK_SCALE_FACTOR)/(lengths[len]*100),
-                        (mins_generic[len-1]/(double)FLINT_CLOCK_SCALE_FACTOR)/(lengths[len]*100),
-                        (mins_precomp[len-1]/(double)FLINT_CLOCK_SCALE_FACTOR)/(lengths[len]*100),
-                        (mins_precomp_lazy[len-1]/(double)FLINT_CLOCK_SCALE_FACTOR)/(lengths[len]*100));
+                        (mins[len]/(double)FLINT_CLOCK_SCALE_FACTOR)/(lengths[len]*100),
+                        (mins_horner[len]/(double)FLINT_CLOCK_SCALE_FACTOR)/(lengths[len]*100),
+                        (mins_precomp[len]/(double)FLINT_CLOCK_SCALE_FACTOR)/(lengths[len]*100),
+                        (mins_precomp_lazy[len]/(double)FLINT_CLOCK_SCALE_FACTOR)/(lengths[len]*100),
+                        (mins_rectangular[len]/(double)FLINT_CLOCK_SCALE_FACTOR)/(lengths[len]*100));
                 flint_printf("\n");
             }
         }
-        else  // i == FLINT_BITS
+        else
         {
-            for (int len = 0; len < 18; ++len)
+            for (int len = 0; len < NUM_LENGTHS; ++len)
             {
-                flint_printf("   len %ld\t%.2lf\t%.2lf\t na \t na",
+                flint_printf("   len %ld\t%.2lf\t%.2lf\t na \t na\t%.2lf",
                         lengths[len],
-                        (mins[len-1]/(double)FLINT_CLOCK_SCALE_FACTOR)/(lengths[len]*100),
-                        (mins_generic[len-1]/(double)FLINT_CLOCK_SCALE_FACTOR)/(lengths[len]*100));
+                        (mins[len]/(double)FLINT_CLOCK_SCALE_FACTOR)/(lengths[len]*100),
+                        (mins_horner[len]/(double)FLINT_CLOCK_SCALE_FACTOR)/(lengths[len]*100),
+                        (mins_rectangular[len]/(double)FLINT_CLOCK_SCALE_FACTOR)/(lengths[len]*100));
                 flint_printf("\n");
             }
         }
diff --git a/src/nmod_poly/test/t-evaluate_nmod.c b/src/nmod_poly/test/t-evaluate_nmod.c
index 95c3196ed2..153311d2c6 100644
--- a/src/nmod_poly/test/t-evaluate_nmod.c
+++ b/src/nmod_poly/test/t-evaluate_nmod.c
@@ -104,7 +104,7 @@ TEST_FUNCTION_START(nmod_poly_evaluate_nmod, state)
 
         // variants
         nmod_poly_set(a, acopy);
-        eval1_v1 = _nmod_poly_evaluate_nmod(a->coeffs, a->length, c, a->mod);
+        eval1_v1 = _nmod_poly_evaluate_nmod_horner(a->coeffs, a->length, c, a->mod);
         eval1_v1 = n_addmod(eval1_v1, _nmod_poly_evaluate_nmod(b->coeffs, b->length, c, b->mod), n);
 
         if (a->mod.norm > 0)
@@ -145,7 +145,7 @@ TEST_FUNCTION_START(nmod_poly_evaluate_nmod, state)
         }
 
         nmod_poly_add(a, a, b);
-        eval2_v1 = _nmod_poly_evaluate_nmod(a->coeffs, a->length, c, a->mod);
+        eval2_v1 = _nmod_poly_evaluate_nmod_horner(a->coeffs, a->length, c, a->mod);
         if (a->mod.norm > 0)
         {
             ulong c_precomp = n_mulmod_precomp_shoup(c, n);