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

//! Provides functions that don't have a numerical solution and must
//! be solved computationally (e.g. evaluation of a polynomial)

/// evaluates a polynomial at `z` where `coeff` are the coeffecients
/// to a polynomial of order `k` where `k` is the length of `coeff` and the
/// coeffecient
/// to the `k`th power is the `k`th element in coeff. E.g. [3,-1,2] equates to
/// `2z^2 - z + 3`
///
/// # Remarks
///
/// Returns 0 for a 0 length coefficient slice
pub fn polynomial(z: f64, coeff: &[f64]) -> f64 {
    let n = coeff.len();
    if n == 0 {
        return 0.0;
    }

    unsafe {
        let mut sum = *coeff.get_unchecked(n - 1);
        for i in (0..n - 1).rev() {
            sum *= z;
            sum += *coeff.get_unchecked(i);
        }
        sum
    }
}

#[cfg_attr(rustfmt, rustfmt_skip)]
#[cfg(test)]
mod test {
    use std::f64;

    // these tests probably could be more robust
    #[test]
    fn test_polynomial() {
        let empty: [f64; 0] = [];
        assert_eq!(super::polynomial(2.0, &empty), 0.0);

        let zero = [0.0];
        assert_eq!(super::polynomial(2.0, &zero), 0.0);

        let mut coeff = [1.0, 0.0, 5.0];
        assert_eq!(super::polynomial(2.0, &coeff), 21.0);

        coeff = [-5.0, -2.0, 3.0];
        assert_eq!(super::polynomial(2.0, &coeff), 3.0);
        assert_eq!(super::polynomial(-2.0, &coeff), 11.0);

        let large_coeff = [-1.35e3, 2.5e2, 8.0, -4.0, 1e2, 3.0];
        assert_eq!(super::polynomial(5.0, &large_coeff), 71475.0);
        assert_eq!(super::polynomial(-5.0, &large_coeff), 51225.0);

        coeff = [f64::INFINITY, -2.0, 3.0];
        assert_eq!(super::polynomial(2.0, &coeff), f64::INFINITY);
        assert_eq!(super::polynomial(-2.0, &coeff), f64::INFINITY);

        coeff = [f64::NEG_INFINITY, -2.0, 3.0];
        assert_eq!(super::polynomial(2.0, &coeff), f64::NEG_INFINITY);
        assert_eq!(super::polynomial(-2.0, &coeff), f64::NEG_INFINITY);

        coeff = [f64::NAN, -2.0, 3.0];
        assert!(super::polynomial(2.0, &coeff).is_nan());
        assert!(super::polynomial(-2.0, &coeff).is_nan());
    }
}