What’s the right way to QuickCheck floating-point routines?

I got a lot of great comments on my previous post about finding roots of polynomials in Haskell. One particularly promising idea I got from commenter Jake was to give up on the idea of having no external dependencies (which, to be fair, in these days of stack and nix and cabal-v2, seems like much less of a big deal than it used to be), and use the hmatrix package to find the eigenvalues of the companion matrix, which are exactly the roots.

So I tried that, and it seems to work great! The only problem is that I still don’t know how to write a reasonable test suite. I started by making a QuickCheck property expressing the fact that if we evaluate a polynomial at the returned roots, we should get something close to zero. I evaluate the polynomial using Horner’s method, which as far as I understand has good numerical stability in addition to being efficient.

polyRoots :: [Double] -> [Double]
polyRoots = ... stuff using hmatrix, see code at end of post ...

horner :: [Double] -> Double -> Double
horner as x = foldl' (\r a -> r*x + a) 0 as

_polyRoots_prop :: [Double] -> Property
_polyRoots_prop as = (length as > 1) ==>
  all ((< 1e-10) . abs . horner as) (polyRoots as)

This property passes 100 tests for quadratic polynomials, but for cubic I get failures; here’s an example. Consider the polynomial

Finding its roots via hmatrix yields three:

[-1.077801388041068, 0.5106483227001805, 150.6238979010295]

Notice that the third root is much bigger in magnitude than the other two, and indeed, that third root is the problematic one. Evaluating the polynomial at these roots via Horner’s method yields

[1.2434497875801753e-14, 1.7763568394002505e-15, -1.1008971512183052e-10]

the third of which is bigger than 1e-10 which I had (very arbitrarily!) chosen as the cutoff for “close enough to zero”. But here’s the thing: after playing around with it a bit, it seems like this is the most accurate possible value for the root that can be represented using Double. That is, if I evaluate the polynomial at any value other than the root that was returned—even if I just change the very last digit by 1 in either direction—I get a result which is farther from zero.

If I make the magic cutoff value bigger—say, 1e-8 instead of 1e-10—then I still get similar counterexamples, but for larger-degree polynomials. I never liked the arbitrary choice of a tolerance anyway, and now it seems to me that saying “evaluating the polynomial at the computed roots should be within this absolute distance from zero” is fundamentally the wrong thing to say; depending on the polynomial, we might have to take what we can get. Some other things I could imagine saying instead include:

• Evaluating the polynomial at the computed roots should be within some relative epsilon of zero, depending on the size/relative size of the coefficients
• The computed roots are as accurate as possible (or close to it) in the sense that evaluating the polynomial at other numbers near the computed roots yields values which are farther from zero

…but, first of all, I don’t know if these are reasonable properties to expect; and even if they were, I’m not sure I know how to express them in Haskell! Any advice is most welcome. Are there any best practices for expressing desirable test properties for floating-point computations?

For completeness, here is the actual code I came up with for finding roots via hmatrix. Notice there is another annoying arbitrary value in there, for deciding when a complex root is close enough to being real that we call it a real root. I’m not really sure what to do about this either.

-- Compute the roots of a polynomial as the eigenvalues of its companion matrix,
polyRoots :: [Double] -> [Double]
polyRoots []     = []
polyRoots (0:as) = polyRoots as
polyRoots (a:as) = mapMaybe toReal eigVals
  where
    n   = length as'
    as' = map (/a) as
    companion = (konst 0 (1,n-1) === ident (n-1)) ||| col (map negate . reverse $ as')

    eigVals = toList . fst . eig $ companion
    toReal (a :+ b)
       | abs b < 1e-10 = Just a   -- This arbitrary value is annoying too!
       | otherwise     = Nothing