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