PolynomialRoots.jl icon indicating copy to clipboard operation
PolynomialRoots.jl copied to clipboard

Wrong results

Open saschatimme opened this issue 4 years ago • 8 comments

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.

saschatimme avatar Mar 05 '20 15:03 saschatimme