[RFC] new function for floating point comparison
Hi,
I have added a couple of utilities for floating point comparison, to be used in unit tests mostly, and would like some comments, especially from people knowledgeable about floating point.
http://github.com/cournape/numpy/tree/new_ulp_comp
The main difference compared to other functions is that they are 'amplitudeindependent', and use IEEE754specific properties. The tolerance is based on ULP, and two numbers x, y are closed depending on how many numbers are representable between x and y at the given precision. The branch contains the following new functions:
* spacing(x): equivalent to the F90 intrinsic. Returns the smallest representable number needed so that spacing(x) + x > x. Spacing(1) is EPS by definition. * assert_array_almost_equal_nulp(x, y, nulp=1): assertion is defined as abs(x  y) <= nulps * spacing(max(abs(x), abs(y))). * assert_array_max_ulp(a, b, maxulp=1, dtype=None): given two numbers a and b, raise an assertion if there are more than maxulp representable numbers between a and b.
They only support single and double precision  for complex number, one could arbitrarily define a distance between numbers based on nulps, say max of number of representable number for real and imag parts. Extended precision would be a bit more painful, because of the variety of implementations.
I hope that they can give more robust/meaningful comparison for most of our unit tests,
cheers,
David
On Thu, Oct 29, 2009 at 02:17, David Cournapeau david@ar.media.kyotou.ac.jp wrote:
Hi,
I have added a couple of utilities for floating point comparison, to be used in unit tests mostly, and would like some comments, especially from people knowledgeable about floating point.
http://github.com/cournape/numpy/tree/new_ulp_comp
The main difference compared to other functions is that they are 'amplitudeindependent', and use IEEE754specific properties. The tolerance is based on ULP, and two numbers x, y are closed depending on how many numbers are representable between x and y at the given precision. The branch contains the following new functions:
* spacing(x): equivalent to the F90 intrinsic. Returns the smallest representable number needed so that spacing(x) + x > x. Spacing(1) is EPS by definition. * assert_array_almost_equal_nulp(x, y, nulp=1): assertion is defined as abs(x  y) <= nulps * spacing(max(abs(x), abs(y))). * assert_array_max_ulp(a, b, maxulp=1, dtype=None): given two numbers a and b, raise an assertion if there are more than maxulp representable numbers between a and b.
That sounds good. Another worthwhile addition would be nextafter().
http://www.opengroup.org/onlinepubs/000095399/functions/nextafter.html
With a little bit of care, a nextafter ufunc can be used to generate a dense grid of floating point values around a given center. This can be used to explore the error characteristics of a function at a very fine level of detail that is otherwise unavailable.
Robert Kern wrote:
That sounds good. Another worthwhile addition would be nextafter().
http://www.opengroup.org/onlinepubs/000095399/functions/nextafter.html
Ah, I did not know about this one. I have implemented it and committed it. One issue is that it will cause failures on platforms without nextafterl and where long double != double, but I don't think we have a lot of those, if any,
David
participants (2)

David Cournapeau

Robert Kern