In a program I wrote a while ago I needed to solve a cubic equation. The algorithm I wrote then used a binary search to find one solution then solved the resulting quadratic equation directly. While I was looking at the code today I wondered if somewhere on wiki there was code to solve a cubic equation directly, and I found one on Rodney Stephenson's page under a section he has for testing tcl performance. I thought the code deserved a page of its own.
# assumes that the first coefficient is always 1
proc cubic { a1 a2 a3 } {
set PI 3.1415926535897932385
set a12 [expr {$a1*$a1}]
set Q [expr {($a12-3*$a2)/9.}]
set R [expr {(2*$a1*$a12-9*$a1*$a2+27*$a3)/54.}]
set Z [expr {$Q*$Q*$Q - $R*$R}]
if {$Z == 0.0} {
if {$Q==0.0} {
set x0 [expr {-$a1/3.}]
return [list $x0]
} else {
if {$R<0} {
set x0 [expr {2*sqrt($Q)-$a1/3.}]
set x1 [expr {-sqrt($Q)-$a1/3.}]
} else {
set x0 [expr {sqrt($Q)- $a1/3.}]
set x1 [expr {-2*sqrt($Q)-$a1/3.}]
}
return [list $x0 $x1]
}
} elseif {$Z<0.0} {
set z2 [expr {pow(sqrt(-$Z) + abs($R), 1.0/3.0)}]
set z1 [expr {$z2 + $Q/$z2}]
if {$R>0} {set z1 [expr {-$z1}]}
set x0 [expr {$z1 - $a1/3.}]
return [list $x0]
} else {
set theta [expr {$R/sqrt($Q*$Q*$Q)}]
set Q2 [expr {-2.0*sqrt($Q)}]
set theta [expr {acos($theta)}]
set x0 [expr {$Q2*cos($theta/3.0)-$a1/3.}]
set x1 [expr {$Q2*cos(($theta+2*$PI)/3.0)-$a1/3.}]
set x2 [expr {$Q2*cos(($theta+4*$PI)/3.0)-$a1/3.}]
return [list $x0 $x1 $x2]
}
}