Evaluation

fmpz_poly_evaluate_fmpz computes f(a) for an integer a, fmpz_poly_evaluate_fmpq for a rational a (in lowest terms), and fmpz_poly_evaluate_mod the value modulo a word-size integer. Values at many points are computed one at a time by fmpz_poly_evaluate_fmpz_vec. For numerical evaluation with error bounds at real or complex points, use arb_fmpz_poly_evaluate_arb and arb_fmpz_poly_evaluate_acb.

Horner's rule and divide and conquer#

Horner's rule f(a) = (\cdots(f_n a + f_{n-1}) a + \cdots) a + f_0 uses n multiplications, but the intermediate results grow: for a of k bits, the i-th product multiplies a number of about ik bits by one of k bits, so the cost is O(n^2 k^2) with classical integer multiplication and still quadratic in n with fast multiplication, since one operand is always small. The divide-and-conquer algorithm writes f = f_0 + x^m f_1 with m = n/2, evaluates the halves recursively and returns f_0(a) + a^m f_1(a), with the powers a^{2^j} computed once. Its products are balanced, so with fast integer multiplication the total cost is O(M(nk) \log n) for the cost M of multiplying integers of nk bits. fmpz_poly_evaluate_fmpz switches from Horner's rule at length 50.

At a rational point a = p/q, both algorithms compute the integer q^n f(p/q) = \sum_i f_i p^i q^{n-i} and divide out the common factor at the end. fmpz_poly_evaluate_fmpq uses Horner's rule for length less than 40 or when the denominator has more than 0.003 n^2 bits, and divide and conquer otherwise.

Evaluation modulo n reduces each coefficient once and runs Horner's rule in word arithmetic. To evaluate a polynomial modulo n at many points, convert it with fmpz_poly_get_nmod_poly and use nmod_poly_evaluate_nmod_vec, which is quasi-linear.

Examples#

#include <flint/flint.h>
#include <flint/fmpz.h>
#include <flint/fmpq.h>
#include <flint/fmpz_poly.h>

int main(void)
{
    fmpz_poly_t f;
    fmpz_t a, y;
    fmpq_t b, z;

    fmpz_poly_init(f);
    fmpz_init(a);
    fmpz_init(y);
    fmpq_init(b);
    fmpq_init(z);

    fmpz_poly_set_str(f, "4  -1 -1 0 1");     /* x^3 - x - 1 */

    fmpz_set_si(a, -3);
    fmpz_poly_evaluate_fmpz(y, f, a);
    flint_printf("f(-3) = "); fmpz_print(y); flint_printf("\n");

    fmpq_set_si(b, 4, 3);
    fmpz_poly_evaluate_fmpq(z, f, b);
    flint_printf("f(4/3) = "); fmpq_print(z); flint_printf("\n");

    flint_printf("f(10^6) mod 1000003 = %wu\n",
                 fmpz_poly_evaluate_mod(f, 1000000, 1000003));

    fmpz_poly_clear(f);
    fmpz_clear(a);
    fmpz_clear(y);
    fmpq_clear(b);
    fmpq_clear(z);
    return 0;
}

f(-3) = -25
f(4/3) = 1/27
f(10^6) mod 1000003 = 999978

The divide-and-conquer algorithm and Horner's rule give the same value. The value of a polynomial of degree 999 with 10-bit coefficients at 3^{100} has about 160000 bits:

#include <flint/flint.h>
#include <flint/fmpz.h>
#include <flint/fmpz_poly.h>

int main(void)
{
    flint_rand_t state;
    fmpz_poly_t f;
    fmpz_t a, y1, y2;

    flint_rand_init(state);
    fmpz_poly_init(f);
    fmpz_init(a);
    fmpz_init(y1);
    fmpz_init(y2);

    fmpz_poly_randtest(f, state, 1000, 10);
    fmpz_poly_set_coeff_ui(f, 999, 1);
    fmpz_ui_pow_ui(a, 3, 100);

    fmpz_poly_evaluate_horner_fmpz(y1, f, a);
    fmpz_poly_evaluate_divconquer_fmpz(y2, f, a);
    flint_printf("equal: %d, bits: %wd\n", fmpz_equal(y1, y2), fmpz_bits(y1));

    fmpz_poly_clear(f);
    fmpz_clear(a);
    fmpz_clear(y1);
    fmpz_clear(y2);
    flint_rand_clear(state);
    return 0;
}

equal: 1, bits: 158338