So where is the problem/bug?
I did construct the companion matrix manually, starting from the vector of
coefficients of the Wilkinson polynomial:
julia> pW = [
1.0, -210, 20615,
-1256850, 53327946, -1672280820,
40171771630, -756111184500, 11310276995381,
-135585182899530, 1307535010540395, -10142299865511450,
63030812099294896, -311333643161390656, 1206647803780373248,
-3599979517947607040, 8037811822645052416, -12870931245150988288,
13803759753640704000, -8752948036761600000, 2432902008176640000];
julia> pw = pW[2:21] / pW[1]; # normalize 1st
coefficient
julia> A = vcat(-pw', hcat(eye(19), zeros(19)); # companion matrix
julia> eigvals(A)
20-element Array{Float64,1}:
20.0003
18.9979
18.0072
16.9854
16.0183
14.9856
14.0067
12.9987
12.0002
10.9992
10.0009
8.99944
8.0002
6.99995
6.00001
5.0
4.0
3.0
2.0
1.0
This is almost the same result as with Matlab above. Does this mean,
Polynomial is constructing a different / wrong companion matrix?
On Monday, May 5, 2014 11:24:47 PM UTC+2, 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.
>
>
>
