Re: [boost] test_fp_comparisons and rounding errors

From: Guillaume Melquiond (gmelquio_at_[hidden])
Date: 2003-07-05 11:18:24

On Sat, 5 Jul 2003, Beman Dawes wrote:

> test_fp_comparisions has been failing for a long time. The issue has to do
> with how much tolerance to allow for rounding errors.
>
> The close_at_tolerance algorithm calculates tolerance as follows:
>
>     n*std::numeric_limits<T>::epsilon()/2
>
> where n is the number of possible floating-point rounding errors.
>
> The particular test case that is failing calls the comparison algorithm
> with a rounding error argument of 4. That allows for up to 2 epsilons
> tolerance.
>
> But if you step through the code, the actual error test values computed (d1
> and d2, in the code) is 3 epsilons for float, and 4 epsilons for double and
> long double.
>
> What is happening here? It seems to me that the error checking code itself
> that computes the values to be checked (d1 and d2) introduces one to four
> possible additional rounding errors. (The error checking code always does
> one subtraction, possibly one subtraction in an abs, possibly another
> subtraction in an abs, and possibly one division). That would account for
> the observed behavior.

Sorry, I don't see all the substractions you are speaking about, I only
see these three floating-point operations in close_at_tolerance:

  FPT diff = detail::fpt_abs( left - right );
  FPT d1   = detail::safe_fpt_division( diff, detail::fpt_abs( right ) );
  FPT d2   = detail::safe_fpt_division( diff, detail::fpt_abs( left ) );

We can suppose the user doesn't want to compare numbers with a relative
error bigger than 1/2 (they are no more "close" numbers). Consequently,
the substraction should not produce any rounding error (thanks to Sterbenz
law). However the divisions will produce rounding errors.

Consequently, the relative errors d1 and d2 are off by ulp/2.

> Is this analysis correct? I know almost nothing about floating-point
> arithmetic, so it is probably flawed. But it looks to me as if the correct
> tolerance formula is  (n+m)*std::numeric_limits<T>::epsilon()/2, where m is
> 1,2,3, or 4 depending on the exact logic path through

m should then be 1, unless I am mistaken.

> close_at_tolerance::operator(). It would also be easy to change the code to
> eliminate some of the operations which add additional possible rounding
> errors.

If the number of ulps is a power of 2 (and consequently tolerance is
2**(-k)), this code doesn't produce any rounding error:

  FPT diff = abs(a - b);
  bool c1 = diff <= tolerance * abs(a);
  bool c2 = diff <= tolerance * abs(b);

(if a and b are not too small)

> If the above is all wrong, perhaps someone can offer the correct analysis
> for why the tests fail.
>
> --Beman

There are also some strange things.

For example, what is the meaning of
  BOOST_CHECK_CLOSE_SHOULD_PASS_N
       (1, 1+std::numeric_limits<FPT>::epsilon() / 2, 1); ?
Indeed, 1+epsilon/2 is equal to 1 (by definition of epsilon) so this test
doesn't have any interest, does it?

Regards,

Guillaume

PS: I find the logging of test_fp_comparisons really ugly. I hope it is
not intended.