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
// 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 cubic equation x^3 + a1*x + a0 = 0.
///
/// In case more than one roots are present, they are returned in the increasing order.
///
/// # Examples
///
/// ```
/// use roots::find_roots_cubic_depressed;
///
/// let one_root = find_roots_cubic_depressed(0f64, 0f64);
/// // Returns Roots::One([0f64]) as 'x^3 = 0' has one root 0
///
/// let three_roots = find_roots_cubic_depressed(-1f32, 0f32);
/// // Returns Roots::Three([-1f32, -0f32, 1f32]) as 'x^3 - x = 0' has roots -1, 0, and 1
/// ```
pub fn find_roots_cubic_depressed<F: FloatType>(a1: F, a0: F) -> Roots<F> {
let _2 = F::from(2i16);
let _3 = F::from(3i16);
let _4 = F::from(4i16);
let _9 = F::from(9i16);
let _18 = F::from(18i16);
let _27 = F::from(27i16);
let _54 = F::from(54i16);
if a1 == F::zero() {
Roots::One([-a0.cbrt()])
} else if a0 == F::zero() {
super::quadratic::find_roots_quadratic(F::one(), F::zero(), a1).add_new_root(F::zero())
} else {
let d = a0 * a0 / _4 + a1 * a1 * a1 / _27;
if d < F::zero() {
// n*a0^2 + m*a1^3 < 0 => a1 < 0
let a = (-_4 * a1 / _3).sqrt();
let phi = (-_4 * a0 / (a * a * a)).acos() / _3;
Roots::One([a * phi.cos()])
.add_new_root(a * (phi + F::two_third_pi()).cos())
.add_new_root(a * (phi - F::two_third_pi()).cos())
} else {
let sqrt_d = d.sqrt();
let a0_div_2 = a0 / _2;
let x1 = (sqrt_d - a0_div_2).cbrt() - (sqrt_d + a0_div_2).cbrt();
if d == F::zero() {
// one real root and one double root
Roots::One([x1]).add_new_root(a0_div_2)
} else {
// one real root
Roots::One([x1])
}
}
}
}
#[cfg(test)]
mod test {
use super::super::super::*;
#[test]
fn test_find_roots_cubic_depressed() {
assert_eq!(find_roots_cubic_depressed(0f32, 0f32), Roots::One([0f32]));
assert_eq!(find_roots_cubic_depressed(-1f64, 0f64), Roots::Three([-1f64, 0f64, 1f64]));
match find_roots_cubic_depressed(-2f64, 2f64) {
Roots::One(x) => {
assert_float_array_eq!(1e-15, x, [-1.769292354238631415240409f64]);
}
_ => {
assert!(false);
}
}
match find_roots_cubic_depressed(-3f64, 2f64) {
Roots::Two(x) => {
assert_float_array_eq!(1e-15, x, [-2f64, 1f64]);
}
_ => {
assert!(false);
}
}
match find_roots_cubic_depressed(-2f64, 1f64) {
Roots::Three(x) => {
assert_float_array_eq!(1e-15, x, [(-1f64 - 5f64.sqrt()) / 2f64, (-1f64 + 5f64.sqrt()) / 2f64, 1f64]);
}
_ => {
assert!(false);
}
}
}
}