| 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 | } |