Merge pull request #1542 from David-Berghaus/Add-modular-splitting

Add modular splitting

doc/source/gr_poly.rst

@@ -355,6 +355,9 @@ Evaluation
 .. function:: int _gr_poly_evaluate_rectangular(gr_ptr res, gr_srcptr poly, slong len, gr_srcptr x, gr_ctx_t ctx)
               int gr_poly_evaluate_rectangular(gr_ptr res, const gr_poly_t poly, gr_srcptr x, gr_ctx_t ctx)
 
+.. function:: int _gr_poly_evaluate_modular(gr_ptr res, gr_srcptr poly, slong len, gr_srcptr x, gr_ctx_t ctx)
+              int gr_poly_evaluate_modular(gr_ptr res, const gr_poly_t poly, gr_srcptr x, gr_ctx_t ctx)
+
 .. function:: int _gr_poly_evaluate_horner(gr_ptr res, gr_srcptr poly, slong len, gr_srcptr x, gr_ctx_t ctx)
               int gr_poly_evaluate_horner(gr_ptr res, const gr_poly_t poly, gr_srcptr x, gr_ctx_t ctx)
 

src/gr_poly.h

@@ -269,6 +269,9 @@ WARN_UNUSED_RESULT int gr_poly_rsqrt_series(gr_poly_t res, const gr_poly_t f, sl
 WARN_UNUSED_RESULT int _gr_poly_evaluate_rectangular(gr_ptr res, gr_srcptr poly, slong len, gr_srcptr x, gr_ctx_t ctx);
 WARN_UNUSED_RESULT int gr_poly_evaluate_rectangular(gr_ptr res, const gr_poly_t poly, gr_srcptr x, gr_ctx_t ctx);
 
+WARN_UNUSED_RESULT int _gr_poly_evaluate_modular(gr_ptr res, gr_srcptr poly, slong len, gr_srcptr x, gr_ctx_t ctx);
+WARN_UNUSED_RESULT int gr_poly_evaluate_modular(gr_ptr res, const gr_poly_t poly, gr_srcptr x, gr_ctx_t ctx);
+
 WARN_UNUSED_RESULT int _gr_poly_evaluate_horner(gr_ptr res, gr_srcptr poly, slong len, gr_srcptr x, gr_ctx_t ctx);
 WARN_UNUSED_RESULT int gr_poly_evaluate_horner(gr_ptr res, const gr_poly_t poly, gr_srcptr x, gr_ctx_t ctx);
 

src/gr_poly/evaluate_modular.c

@@ -0,0 +1,88 @@
+/*
+    Copyright (C) 2023 David Berghaus
+
+    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 2.1 of the License, or
+    (at your option) any later version.  See <http://www.gnu.org/licenses/>.
+*/
+
+#include "ulong_extras.h"
+#include "gr_vec.h"
+#include "gr_poly.h"
+
+int
+_gr_poly_evaluate_modular(gr_ptr y, gr_srcptr poly,
+    slong len, gr_srcptr x, gr_ctx_t ctx)
+{
+    slong sz = ctx->sizeof_elem;
+    int status = GR_SUCCESS;
+    gr_method_void_unary_op set_shallow = GR_VOID_UNARY_OP(ctx, SET_SHALLOW);
+
+    if (len <= 2)
+    {
+        if (len == 0)
+            return gr_zero(y, ctx);
+
+        if (len == 1)
+            return gr_set(y, poly, ctx);
+
+        status |= gr_mul(y, x, GR_ENTRY(poly, 1, sz), ctx);
+        status |= gr_add(y, y, GR_ENTRY(poly, 0, sz), ctx);
+    }
+    else
+    {
+        slong i, j, k, coeff_index, l, m;
+        gr_ptr xs, ys, tmp, partial_results;
+        gr_ptr tmp_gr;
+        k = n_sqrt(len)+1;
+        j = (len + k - 1) / k;
+
+        TMP_INIT;
+
+        TMP_START;
+        tmp = TMP_ALLOC(sz * k);
+        GR_TMP_INIT(tmp_gr, ctx);
+        GR_TMP_INIT_VEC(xs, j, ctx);
+        GR_TMP_INIT_VEC(ys, k, ctx);
+        GR_TMP_INIT_VEC(partial_results, j, ctx);
+
+        status |= _gr_vec_set_powers(xs, x, j, ctx);
+        status |= gr_mul(tmp_gr, x, GR_ENTRY(xs, j-1, sz), ctx); /* This is x^j */
+        status |= _gr_vec_set_powers(ys, tmp_gr, k, ctx);
+
+        for (l = 0; l < j; l++){
+            i = 0; /* Count number of coeffs in this row */
+            for (m = 0; m < k; m++)
+            {
+                coeff_index = j*m+l;
+                if (coeff_index < len)
+                { 
+                    set_shallow(GR_ENTRY(tmp, m, sz), GR_ENTRY(poly, coeff_index, sz), ctx);
+                    i++;
+                }
+                else
+                {
+                    break;
+                }
+            }
+            status |= _gr_vec_dot(GR_ENTRY(partial_results, l, sz), NULL, 0, tmp, ys, i, ctx);
+        }
+        status |=_gr_vec_dot(y, NULL, 0, partial_results, xs, j, ctx);
+        GR_TMP_CLEAR(tmp_gr, ctx);
+        GR_TMP_CLEAR_VEC(xs, j, ctx);
+        GR_TMP_CLEAR_VEC(ys, k, ctx);
+        GR_TMP_CLEAR_VEC(partial_results, j, ctx);
+        TMP_END;
+    }
+    return status;
+}
+
+
+int
+gr_poly_evaluate_modular(gr_ptr res, const gr_poly_t f, gr_srcptr x, gr_ctx_t ctx)
+{
+    return _gr_poly_evaluate_modular(res, f->coeffs, f->length, x, ctx);
+}

src/gr_poly/test/t-evaluate_modular.c

@@ -0,0 +1,88 @@
+/*
+    Copyright (C) 2023 Fredrik Johansson
+
+    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 2.1 of the License, or
+    (at your option) any later version.  See <http://www.gnu.org/licenses/>.
+*/
+
+#include "ulong_extras.h"
+#include "gr_poly.h"
+
+int main()
+{
+    slong iter;
+    flint_rand_t state;
+
+    flint_printf("evaluate_modular....");
+    fflush(stdout);
+
+    flint_randinit(state);
+
+    for (iter = 0; iter < 1000; iter++)
+    {
+        int status;
+        gr_ctx_t ctx;
+        gr_poly_t F, G, FG;
+        gr_ptr x, Fx, Gx, FxGx, FGx;
+
+        /* Test F(x) + G(x) = (F + G)(x) */
+        gr_ctx_init_random(ctx, state);
+
+        gr_poly_init(F, ctx);
+        gr_poly_init(G, ctx);
+        gr_poly_init(FG, ctx);
+
+        x = gr_heap_init(ctx);
+        Fx = gr_heap_init(ctx);
+        Gx = gr_heap_init(ctx);
+        FxGx = gr_heap_init(ctx);
+        FGx = gr_heap_init(ctx);
+
+        status = GR_SUCCESS;
+
+        status |= gr_poly_randtest(F, state, 1 + n_randint(state, 10), ctx);
+        status |= gr_poly_randtest(G, state, 1 + n_randint(state, 10), ctx);
+        status |= gr_randtest(x, state, ctx);
+        status |= gr_randtest(Fx, state, ctx);
+
+        status |= gr_poly_add(FG, F, G, ctx);
+
+        status |= gr_poly_evaluate_modular(Fx, F, x, ctx);
+        status |= gr_set(Gx, x, ctx); /* test aliasing */
+        status |= gr_poly_evaluate_modular(Gx, G, Gx, ctx);
+        status |= gr_add(FxGx, Fx, Gx, ctx);
+        status |= gr_poly_evaluate_modular(FGx, FG, x, ctx);
+
+        if (status == GR_SUCCESS && gr_equal(FxGx, FGx, ctx) == T_FALSE)
+        {
+            flint_printf("FAIL\n\n");
+            flint_printf("F = "); gr_poly_print(F, ctx); flint_printf("\n");
+            flint_printf("G = "); gr_poly_print(G, ctx); flint_printf("\n");
+            flint_printf("x = "); gr_print(x, ctx); flint_printf("\n");
+            flint_printf("FxGx = "); gr_print(FxGx, ctx); flint_printf("\n");
+            flint_printf("FGx = "); gr_print(FGx, ctx); flint_printf("\n");
+            flint_abort();
+        }
+
+        gr_poly_clear(F, ctx);
+        gr_poly_clear(G, ctx);
+        gr_poly_clear(FG, ctx);
+
+        gr_heap_clear(x, ctx);
+        gr_heap_clear(Fx, ctx);
+        gr_heap_clear(Gx, ctx);
+        gr_heap_clear(FxGx, ctx);
+        gr_heap_clear(FGx, ctx);
+
+        gr_ctx_clear(ctx);
+    }
+
+    flint_randclear(state);
+    flint_cleanup();
+    flint_printf("PASS\n");
+    return 0;
+}