1 | module math |
2 | |
3 | import math.internal |
4 | |
5 | fn poly_n_eval(c []f64, n int, x f64) f64 { |
6 | if c.len == 0 { |
7 | panic('coeficients can not be empty') |
8 | } |
9 | len := int(min(c.len, n)) |
10 | mut ans := c[len - 1] |
11 | for e in c[..len - 1] { |
12 | ans = e + x * ans |
13 | } |
14 | return ans |
15 | } |
16 | |
17 | fn poly_n_1_eval(c []f64, n int, x f64) f64 { |
18 | if c.len == 0 { |
19 | panic('coeficients can not be empty') |
20 | } |
21 | len := int(min(c.len, n)) - 1 |
22 | mut ans := c[len - 1] |
23 | for e in c[..len - 1] { |
24 | ans = e + x * ans |
25 | } |
26 | return ans |
27 | } |
28 | |
29 | [inline] |
30 | fn poly_eval(c []f64, x f64) f64 { |
31 | return poly_n_eval(c, c.len, x) |
32 | } |
33 | |
34 | [inline] |
35 | fn poly_1_eval(c []f64, x f64) f64 { |
36 | return poly_n_1_eval(c, c.len, x) |
37 | } |
38 | |
39 | // data for a Chebyshev series over a given interval |
40 | struct ChebSeries { |
41 | pub: |
42 | c []f64 // coefficients |
43 | order int // order of expansion |
44 | a f64 // lower interval point |
45 | b f64 // upper interval point |
46 | } |
47 | |
48 | fn (cs ChebSeries) eval_e(x f64) (f64, f64) { |
49 | mut d := 0.0 |
50 | mut dd := 0.0 |
51 | y := (2.0 * x - cs.a - cs.b) / (cs.b - cs.a) |
52 | y2 := 2.0 * y |
53 | mut e_ := 0.0 |
54 | mut temp := 0.0 |
55 | for j := cs.order; j >= 1; j-- { |
56 | temp = d |
57 | d = y2 * d - dd + cs.c[j] |
58 | e_ += abs(y2 * temp) + abs(dd) + abs(cs.c[j]) |
59 | dd = temp |
60 | } |
61 | temp = d |
62 | d = y * d - dd + 0.5 * cs.c[0] |
63 | e_ += abs(y * temp) + abs(dd) + 0.5 * abs(cs.c[0]) |
64 | return d, f64(internal.f64_epsilon) * e_ + abs(cs.c[cs.order]) |
65 | } |