v

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