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]
    }
 }

----

[Category Mathematics]