#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)
{
return _gr_poly_evaluate_horner(y, poly, len, x, ctx);
}
else
{
slong i, j, k, coeff_index, l, m;
gr_ptr xs, ys, tmp, partial_results;
k = n_sqrt(len)+1;
j = (len + k - 1) / k;
TMP_INIT;
TMP_START;
tmp = TMP_ALLOC(sz * k);
GR_TMP_INIT_VEC(xs, j + 1, ctx);
GR_TMP_INIT_VEC(ys, k, ctx);
GR_TMP_INIT_VEC(partial_results, j, ctx);
status |= _gr_vec_set_powers(xs, x, j + 1, ctx);
status |= _gr_vec_set_powers(ys, GR_ENTRY(xs, j, sz), k, ctx);
for (l = 0; l < j; l++)
{
i = 0;
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_VEC(xs, j + 1, 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);
}