@@ -45,16 +45,16 @@ fn find_roots_via_depressed_quartic<F: FloatType>(a4: F, a3: F, a2: F, a1: F, a0
45
45
let a4_pow_3 = a4_pow_2 * a4;
46
46
let a4_pow_4 = a4_pow_2 * a4_pow_2;
47
47
// Re-use pre-calculated values
48
- let p = pp / ( _8 * dbg ! ( a4_pow_2) ) ;
49
- let q = rr / ( _8 * dbg ! ( a4_pow_3) ) ;
50
- let r = ( dd + _16 * a4_pow_2 * ( _12 * a0 * a4 - _3 * a1 * a3 + a2 * a2) ) / ( _256 * dbg ! ( a4_pow_4) ) ;
48
+ let p = pp / ( _8 * a4_pow_2) ;
49
+ let q = rr / ( _8 * a4_pow_3) ;
50
+ let r = ( dd + _16 * a4_pow_2 * ( _12 * a0 * a4 - _3 * a1 * a3 + a2 * a2) ) / ( _256 * a4_pow_4) ;
51
51
52
52
let mut roots = Roots :: No ( [ ] ) ;
53
- for y in super :: quartic_depressed:: find_roots_quartic_depressed ( dbg ! ( p ) , dbg ! ( q ) , dbg ! ( r ) )
53
+ for y in super :: quartic_depressed:: find_roots_quartic_depressed ( p , q , r )
54
54
. as_ref ( )
55
55
. iter ( )
56
56
{
57
- roots = roots. add_new_root ( dbg ! ( * y ) - a3 / ( _4 * a4) ) ;
57
+ roots = roots. add_new_root ( * y - a3 / ( _4 * a4) ) ;
58
58
}
59
59
roots
60
60
}
@@ -75,10 +75,10 @@ fn find_roots_via_depressed_quartic<F: FloatType>(a4: F, a3: F, a2: F, a1: F, a0
75
75
///
76
76
/// let two_roots = find_roots_quartic(1f32, 0f32, 0f32, 0f32, -1f32);
77
77
/// // Returns Roots::Two([-1f32, 1f32]) as 'x^4 - 1 = 0' has roots -1 and 1
78
- ///
78
+ ///
79
79
/// let multiple_roots = find_roots_quartic(-14.0625f64, -3.75f64, 29.75f64, 4.0f64, -16.0f64);
80
80
/// // Returns Roots::Two([-1.1016116464173349f64, 0.9682783130840016f64])
81
- ///
81
+ ///
82
82
/// let multiple_roots_not_found = find_roots_quartic(-14.0625f32, -3.75f32, 29.75f32, 4.0f32, -16.0f32);
83
83
/// // Returns Roots::No([]) because of f32 rounding errors when trying to calculate the discriminant
84
84
/// ```
@@ -114,15 +114,13 @@ pub fn find_roots_quartic<F: FloatType>(a4: F, a3: F, a2: F, a1: F, a0: F) -> Ro
114
114
// Discriminant
115
115
// https://en.wikipedia.org/wiki/Quartic_function#Nature_of_the_roots
116
116
// Partially simplifed to keep intermediate values smaller (to minimize rounding errors).
117
- let discriminant = dbg ! ( a4 * a0 * a4 * ( _256 * a4 * a0 * a0 + a1 * ( _144 * a2 * dbg!( a1) - _192 * dbg!( a3) * a0) ) )
118
- + dbg ! ( a4 * a0 * a2 * a2 * ( _16 * a2 * a2 - _80 * a3 * a1 - _128 * a4 * a0) )
119
- + dbg ! (
120
- a3 * a3
121
- * ( a4 * a0 * ( _144 * a2 * a0 - _6 * a1 * a1)
122
- + ( a0 * ( _18 * a3 * a2 * a1 - _27 * a3 * a3 * a0 - _4 * a2 * a2 * a2)
123
- + a1 * a1 * ( a2 * a2 - _4 * a3 * a1) ) )
124
- )
125
- + dbg ! ( a4 * a1 * a1 * ( _18 * a3 * a2 * a1 - _27 * a4 * a1 * a1 - _4 * a2 * a2 * a2) ) ;
117
+ let discriminant = a4 * a0 * a4 * ( _256 * a4 * a0 * a0 + a1 * ( _144 * a2 * a1 - _192 * a3 * a0) )
118
+ + a4 * a0 * a2 * a2 * ( _16 * a2 * a2 - _80 * a3 * a1 - _128 * a4 * a0)
119
+ + ( a3
120
+ * a3
121
+ * ( a4 * a0 * ( _144 * a2 * a0 - _6 * a1 * a1)
122
+ + ( a0 * ( _18 * a3 * a2 * a1 - _27 * a3 * a3 * a0 - _4 * a2 * a2 * a2) + a1 * a1 * ( a2 * a2 - _4 * a3 * a1) ) ) )
123
+ + a4 * a1 * a1 * ( _18 * a3 * a2 * a1 - _27 * a4 * a1 * a1 - _4 * a2 * a2 * a2) ;
126
124
let pp = _8 * a4 * a2 - _3 * a3 * a3;
127
125
let rr = a3 * a3 * a3 + _8 * a4 * a4 * a1 - _4 * a4 * a3 * a2;
128
126
let delta0 = a2 * a2 - _3 * a3 * a1 + _12 * a4 * a0;
@@ -131,15 +129,15 @@ pub fn find_roots_quartic<F: FloatType>(a4: F, a3: F, a2: F, a1: F, a0: F) -> Ro
131
129
- _3 * a3 * a3 * a3 * a3;
132
130
133
131
// Handle special cases
134
- let double_root = dbg ! ( discriminant) == F :: zero ( ) ;
135
- if dbg ! ( double_root) {
136
- let triple_root = double_root && dbg ! ( delta0) == F :: zero ( ) ;
137
- let quadruple_root = triple_root && dbg ! ( dd ) == F :: zero ( ) ;
138
- let no_roots = dd == F :: zero ( ) && dbg ! ( pp ) > F :: zero ( ) && dbg ! ( rr ) == F :: zero ( ) ;
139
- if dbg ! ( quadruple_root) {
132
+ let double_root = discriminant == F :: zero ( ) ;
133
+ if double_root {
134
+ let triple_root = double_root && delta0 == F :: zero ( ) ;
135
+ let quadruple_root = triple_root && dd == F :: zero ( ) ;
136
+ let no_roots = dd == F :: zero ( ) && pp > F :: zero ( ) && rr == F :: zero ( ) ;
137
+ if quadruple_root {
140
138
// Wiki: all four roots are equal
141
139
Roots :: One ( [ -a3 / ( _4 * a4) ] )
142
- } else if dbg ! ( triple_root) {
140
+ } else if triple_root {
143
141
// Wiki: At least three roots are equal to each other
144
142
// x0 is the unique root of the remainder of the Euclidean division of the quartic by its second derivative
145
143
//
@@ -156,21 +154,21 @@ pub fn find_roots_quartic<F: FloatType>(a4: F, a3: F, a2: F, a1: F, a0: F) -> Ro
156
154
// (−72*a^2*e+10*a*c^2−3*b^2*c)/(9*(8*a^2*d−4*a*b*c+b^3))
157
155
let x0 = ( -_72 * a4 * a4 * a0 + _10 * a4 * a2 * a2 - _3 * a3 * a3 * a2)
158
156
/ ( _9 * ( _8 * a4 * a4 * a1 - _4 * a4 * a3 * a2 + a3 * a3 * a3) ) ;
159
- let roots = dbg ! ( Roots :: One ( [ x0] ) ) ;
160
- roots. add_new_root ( dbg ! ( -( a3 / a4 + _3 * x0) ) )
161
- } else if dbg ! ( no_roots) {
157
+ let roots = Roots :: One ( [ x0] ) ;
158
+ roots. add_new_root ( -( a3 / a4 + _3 * x0) )
159
+ } else if no_roots {
162
160
// Wiki: two complex conjugate double roots
163
161
Roots :: No ( [ ] )
164
162
} else {
165
- find_roots_via_depressed_quartic ( a4, a3, a2, a1, a0, pp, dbg ! ( rr ) , dd)
163
+ find_roots_via_depressed_quartic ( a4, a3, a2, a1, a0, pp, rr , dd)
166
164
}
167
165
} else {
168
- let no_roots = dbg ! ( pp ) > F :: zero ( ) || dbg ! ( dd ) > F :: zero ( ) ;
169
- if dbg ! ( no_roots) {
166
+ let no_roots = pp > F :: zero ( ) || dd > F :: zero ( ) ;
167
+ if no_roots {
170
168
// Wiki: two pairs of non-real complex conjugate roots
171
169
Roots :: No ( [ ] )
172
170
} else {
173
- find_roots_via_depressed_quartic ( a4, a3, a2, a1, a0, pp, dbg ! ( rr ) , dd)
171
+ find_roots_via_depressed_quartic ( a4, a3, a2, a1, a0, pp, rr , dd)
174
172
}
175
173
}
176
174
}
@@ -222,7 +220,13 @@ mod test {
222
220
) ;
223
221
// ... even after normalizing
224
222
assert_eq ! (
225
- find_roots_quartic( 1f32 , -3.75f32 /-14.0625f32 , 29.75f32 /-14.0625f32 , 4.0f32 /-14.0625f32 , -16.0f32 /-14.0625f32 ) ,
223
+ find_roots_quartic(
224
+ 1f32 ,
225
+ -3.75f32 / -14.0625f32 ,
226
+ 29.75f32 / -14.0625f32 ,
227
+ 4.0f32 / -14.0625f32 ,
228
+ -16.0f32 / -14.0625f32
229
+ ) ,
226
230
Roots :: No ( [ ] )
227
231
) ;
228
232
// assert_eq!(
0 commit comments