1"""Evaluate polynomials.
2
3All of the coefficients are stored in reverse order, i.e. if the
4polynomial is
5
6    u_n x^n + u_{n - 1} x^{n - 1} + ... + u_0,
7
8then coeffs[0] = u_n, coeffs[1] = u_{n - 1}, ..., coeffs[n] = u_0.
9
10References
11----------
12[1] Knuth, "The Art of Computer Programming, Volume II"
13
14"""
15from ._complexstuff cimport zabs
16
17cdef extern from "_c99compat.h":
18    double sc_fma(double x, double y, double z) nogil
19
20
21cdef inline double complex cevalpoly(double *coeffs, int degree,
22                                     double complex z) nogil:
23    """Evaluate a polynomial with real coefficients at a complex point.
24
25    Uses equation (3) in section 4.6.4 of [1]. Note that it is more
26    efficient than Horner's method.
27
28    """
29    cdef:
30        int j
31        double a = coeffs[0]
32        double b = coeffs[1]
33        double r = 2*z.real
34        double s = z.real*z.real + z.imag*z.imag
35        double tmp
36
37    for j in range(2, degree + 1):
38        tmp = b
39        b = sc_fma(-s, a, coeffs[j])
40        a = sc_fma(r, a, tmp)
41    return z*a + b
42