fix for wrong rounding of boost::math::round in classes of corner cases #8
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
This pull request fixes two classes of rounding bugs of
boost::math::round
occuring in two corner cases.The first corner case are the direct predecessor/successors of 0.5/-0.5. For
float
anddouble
types these are the following numbers:These numbers are nearer to zero than -0.5/0.5 and must be rounded to zero according to the documentation of
boost::math::round
, but are rounded to -1.0/1.0.To explain this, lets assume an imaginary floating-point type resembling IEEE 754 types with a four-bit mantissa (the exponent width is not important for this). As the positive and negative cases are symmetric, lets only look at the positive case. The significant expression of the original implementation to get the away-from-zero rounded number for
v
for this case is:0.5 is represented in binary as
0.10000b
(the1
is implicit and there are four0
bits in the mantissa). The direct predecessor (the largest representable number that is smaller than 0.5) is0.011111b
= 0.484375. The original implementation now adds 0.5 to this number:The result is too wide for the mantissa of the data type and is rounded by the precessor to
1.0000b
(actually, this is dependent on the rounding flags). Taking the floor from this number yields 1.0 as the result instead of the correct 0.0.The second corner case are numbers of the format
v = +/-1xyz1b
for arbitraryx
,y
andz
. In this case, a similar rounding problem occurs:Again, the last digit is not representable anymore and the processor can round this up
v + 1
(v - 1
for negativev
), yielding that as the result instead ofv
. This is especially problematic, as an already "round" (integral) number is rounded to another integral number.This pull request fixes the two corner cases by handling the the numbers
-0.5 < v < 0.5
as a special case to fix the first cornercase and by changing the calculation order for the second corner case. For positivev
, the significant expression now is:This calulcates
ceil(v)
first, avoiding any rounding asceil
doesn't need to round andv
is the raw input value. For all integralv
,ceil
yields exactlyv
, so the substraction result is zero without any rounding. For non-integralv
the substraction is alwaysexact and representable in the data type as the value of the result is always smaller thanv
. This guarantees that the comparison with 0.5 is correct. The substraction of 1 in the rounding-down raises no problem, as it is either exact or truncated if the precision ofv
can't represent 1 anymore.For negative
v
, the implementation/explanation is symmetric.The error probably won't occur on the x86 32 bit platform not using SSE, but using extended precision in the 387 FPU registers, I haven't tested this. I also didn't test the probable performance penalty of the correction (especially in a real-world application). Maybe there also exists a better fix than in this pull request (I didn't yet come up with one). As I think correctness (even in such admittedly (probably) rare corner cases on the edge of the precision) is more important than performance in most applications, a slightly slower implementation that is correct should be preferred.
I ran a (C++11 only) successfull test comparing
std::round
againstboost::math::round
for the complete finite value space offloat
. I did not include this, as it takes a long time to run, but is is in theroundbug-full-scan
branch of my fork of the boost.math repository.I did not look into non IEEE 754 floating point number representation, but I think that the implementation should not induce errors there.