Enum roots::Roots

source · 
pub enum Roots<F: FloatType> {
    No([F; 0]),
    One([F; 1]),
    Two([F; 2]),
    Three([F; 3]),
    Four([F; 4]),
}
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§

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
Hide additional 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)
        }
    }
}

Trait Implementations§

Converts this type into a shared reference of the (usually inferred) input type.
Formats the value using the given formatter. Read more
This method tests for self and other values to be equal, and is used by ==.
This method tests for !=. The default implementation is almost always sufficient, and should not be overridden without very good reason.

Auto Trait Implementations§

Blanket Implementations§

Gets the TypeId of self. Read more
Immutably borrows from an owned value. Read more
Mutably borrows from an owned value. Read more

Returns the argument unchanged.

Calls U::from(self).

That is, this conversion is whatever the implementation of From<T> for U chooses to do.

The type returned in the event of a conversion error.
Performs the conversion.
The type returned in the event of a conversion error.
Performs the conversion.