cheby 0.2.0

Unit-safe Chebyshev approximation and spectral numerics for Rust.
Documentation 
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

//! Clenshaw evaluation for Chebyshev series.
//!
//! For coefficients `a_k`, the represented polynomial is
//! `sum_k a_k T_k(x)`. Clenshaw evaluates that basis recurrence directly:
//! it walks the coefficients backward with
//! `b_k = 2x b_{k+1} - b_{k+2} + a_k`, then returns
//! `a_0 + x b_1 - b_2`. This avoids converting the series to powers of `x`.

use super::ChebyScalar;

/// Evaluate a Chebyshev series at normalized coordinate `x`.
#[inline]
pub fn evaluate<T: ChebyScalar>(coeffs: &[T], x: f64) -> T {
    let n = coeffs.len();
    if n == 0 {
        return T::zero();
    }
    if n == 1 {
        return coeffs[0];
    }

    let two_x = 2.0 * x;
    let mut b1 = T::zero();
    let mut b2 = T::zero();
    for k in (1..n).rev() {
        let b0 = b1 * two_x - b2 + coeffs[k];
        b2 = b1;
        b1 = b0;
    }
    coeffs[0] + b1 * x - b2
}

/// Evaluate `df/dx` for a normalized-coordinate Chebyshev series.
#[inline]
pub fn evaluate_derivative<T: ChebyScalar>(coeffs: &[T], x: f64) -> T {
    evaluate_both(coeffs, x).1
}

/// Evaluate both value and derivative with respect to normalized `x`.
#[inline]
pub fn evaluate_both<T: ChebyScalar>(coeffs: &[T], x: f64) -> (T, T) {
    let n = coeffs.len();
    if n == 0 {
        return (T::zero(), T::zero());
    }
    if n == 1 {
        return (coeffs[0], T::zero());
    }

    let two_x = 2.0 * x;
    let mut b1 = T::zero();
    let mut b2 = T::zero();
    let mut db1 = T::zero();
    let mut db2 = T::zero();

    for k in (1..n).rev() {
        let b0 = b1 * two_x - b2 + coeffs[k];
        let db0 = db1 * two_x - db2 + b1 * 2.0;
        b2 = b1;
        b1 = b0;
        db2 = db1;
        db1 = db0;
    }

    let value = coeffs[0] + b1 * x - b2;
    let derivative = b1 + db1 * x - db2;
    (value, derivative)
}