<p>The package <a href="https://github.com/giordano/PolynomialRoots.jl" rel="nofollow"><code>PolynomialRoots.jl</code></a> provides the function <code>roots()</code> to find all (real and complex) roots of polynomials of any order.  The only mandatory argument is the array with coefficients of the polynomial in ascending order.</p>

<p>For example, in order to find the roots of</p>

<pre><code>6x^5 + 5x^4 + 3x^2 + 2x + 1
</code></pre>

<p>after loading the package (<code>using PolynomialRoots</code>) you can use</p>



<pre class="lang-julia prettyprint-override"><code>julia&gt; roots([1, 2, 3, 4, 5, 6])
5-element Array{Complex{Float64},1}:
           0.294195-0.668367im
 -0.670332+2.77556e-17im
           0.294195+0.668367im
          -0.375695-0.570175im
          -0.375695+0.570175im
</code></pre>

<p>The package is a Julia implementation of the root-finding algorithm described in this paper: <a href="http://arxiv.org/abs/1203.1034" rel="nofollow">http://arxiv.org/abs/1203.1034</a></p>

<p><code>PolynomialRoots.jl</code> has also support for arbitrary precision calculation.  This is useful for solving equation that cannot be solved in double precision.  For example</p>

<pre class="lang-julia prettyprint-override"><code>julia&gt; r = roots([94906268.375, -189812534, 94906265.625]);

julia&gt; (r[1], r[2])
(1.0000000144879793 - 0.0im,1.0000000144879788 + 0.0im)
</code></pre>

<p>gives the wrong result for the polynomial, instead passing the input array in arbitrary precision forces arbitrary precision calculations that provide the right answer (see <a href="https://en.wikipedia.org/wiki/Loss_of_significance" rel="nofollow">https://en.wikipedia.org/wiki/Loss_of_significance</a>):</p>

<pre class="lang-julia prettyprint-override"><code>julia&gt; r = roots([BigFloat(94906268.375), BigFloat(-189812534), BigFloat(94906265.625)]);

julia&gt; (Float64(r[1]), Float64(r[2]))
(1.0000000289759583,1.0)
</code></pre>
Generic function for solving n-order polynomial roots in Julia
<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>