Function roots::find_roots_cubic_normalized

source · 
pub fn find_roots_cubic_normalized<F: FloatType>(a2: F, a1: F, a0: F) -> Roots<F>
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
Hide additional 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)
        }
    }
}