Function roots::find_roots_cubic_normalized
source · Expand description
Solves a normalized cubic equation x^3 + a2x^2 + a1x + a0 = 0.
Trigonometric solution (arccos/cos) is implemented for three roots.
In case more than one roots are present, they are returned in the increasing order.
Examples
use roots::find_roots_cubic_normalized;
let one_root = find_roots_cubic_normalized(0f64, 0f64, 0f64);
// Returns Roots::One([0f64]) as 'x^3 = 0' has one root 0
let three_roots = find_roots_cubic_normalized(0f32, -1f32, 0f32);
// Returns Roots::Three([-1f32, -0f32, 1f32]) as 'x^3 - x = 0' has roots -1, 0, and 1
Examples found in repository?
src/analytical/quartic_depressed.rs (line 50)
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
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([])
}
}
}
More examples
src/analytical/cubic.rs (line 69)
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 124 125
pub fn find_roots_cubic<F: FloatType>(a3: F, a2: F, a1: F, a0: F) -> Roots<F> {
// Handle non-standard cases
if a3 == F::zero() {
// a3 = 0; a2*x^2+a1*x+a0=0; solve quadratic equation
super::quadratic::find_roots_quadratic(a2, a1, a0)
} else if a2 == F::zero() {
// a2 = 0; a3*x^3+a1*x+a0=0; solve depressed cubic equation
super::cubic_depressed::find_roots_cubic_depressed(a1 / a3, a0 / a3)
} else if a3 == F::one() {
// solve normalized cubic expression
super::cubic_normalized::find_roots_cubic_normalized(a2, a1, a0)
} else {
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);
// standard case
let d = _18 * a3 * a2 * a1 * a0 - _4 * a2 * a2 * a2 * a0 + a2 * a2 * a1 * a1
- _4 * a3 * a1 * a1 * a1
- _27 * a3 * a3 * a0 * a0;
let d0 = a2 * a2 - _3 * a3 * a1;
let d1 = _2 * a2 * a2 * a2 - _9 * a3 * a2 * a1 + _27 * a3 * a3 * a0;
if d < F::zero() {
// one real root
let sqrt = (-_27 * a3 * a3 * d).sqrt();
let c = F::cbrt(if d1 < F::zero() { d1 - sqrt } else { d1 + sqrt } / _2);
let x = -(a2 + c + d0 / c) / (_3 * a3);
Roots::One([x])
} else if d == F::zero() {
// multiple roots
if d0 == F::zero() {
// triple root
Roots::One([-a2 / (a3 * _3)])
} else {
// single root and double root
Roots::One([(_9 * a3 * a0 - a2 * a1) / (d0 * _2)])
.add_new_root((_4 * a3 * a2 * a1 - _9 * a3 * a3 * a0 - a2 * a2 * a2) / (a3 * d0))
}
} else {
// three real roots
let c3_img = F::sqrt(_27 * a3 * a3 * d) / _2;
let c3_real = d1 / _2;
let c3_module = F::sqrt(c3_img * c3_img + c3_real * c3_real);
let c3_phase = _2 * F::atan(c3_img / (c3_real + c3_module));
let c_module = F::cbrt(c3_module);
let c_phase = c3_phase / _3;
let c_real = c_module * F::cos(c_phase);
let c_img = c_module * F::sin(c_phase);
let x0_real = -(a2 + c_real + (d0 * c_real) / (c_module * c_module)) / (_3 * a3);
let e_real = -F::one() / _2;
let e_img = F::sqrt(_3) / _2;
let c1_real = c_real * e_real - c_img * e_img;
let c1_img = c_real * e_img + c_img * e_real;
let x1_real = -(a2 + c1_real + (d0 * c1_real) / (c1_real * c1_real + c1_img * c1_img)) / (_3 * a3);
let c2_real = c1_real * e_real - c1_img * e_img;
let c2_img = c1_real * e_img + c1_img * e_real;
let x2_real = -(a2 + c2_real + (d0 * c2_real) / (c2_real * c2_real + c2_img * c2_img)) / (_3 * a3);
Roots::One([x0_real]).add_new_root(x1_real).add_new_root(x2_real)
}
}
}