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
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
/* polevl.c
* p1evl.c
*
* Evaluate polynomial
*
*
*
* SYNOPSIS:
*
* int N;
* double x, y, coef[N+1], polevl[];
*
* y = polevl( x, coef, N );
*
*
*
* DESCRIPTION:
*
* Evaluates polynomial of degree N:
*
* 2 N
* y = C + C x + C x +...+ C x
* 0 1 2 N
*
* Coefficients are stored in reverse order:
*
* coef[0] = C , ..., coef[N] = C .
* N 0
*
* The function p1evl() assumes that c_N = 1.0 so that coefficent
* is omitted from the array. Its calling arguments are
* otherwise the same as polevl().
*
*
* SPEED:
*
* In the interest of speed, there are no checks for out
* of bounds arithmetic. This routine is used by most of
* the functions in the library. Depending on available
* equipment features, the user may wish to rewrite the
* program in microcode or assembly language.
*
*/
/*
* Cephes Math Library Release 2.1: December, 1988
* Copyright 1984, 1987, 1988 by Stephen L. Moshier
* Direct inquiries to 30 Frost Street, Cambridge, MA 02140
*/
/* Sources:
* [1] Holin et. al., "Polynomial and Rational Function Evaluation",
* https://www.boost.org/doc/libs/1_61_0/libs/math/doc/html/math_toolkit/roots/rational.html
*/
/* Scipy changes:
* - 06-23-2016: add code for evaluating rational functions
*/
// from https://github.com/scipy/scipy/blob/c4ce0c4560bc635867512c4d2ea6db6f666d3eeb/scipy/special/cephes/polevl.h#L67
pub fn polevl(x: f64, coef: &[f64], n: usize) -> f64 {
let mut coef_iter = coef.iter().take(n+1);
let mut ans = *coef_iter.next().unwrap();
for i in coef_iter {
ans = ans * x + *i;
}
ans
}
/* p1evl() */
/* N
* Evaluate polynomial when coefficient of x is 1.0.
* That is, C_{N} is assumed to be 1, and that coefficient
* is not included in the input array coef.
* coef must have length N and contain the polynomial coefficients
* stored as
* coef[0] = C_{N-1}
* coef[1] = C_{N-2}
* ...
* coef[N-2] = C_1
* coef[N-1] = C_0
* Otherwise same as polevl.
*/
pub fn p1evl(x: f64, coef: &[f64], n: usize) -> f64
{
let mut coef_iter = coef.iter().take(n);
let mut ans = x + *coef_iter.next().unwrap();
for i in coef_iter {
ans = ans * x + *i;
}
ans
}
/* Evaluate a rational function. See [1]. */
pub fn ratevl(x: f64, num: &[f64], m: isize, denom: &[f64], n: isize) -> f64 {
let absx = x.abs();
let (dir, mut p, y) = if absx > 1.0 {
/* Evaluate as a polynomial in 1/x. */
(-1_isize, m, 1.0 / x)
} else {
(1_isize, 0_isize, x)
};
/* Evaluate the numerator */
let mut num_ans = num[p as usize];
p += dir;
for _ in 1..=m {
num_ans = num_ans * y + num[p as usize];
p += dir;
}
/* Evaluate the denominator */
if absx > 1.0 {
p = n;
} else {
p = 0;
}
let mut denom_ans = denom[p as usize];
p += dir;
for _ in 1..=n {
denom_ans = denom_ans * y + denom[p as usize];
p += dir;
}
if absx > 1.0 {
let i = n - m;
x.powi(i as i32) * num_ans / denom_ans
} else {
num_ans / denom_ans
}
}
#[cfg(test)]
mod polevl_tests {
use super::*;
const TAYLOR0: [f64; 10] = [
-1.0000000009110164892,
-1.0000000057646759799,
-9.9999983138417361078e-1,
-1.0000013011460139596,
-1.000001940896320456,
-9.9987929950057116496e-1,
-1.000785194477042408,
-1.0031782279542924256,
-9.1893853320467274178e-1,
-1.5,
];
// #[test]
// fn polevl_timing() {
// let mut s: f64 = 0.0;
// let now = std::time::Instant::now();
// for i in 0..100000 {
// let x = i as f64 * 0.01 / 100000.0 - 0.01;
// s += polevl(x, &TAYLOR0, 9);
// }
// println!("Time: {}", now.elapsed().as_micros());
// println!("Sum: {}", s);
// let mut s: f64 = 0.0;
// let now = std::time::Instant::now();
// for i in 0..100000 {
// let x = i as f64 * 0.01 / 100000.0 - 0.01;
// s += polevl2(x, &TAYLOR0, 9);
// }
// println!("Time: {}", now.elapsed().as_micros());
// println!("Sum: {}", s);
// }
#[test]
fn polevl_test() {
let coefs = [2.0, 3.0, 4.0];
assert_eq!(polevl(1.0, &coefs, 1), 5.0);
assert_eq!(polevl(1.0, &coefs, 2), 9.0);
assert_eq!(polevl(2.0, &coefs, 2), 18.0);
}
}