In numerical computation, how could increasing the number of bits of precision decrease accuracy?

Offtopic
5

A more appropriate test would seem to be:

  1. generating random Complex64 polynomial coefficients c
  2. convert to BigFloat (an exact conversion) and find roots r = roots(big.(c))
  3. construct new BigFloat polynomial c2 coefficients by multiplying monomials
  4. compare to c.

This should match if you have sufficiently high precision. However, I find that it gives a wildly wrong error with higher and higher probability as the degree increases:

using PolynomialRoots, DynamicPolynomials

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

d = 5
c = [randn(ComplexF64); randn(ComplexF64, d - 1) * 1e-4; 1]
setprecision(8192)
r = roots(big.(c), polish=true)
c2 = ComplexF64.(poly(r))
c2 - c

90% of the time gives (correctly)

6-element Vector{ComplexF64}:
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im

but if we increase the degree to d = 6 it is wildy wrong about 85% of the time (≈ 15% of the time it correctly gives zero), giving answers for c - c2 like

7-element Vector{ComplexF64}:
 3.9439141771536586e6 - 7.66389580196865e6im
 2.2082351722049564e6 - 2.857913579123853e6im
    469737.1655469684 - 420674.6692220386im
    50125.96985319289 - 30578.538953102627im
    2872.684011890461 - 1098.2968892416864im
     84.4854525989318 - 15.599756742987397im
                  0.0 + 0.0im

and if you check with evalpoly you’ll find that these aren’t roots at all, and indeed p(r) seems to have a constant(!!) offset:

julia> ComplexF64.(evalpoly.(r, Ref(c)))
6-element Vector{ComplexF64}:
 -187275.93508116304 + 125277.6562628678im
 -187275.93508116304 + 125277.6562628678im
 -187275.93508116304 + 125277.6562628678im
 -187275.93508116304 + 125277.6562628678im
 -187275.93508116304 + 125277.6562628678im
 -187275.93508116304 + 125277.6562628678im

In fact, if you look at the “roots” they are all the same, which is clearly wrong:

julia> ComplexF64.(r)
6-element Vector{ComplexF64}:
 -0.4426969318472816 - 6.30885464214318im
 -0.4426969318472816 - 6.30885464214318im
 -0.4426969318472816 - 6.30885464214318im
 -0.4426969318472816 - 6.30885464214318im
 -0.4426969318472816 - 6.30885464214318im
 -0.4426969318472816 - 6.30885464214318im

Conclusion: PolynomialRoots appears to be returning garbage roots some of the time. (This also falls under the category of “you aren’t computing what you think you are”.)

3 Likes