Expand description
Sorted and unique list of roots of an equation.
Variants§
No([F; 0])
Equation has no roots
One([F; 1])
Equation has one root (or all roots of the equation are the same)
Two([F; 2])
Equation has two roots
Three([F; 3])
Equation has three roots
Four([F; 4])
Equation has four roots
Implementations§
source§impl<F: FloatType> Roots<F>
impl<F: FloatType> Roots<F>
sourcepub fn add_new_root(self, new_root: F) -> Self
pub fn add_new_root(self, new_root: F) -> Self
Add a new root to existing ones keeping the list of roots ordered and unique.
Examples found in repository?
src/analytical/biquadratic.rs (line 53)
46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67
pub fn find_roots_biquadratic<F: FloatType>(a4: F, a2: F, a0: F) -> Roots<F> {
// Handle non-standard cases
if a4 == F::zero() {
// a4 = 0; a2*x^2 + a0 = 0; solve quadratic equation
super::quadratic::find_roots_quadratic(a2, F::zero(), a0)
} else if a0 == F::zero() {
// a0 = 0; a4*x^4 + a2*x^2 = 0; solve quadratic equation and add zero root
super::quadratic::find_roots_quadratic(a4, F::zero(), a2).add_new_root(F::zero())
} else {
// solve the corresponding quadratic equation and order roots
let mut roots = Roots::No([]);
for x in super::quadratic::find_roots_quadratic(a4, a2, a0).as_ref().iter() {
if *x > F::zero() {
let sqrt_x = x.sqrt();
roots = roots.add_new_root(-sqrt_x).add_new_root(sqrt_x);
} else if *x == F::zero() {
roots = roots.add_new_root(F::zero());
}
}
roots
}
}
More examples
src/analytical/quartic.rs (line 57)
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 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175
fn find_roots_via_depressed_quartic<F: FloatType>(a4: F, a3: F, a2: F, a1: F, a0: F, pp: F, rr: F, dd: F) -> Roots<F> {
// Depressed quartic
// https://en.wikipedia.org/wiki/Quartic_function#Converting_to_a_depressed_quartic
let _2 = F::from(2i16);
let _3 = F::from(3i16);
let _4 = F::from(4i16);
let _6 = F::from(6i16);
let _8 = F::from(8i16);
let _12 = F::from(12i16);
let _16 = F::from(16i16);
let _256 = F::from(256i16);
// a4*x^4 + a3*x^3 + a2*x^2 + a1*x + a0 = 0 => y^4 + p*y^2 + q*y + r.
let a4_pow_2 = a4 * a4;
let a4_pow_3 = a4_pow_2 * a4;
let a4_pow_4 = a4_pow_2 * a4_pow_2;
// Re-use pre-calculated values
let p = pp / (_8 * a4_pow_2);
let q = rr / (_8 * a4_pow_3);
let r = (dd + _16 * a4_pow_2 * (_12 * a0 * a4 - _3 * a1 * a3 + a2 * a2)) / (_256 * a4_pow_4);
let mut roots = Roots::No([]);
for y in super::quartic_depressed::find_roots_quartic_depressed(p, q, r)
.as_ref()
.iter()
{
roots = roots.add_new_root(*y - a3 / (_4 * a4));
}
roots
}
/// Solves a quartic equation a4*x^4 + a3*x^3 + a2*x^2 + a1*x + a0 = 0.
///
/// Returned roots are ordered.
/// Precision is about 5e-15 for f64, 5e-7 for f32.
/// WARNING: f32 is often not enough to find multiple roots.
///
/// # Examples
///
/// ```
/// use roots::find_roots_quartic;
///
/// let one_root = find_roots_quartic(1f64, 0f64, 0f64, 0f64, 0f64);
/// // Returns Roots::One([0f64]) as 'x^4 = 0' has one root 0
///
/// let two_roots = find_roots_quartic(1f32, 0f32, 0f32, 0f32, -1f32);
/// // Returns Roots::Two([-1f32, 1f32]) as 'x^4 - 1 = 0' has roots -1 and 1
///
/// let multiple_roots = find_roots_quartic(-14.0625f64, -3.75f64, 29.75f64, 4.0f64, -16.0f64);
/// // Returns Roots::Two([-1.1016116464173349f64, 0.9682783130840016f64])
///
/// let multiple_roots_not_found = find_roots_quartic(-14.0625f32, -3.75f32, 29.75f32, 4.0f32, -16.0f32);
/// // Returns Roots::No([]) because of f32 rounding errors when trying to calculate the discriminant
/// ```
pub fn find_roots_quartic<F: FloatType>(a4: F, a3: F, a2: F, a1: F, a0: F) -> Roots<F> {
// Handle non-standard cases
if a4 == F::zero() {
// a4 = 0; a3*x^3 + a2*x^2 + a1*x + a0 = 0; solve cubic equation
super::cubic::find_roots_cubic(a3, a2, a1, a0)
} else if a0 == F::zero() {
// a0 = 0; x^4 + a2*x^2 + a1*x = 0; reduce to cubic and arrange results
super::cubic::find_roots_cubic(a4, a3, a2, a1).add_new_root(F::zero())
} else if a1 == F::zero() && a3 == F::zero() {
// a1 = 0, a3 =0; a4*x^4 + a2*x^2 + a0 = 0; solve bi-quadratic equation
super::biquadratic::find_roots_biquadratic(a4, a2, a0)
} else {
let _3 = F::from(3i16);
let _4 = F::from(4i16);
let _6 = F::from(6i16);
let _8 = F::from(8i16);
let _9 = F::from(9i16);
let _10 = F::from(10i16);
let _12 = F::from(12i16);
let _16 = F::from(16i16);
let _18 = F::from(18i16);
let _27 = F::from(27i16);
let _64 = F::from(64i16);
let _72 = F::from(72i16);
let _80 = F::from(80i16);
let _128 = F::from(128i16);
let _144 = F::from(144i16);
let _192 = F::from(192i16);
let _256 = F::from(256i16);
// Discriminant
// https://en.wikipedia.org/wiki/Quartic_function#Nature_of_the_roots
// Partially simplifed to keep intermediate values smaller (to minimize rounding errors).
let discriminant = a4 * a0 * a4 * (_256 * a4 * a0 * a0 + a1 * (_144 * a2 * a1 - _192 * a3 * a0))
+ a4 * a0 * a2 * a2 * (_16 * a2 * a2 - _80 * a3 * a1 - _128 * a4 * a0)
+ (a3
* a3
* (a4 * a0 * (_144 * a2 * a0 - _6 * a1 * a1)
+ (a0 * (_18 * a3 * a2 * a1 - _27 * a3 * a3 * a0 - _4 * a2 * a2 * a2) + a1 * a1 * (a2 * a2 - _4 * a3 * a1))))
+ a4 * a1 * a1 * (_18 * a3 * a2 * a1 - _27 * a4 * a1 * a1 - _4 * a2 * a2 * a2);
let pp = _8 * a4 * a2 - _3 * a3 * a3;
let rr = a3 * a3 * a3 + _8 * a4 * a4 * a1 - _4 * a4 * a3 * a2;
let delta0 = a2 * a2 - _3 * a3 * a1 + _12 * a4 * a0;
let dd = _64 * a4 * a4 * a4 * a0 - _16 * a4 * a4 * a2 * a2 + _16 * a4 * a3 * a3 * a2
- _16 * a4 * a4 * a3 * a1
- _3 * a3 * a3 * a3 * a3;
// Handle special cases
let double_root = discriminant == F::zero();
if double_root {
let triple_root = double_root && delta0 == F::zero();
let quadruple_root = triple_root && dd == F::zero();
let no_roots = dd == F::zero() && pp > F::zero() && rr == F::zero();
if quadruple_root {
// Wiki: all four roots are equal
Roots::One([-a3 / (_4 * a4)])
} else if triple_root {
// Wiki: At least three roots are equal to each other
// x0 is the unique root of the remainder of the Euclidean division of the quartic by its second derivative
//
// Solved by SymPy (ra is the desired reminder)
// a, b, c, d, e = symbols('a,b,c,d,e')
// f=a*x**4+b*x**3+c*x**2+d*x+e // Quartic polynom
// g=6*a*x**2+3*b*x+c // Second derivative
// q, r = div(f, g) // SymPy only finds the highest power
// simplify(f-(q*g+r)) == 0 // Verify the first division
// qa, ra = div(r/a,g/a) // Workaround to get the second division
// simplify(f-((q+qa)*g+ra*a)) == 0 // Verify the second division
// solve(ra,x)
// ----- yields
// (−72*a^2*e+10*a*c^2−3*b^2*c)/(9*(8*a^2*d−4*a*b*c+b^3))
let x0 = (-_72 * a4 * a4 * a0 + _10 * a4 * a2 * a2 - _3 * a3 * a3 * a2)
/ (_9 * (_8 * a4 * a4 * a1 - _4 * a4 * a3 * a2 + a3 * a3 * a3));
let roots = Roots::One([x0]);
roots.add_new_root(-(a3 / a4 + _3 * x0))
} else if no_roots {
// Wiki: two complex conjugate double roots
Roots::No([])
} else {
find_roots_via_depressed_quartic(a4, a3, a2, a1, a0, pp, rr, dd)
}
} else {
let no_roots = discriminant > F::zero() && (pp > F::zero() || dd > F::zero());
if no_roots {
// Wiki: two pairs of non-real complex conjugate roots
Roots::No([])
} else {
find_roots_via_depressed_quartic(a4, a3, a2, a1, a0, pp, rr, dd)
}
}
}
}
src/analytical/cubic_normalized.rs (line 65)
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
pub fn find_roots_cubic_normalized<F: FloatType>(a2: F, 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);
let q = (_3 * a1 - a2 * a2) / _9;
let r = (_9 * a2 * a1 - _27 * a0 - _2 * a2 * a2 * a2) / _54;
let q3 = q * q * q;
let d = q3 + r * r;
let a2_div_3 = a2 / _3;
if d < F::zero() {
let phi_3 = (r / (-q3).sqrt()).acos() / _3;
let sqrt_q_2 = _2 * (-q).sqrt();
Roots::One([sqrt_q_2 * phi_3.cos() - a2_div_3])
.add_new_root(sqrt_q_2 * (phi_3 - F::two_third_pi()).cos() - a2_div_3)
.add_new_root(sqrt_q_2 * (phi_3 + F::two_third_pi()).cos() - a2_div_3)
} else {
let sqrt_d = d.sqrt();
let s = (r + sqrt_d).cbrt();
let t = (r - sqrt_d).cbrt();
if s == t {
if s + t == F::zero() {
Roots::One([s + t - a2_div_3])
} else {
Roots::One([s + t - a2_div_3]).add_new_root(-(s + t) / _2 - a2_div_3)
}
} else {
Roots::One([s + t - a2_div_3])
}
}
}
src/analytical/cubic_depressed.rs (line 55)
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
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])
}
}
}
}
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([])
}
}
}
src/analytical/cubic.rs (line 98)
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)
}
}
}