This solves a cubic equation of the form:
f(x) = ax^3 + bx^2 + cx + d
The result is returned in a dict containing the three solutions, root1, root2 and root3, consisting of a real part and an imaginary part.
Was transcoded from javascript from http://www.random-science-tools.com/maths/cubic.htm
KPV see also Cubic equation
proc cubic_root x {
if {$x < 0} {
return [expr {-pow(- $x,(1.0/3.0))}]
} else {
return [expr {pow($x,(1.0/3.0))}]
}
}
proc solveCubic {a b c d} {
# solve a cubic equation
set rt {root1 {re 0 im 0} root2 {re 0 im 0} root3 {re 0 im 0}}
set q [expr {double(3 * $a * $c - $b * $b) / (9 * $a * $a)}]
set r [expr {double(9 * $a * $b * $c - 27 * $a * $a * $d - 2 * $b * $b * $b) / (54 * $a * $a * $a)}]
set qqq_plus_rr [expr { $q * $q * $q + $r * $r}]
if {$qqq_plus_rr < 0} {
# All three roots real and unequal
set j [expr {sqrt(- $q)}]
set i [expr {$j * $j * $j}]
set k [expr {acos(double($r) / $i)}]
set M [expr {cos($k/3.0)}]
set N [expr {sqrt(3) * sin($k/3.0)}]
set b_over_3a [expr {$b/(3.0 * $a)}]
dict set rt root1 re [expr {2 * $j * cos($k / 3.0) - $b_over_3a}]
dict set rt root2 re [expr {- $j * ($M + $N) - $b_over_3a}]
dict set rt root3 re [expr {$j * ($N - $M) - $b_over_3a}]
} else {
if {$qqq_plus_rr > 0} {
set root_qqq_plus_rr [expr {sqrt($qqq_plus_rr)}]
set s [cubic_root [expr {$r + $root_qqq_plus_rr}]]
set t [cubic_root [expr {$r - $root_qqq_plus_rr}]]
dict set rt root1 re [expr {$s + $t - $b / (3.0 * $a)}]
dict set rt root2 re [expr {-0.5 * ($s + $t) - $b / (3.0 * $a)}]
dict set rt root2 im [expr {sqrt(3) / 2.0 * ($s - $t)}]
dict set rt root3 re [dict get $rt root2 re]
dict set rt root3 im [expr {-[dict get $rt root2 im]}]
} else {
# ppp == rr so three real equal roots
set res [expr {-[cubic_root [expr {double($d) / $a}]]}]
dict set rt root1 re $res
dict set rt root2 re $res
dict set rt root3 re $res
}
}
return $rt;
} # usage:
solveCubic 28 5 -5 11
# -> root1 {re -0.8837041642689161 im 0} root2 {re 0.35256636784874373 im 0.5659101188251707} root3 {re 0.35256636784874373 im -0.5659101188251707}