Keith Vetter 2003-02-26 -- 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]
}
}escargo 26 Feb 2003 - Is it worthwhile to manually factor out common subexpressions like sqrt($Q) and $a1/3. that occur in multiple places? Is it possible for the byte compiler to do optimizations across multiple expr evaluations?
DKF - It's theoretically possible, but we don't. And it's not in general possible to factor out any read of a variable due to traces.
KPV - I tried comparing the speed of the above code with tweaked code that factors out the common subexpressions noted above. The result on my 800mhz Pentium was that both took exactly the same amount of time. Looking at the code explains why: factoring out sqrt($Q) only saves you 1 sqrt call and only when there are two solutions; factoring out $a1/3. only saves you one division when there are two solutions and two divisions when there are three solutions. Since each call to cubic has about forty floating point operations, including several transcendental evaluations, this optimization is lost in the noise.
CL points out to casual readers that, as Keith and Donal know quite well, performance is far from the only reason to factor common sub-expressions. It's at least as importance for maintainability. [Find reference on refactoring.]
DKF points out that in this particular case it doesn't help too much anyway; the code is only truly clear if you understand how to go about solving cubic equations in the first place... ;^)