Not just 32 and 64, these coefficients, as produced through poly([1:20]) or poly([1.0:20.0]), aren't accurate as integers or floating point at 64 bit. You can see this by comparing `poly(big(1):big(20)).a` with `big(1.0:20.0).a`, say. Perhaps it is interesting, the differences in the coefficients found here (http://blogs.mathworks.com/cleve/2013/03/04/wilkinsons-polynomials/) differ in one term (the x^16 coefficient) from those of julia on my macbook pro. This small perturbation might be part of the difference you see as well.
Now just for fun, the bisection method can get this right using big integers,as you can see with `using Roots; fzeros(x->polyval(poly(big(1):big(20)),x), 0, 21)`. (It kinda gets lucky, fzeros just splits the interval into a few hundred pieces and looks for bracketing intervals.) On Monday, May 5, 2014 5:24:47 PM UTC-4, Tony Kelman wrote: > > Good catch. This looks almost entirely due to the order in which > Polynomial.jl is creating the companion matrix. The modified version at > https://github.com/loladiro/Polynomials.jl which stores the coefficients > in the mathematical order instead of the Matlab order does a little better, > but worryingly I'm getting noticeably different results between 32 bit and > 64 bit. > > Just transposing the companion matrix gives much better-behaved results > for this polynomial, not exactly the same as Matlab but pretty close. > > > > On Monday, May 5, 2014 1:14:44 PM UTC-7, Hans W Borchers wrote: >> >> The following is the so-called Wilkinson polynomial, that has zeros at 1, >> 2, ..., 20 : >> >> julia> using Polynomial >> >> >> julia> pW = poly([1.0:20.0]) >> >> >> julia> roots(pW) >> 20-element Array{Complex{Float64},1}: >> 1.0-0.0im >> 2.0-0.0im >> 3.0-0.0im >> 3.99999-0.0im >> 5.00006-0.0im >> 6.00181-0.0im >> 6.95031-0.0im >> 8.09172-0.579086im >> 8.09172+0.579086im >> 9.70891-1.64679im >> 9.70891+1.64679im >> 11.7999-2.59604im >> 11.7999+2.59604im >> 14.3615-3.12883im >> 14.3615+3.12883im >> 17.1543-2.89387im >> 17.1543+2.89387im >> 19.5758-1.70785im >> 19.5758+1.70785im >> 20.6635-0.0im >> >> >> Compare this with the same computation in Matlab: >> >> >> pw = poly(1:20); >> >> roots(pw) >> >> ans = >> >> 20.0003 >> 18.9972 >> 18.0112 >> 16.9711 >> 16.0483 >> 14.9354 >> 14.0653 >> 12.9491 >> 12.0334 >> 10.9840 >> 10.0061 >> 8.9984 >> 8.0003 >> 7.0000 >> 6.0000 >> 5.0000 >> 4.0000 >> 3.0000 >> 2.0000 >> 1.0000 >> >> >> The roots of this polynomial are known to be difficult to determine, but >> the results in Julia's Polynomial package seem really far off. >> >
