roots 0.0.8

Library of well known algorithms for numerical root finding.
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
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123

// Copyright (c) 2015, Mikhail Vorotilov
// All rights reserved.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are met:
//
// * Redistributions of source code must retain the above copyright notice, this
//   list of conditions and the following disclaimer.
//
// * Redistributions in binary form must reproduce the above copyright notice,
//   this list of conditions and the following disclaimer in the documentation
//   and/or other materials provided with the distribution.
//
// THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
// AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
// DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
// FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
// DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
// SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
// CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
// OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
// OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

use super::super::FloatType;
use super::super::Roots;

/// Solves a depressed quartic equation x^4 + a2*x^2 + a1*x + a0 = 0.
///
/// Returned roots are ordered. Precision is about 1e-14 for f64.
///
/// # Examples
///
/// ```
/// use roots::find_roots_quartic_depressed;
///
/// let one_root = find_roots_quartic_depressed(1f64, 0f64, 0f64);
/// // Returns Roots::One([0f64]) as 'x^4 = 0' has one root 0
///
/// let two_roots = find_roots_quartic_depressed(1f32, 0f32, -1f32);
/// // Returns Roots::Two([-1f32, 1f32]) as 'x^4 - 1 = 0' has roots -1 and 1
/// ```
pub fn find_roots_quartic_depressed<F: FloatType>(a2: F, a1: F, a0: F) -> Roots<F> {
    // Handle non-standard cases
    if a1 == F::zero() {
        // a1 = 0; x^4 + a2*x^2 + a0 = 0; solve biquadratic equation
        super::biquadratic::find_roots_biquadratic(F::one(), a2, a0)
    } else if a0 == F::zero() {
        // a0 = 0; x^4 + a2*x^2 + a1*x = 0; reduce to normalized cubic and add zero root
        super::cubic_normalized::find_roots_cubic_normalized(F::zero(), a2, a1).add_new_root(F::zero())
    } else {
        let _2 = F::from(2i16);
        let _5 = F::from(5i16);

        // Solve the auxiliary equation y^3 + (5/2)*a2*y^2 + (2*a2^2-a0)*y + (a2^3/2 - a2*a0/2 - a1^2/8) = 0
        let a2_pow_2 = a2 * a2;
        let a1_div_2 = a1 / _2;
        let b2 = a2 * _5 / _2;
        let b1 = _2 * a2_pow_2 - a0;
        let b0 = (a2_pow_2 * a2 - a2 * a0 - a1_div_2 * a1_div_2) / _2;

        // At least one root always exists. The last root is the maximal one.
        let resolvent_roots = super::cubic_normalized::find_roots_cubic_normalized(b2, b1, b0);
        let y = resolvent_roots.as_ref().iter().last().unwrap();

        let _a2_plus_2y = a2 + _2 * *y;
        if _a2_plus_2y > F::zero() {
            let sqrt_a2_plus_2y = _a2_plus_2y.sqrt();
            let q0a = a2 + *y - a1_div_2 / sqrt_a2_plus_2y;
            let q0b = a2 + *y + a1_div_2 / sqrt_a2_plus_2y;

            let mut roots = super::quadratic::find_roots_quadratic(F::one(), sqrt_a2_plus_2y, q0a);
            for x in super::quadratic::find_roots_quadratic(F::one(), -sqrt_a2_plus_2y, q0b)
                .as_ref()
                .iter()
            {
                roots = roots.add_new_root(*x);
            }
            roots
        } else {
            Roots::No([])
        }
    }
}

#[cfg(test)]
mod test {
    use super::super::super::*;

    #[test]
    fn test_find_roots_quartic_depressed() {
        assert_eq!(find_roots_quartic_depressed(0f32, 0f32, 0f32), Roots::One([0f32]));
        assert_eq!(find_roots_quartic_depressed(1f32, 1f32, 1f32), Roots::No([]));

        // Thanks WolframAlpha for the test data
        match find_roots_quartic_depressed(1f64, 1f64, -1f64) {
            Roots::Two(x) => {
                assert_float_array_eq!(1e-15, x, [-1f64, 0.5698402909980532659114f64]);
            }
            _ => {
                assert!(false);
            }
        }

        match find_roots_quartic_depressed(-10f64, 5f64, 1f64) {
            Roots::Four(x) => {
                assert_float_array_eq!(
                    1e-15,
                    x,
                    [
                        -3.3754294311910698f64,
                        -0.1531811728532153f64,
                        0.67861075799846644f64,
                        2.84999984604581877f64
                    ]
                );
            }
            _ => {
                assert!(false);
            }
        }
    }
}