This is the mail archive of the
gcc@gcc.gnu.org
mailing list for the GCC project.
Re: Using crlibm as the default math library in GCC sources
On 2007-11-16 19:33:38 -0500, Geert Bosch wrote:
> If a real number is rounded to the nearest number in a binary floating
> point format with 2*p+1 bits and that result rounded again to
> a p-bit float, the result is the same as if the original value was
> rounded to a p-bit float directly.
No, this is wrong: if you round a number x to 2p+1 bits and get
a number x' which is a rounding boundary for p-bit floats (i.e.
a p+1-bit number whose bit p+1 is 1), then x' will always be
rounded to some p-bit number x'', whether x was above or below x'.
For instance, consider
log(1.1000101001001101000000010101000110010110101100010101 * 2^37).
The significand y of the result is:
1.1010000101000001000011011100111011111110111101101101 0 1^56 0100...
\_____________________ 53 bits ______________________/
where 1^56 means: a sequence of 56 bits 1. One has:
rnd_107(y) =
1.1010000101000001000011011100111011111110111101101101 1
rnd_53(rnd_107(y)) =
1.1010000101000001000011011100111011111110111101101110
due to the even-rounding rule.
But rnd_53(y) =
1.1010000101000001000011011100111011111110111101101101
Of course, it is possible to prove that such a case doesn't happen
for some simple functions (perhaps the square root? I also have a
similar result for floor(x/y) [*]), but in general, no proofs are
known.
[*] http://www.vinc17.org/research/papers/rr_intdiv.pdf Section 5.2.
Alternatively, you need the ternary value for the first rounding (what
MPFR gives, but CRlibm doesn't give it).
--
Vincent Lefèvre <vincent@vinc17.org> - Web: <http://www.vinc17.org/>
100% accessible validated (X)HTML - Blog: <http://www.vinc17.org/blog/>
Work: CR INRIA - computer arithmetic / Arenaire project (LIP, ENS-Lyon)