PolynomialRoots.jl
PolynomialRoots.jl copied to clipboard
Wrong results
I did some test computations and found the following troubling result:
julia> using DynamicPolynomials
julia> using PolynomialRoots
julia> function poly(roots)
@polyvar x
d = length(roots)
f = prod(x - xi for xi in roots)
[coefficient(f, x^i) for i in 0:d]
end
poly (generic function with 1 method)
julia> d = 20
20
julia> c = [randn(ComplexF64); randn(ComplexF64, d - 1) * 1e-4; 1]
21-element Array{Complex{Float64},1}:
0.6038117420721143 - 0.4692317923460019im
8.677470319492543e-5 - 1.9856488102339103e-5im
-2.8145901548439397e-5 + 2.5958816494607353e-6im
-9.933810484391282e-5 - 0.0001480730021139855im
-8.820204766296887e-5 - 0.0001042484315957653im
7.88021087360911e-5 + 2.5513881910383407e-5im
-5.007828702552179e-5 - 4.850423221617593e-7im
3.788521654994467e-5 + 2.0944404532776053e-5im
-1.1315308780483735e-5 - 1.8001817822660527e-5im
0.00011608934805185687 + 0.0001503409434982439im
4.125906444101907e-5 + 9.978779625362144e-5im
-5.52502458307844e-5 + 1.244622053838063e-5im
3.170867062110031e-5 - 1.9111588344365587e-5im
9.754810243956861e-5 + 6.369172830914271e-5im
7.028073862086189e-5 - 3.3819410763212794e-5im
0.0001238375611115403 - 5.144632426358618e-5im
-5.263257681962449e-5 + 6.227546399484231e-5im
4.6728130048484104e-5 - 0.00011791878030739642im
0.00010423858969876596 - 0.00011299961775497805im
-8.557623396259921e-5 - 2.825091927864975e-5im
1.0 + 0.0im
julia> r = roots(c)
20-element Array{Complex{Float64},1}:
7.802412795926557 - 4.178757163336667im
-1.2054492306931206 + 8.768496484683855im
-8.768489889160787 - 1.2054534486770816im
-4.178754803019869 - 7.80240884258056im
8.768498428644229 + 1.2054564365850402im
1.2054578038795405 - 8.768493818879023im
6.385317333146652 + 6.129227118201205im
-6.129219730795856 + 6.38531432406573im
3.856072372001663 - 7.966827038663316im
-7.9668235583329405 - 3.856068311866256im
-7.802404396639403 + 4.1787600200142485im
4.178763518989671 + 7.802411631716743im
6.129228142519897 - 6.385311565395753im
-6.385308625464333 - 6.129224228027822im
8.71184335999314 - 1.5631567624600453im
1.563163615363657 + 8.711841873146625im
-3.856063894157322 + 7.966829725325496im
7.966832198295311 + 3.856071278929593im
-1.5631549515439616 - 8.711839169108076im
-8.711834912718764 + 1.563159707245344im
julia> c - poly(r)
21-element Array{Complex{Float64},1}:
-7.987123734845452e18 + 3.4639346941464863e18im
896.0000867747032 + 1151.999980143512im
-16.00002814590155 + 18.00000259588165im
-8.000099338104844 - 14.500148073002114im
0.24991179795233703 - 0.18760424843159576im
0.1563288021087361 + 0.17190051388191038im
-0.007862578287025522 - 4.850423221617593e-7im
-0.00020625540845005532 + 0.000753366279532776im
-1.1315308780483735e-5 - 0.00020110728657266052im
2.4536613676856873e-5 - 2.513513072050611e-5im
2.6352546265659437e-6 + 1.5593416637776903e-6im
-2.6542446384986325e-8 + 1.0805906987477109e-7im
-1.7764205750584404e-8 + 6.601467340955743e-9im
-9.523320931490147e-10 - 4.682197821015313e-9im
5.2134349463363316e-11 - 3.523188618357654e-10im
5.818546518080239e-12 + 2.4179138238041714e-13im
5.263254309923089e-12 + 7.162041142441677e-13im
1.88754143337546e-13 - 2.443320881988925e-13im
-4.330127294054076e-14 - 1.3839515471125718e-14im
-2.160977560089483e-15 - 1.5989373182683647e-15im
0.0 + 0.0im
I am not sure whether this is a problem with the algorithm or the specific implementation.