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