title: Generic function for solving n-order polynomial roots in Julia
body: <p>All, </p>
  
  <p>I've just been starting to play around with the Julia language and am enjoying it quite a bit. At the end of the 3rd tutorial there's an interesting problem: genericize the quadratic formula such that it solves for the roots of any <a href="http://forio.com/products/julia-studio/tutorials/beginner/3/" rel="nofollow">n-order polynomial equation</a>.</p>
  
  <p>This struck me as (a) an interesting programming problem and (b) an interesting Julia problem. Has anyone out there solved this one? For reference, here is the Julia code with a couple toy examples. Again, the idea is to make this generic for any n-order polynomial. </p>
  
  <p>Cheers,</p>
  
  <p>Aaron </p>
  
  <pre><code>function derivative(f)
      return function(x)
          # pick a small value for h
          h = x == 0 ? sqrt(eps(Float64)) : sqrt(eps(Float64)) * x
  
          # floating point arithmetic gymnastics
          xph = x + h
          dx = xph - x
  
          # evaluate f at x + h
          f1 = f(xph)
  
          # evaluate f at x
          f0 = f(x)
  
          # divide the difference by h
          return (f1 - f0) / dx
      end
  end
  
  
  function quadratic(f)
  
      f1 = derivative(f)
  
      c = f(0.0)
  
      b = f1(0.0)
  
      a = f(1.0) - b - c
  
      return (-b + sqrt(b^2 - 4a*c + 0im))/2a, (-b - sqrt(b^2 - 4a*c + 0im))/2a
  end
  
  quadratic((x) -&gt; x^2 - x - 2)
  quadratic((x) -&gt; x^2 + 2)
  </code></pre>