rsl-polynomials 0.1.5

A re-write of GSL's Polynomials Routines in 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
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94

use is_close::is_close;
use num::complex::Complex64;

use crate::{PolyError, Polynomial};

// GSL's tests use this tolerance
const EPS: f64 = 100.0 * f64::EPSILON;

#[test]
fn test_build_polynomial() {
    assert!(Polynomial::build(&Vec::<f64>::new()).is_ok());
    assert!(Polynomial::build(&[1.0]).is_ok());
    assert!(Polynomial::build(&[1.0, 2.0, 3.0]).is_ok());
}

#[test]
fn test_new_polynomial() {
    let float_poly = Polynomial::<f64>::new();
    let complex_poly = Polynomial::<Complex64>::new();

    assert_eq!(float_poly.coef, [0.0]);
    assert_eq!(complex_poly.coef, [Complex64::new(0.0, 0.0)]);
}

#[test]
fn test_default_polynomial() {
    let default_float_poly = Polynomial::<f64>::default();
    let default_complex_poly = Polynomial::<Complex64>::default();

    assert_eq!(default_float_poly.coef, [0.0]);
    assert_eq!(default_complex_poly.coef, [Complex64::new(0.0, 0.0)]);
}

#[test]
fn test_build_polynomial_invalid() {
    let poly1 = Polynomial::build(&[1.0, 2.0, std::f64::NAN]);
    let poly2 = Polynomial::build(&[1.0, 2.0, std::f64::INFINITY]);

    assert!(matches!(poly1.unwrap_err(), PolyError::InvalidCoefficients));
    assert!(matches!(poly2.unwrap_err(), PolyError::InvalidCoefficients));
}

#[test]
#[rustfmt::skip]
fn test_trim_trailing_zeros() {
    let poly0 = Polynomial::build(&[0.0]).unwrap().to_trimmed();
    let poly1 = Polynomial::build(&[0.0, 1.0, 2.0]).unwrap().to_trimmed();
    let poly2 = Polynomial::build(&[0.0, 1.0, 2.0, 0.0, 0.0]).unwrap().to_trimmed();
    let poly3 = Polynomial::build(&[1.0, 2.0]).unwrap().to_trimmed();
    let poly4 = Polynomial::build(&[1.0, 2.0, 0.0, 0.0]).unwrap().to_trimmed();
    let poly5 = Polynomial::build(&[1.0, 0.0, 2.0]).unwrap().to_trimmed();

    assert_eq!(poly0.coef, [0.0]);
    assert_eq!(poly1.coef, [0.0, 1.0, 2.0]);
    assert_eq!(poly2.coef, [0.0, 1.0, 2.0]);
    assert_eq!(poly3.coef, [1.0, 2.0]);
    assert_eq!(poly4.coef, [1.0, 2.0]);
    assert_eq!(poly5.coef, [1.0, 0.0, 2.0]);
}

#[test]
fn test_debug() {
    let poly = Polynomial::build(&[0.0]).unwrap();
    let _ = format!("{:?}", poly);
    let _ = format!("{:#?}", poly);
}

#[test]
#[rustfmt::skip]
fn test_monic() {
    let poly0 = Polynomial::build(&[0.0]).unwrap().to_monic();
    let poly1 = Polynomial::build(&[0.0, 8.0, 4.0]).unwrap().to_monic();
    let poly2 = Polynomial::build(&[6.0, 0.0, 3.0]).unwrap().to_monic();
    let poly3 = Polynomial::build(&[0.0, 2.0, 2.0, 0.0, 0.0]).unwrap().to_monic();

    assert_eq!(poly0.coef, [0.0]);
    assert_eq!(poly1.coef, [0.0, 2.0, 1.0]);
    assert_eq!(poly2.coef, [2.0, 0.0, 1.0]);
    assert_eq!(poly3.coef, [0.0, 1.0, 1.0]);
}

#[test]
fn test_to_depressed_cubic() {
    // Example: https://www.johndcook.com/blog/2022/11/19/how-to-depress-a-cubic/
    let poly1 = Polynomial::build(&[22.0, 20.0, 19.0, 11.0])
        .unwrap()
        .to_depressed_cubic()
        .unwrap();

    assert!(is_close!(poly1.coef[0], 47972.0 / 35937.0, rel_tol = EPS));
    assert!(is_close!(poly1.coef[1], 299.0 / 363.0, rel_tol = EPS));
    assert!(is_close!(poly1.coef[2], 0.0, rel_tol = EPS));
    assert!(is_close!(poly1.coef[3], 1.0, rel_tol = EPS));
}