Updated 2013-07-17 20:31:13 by dzach

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}