Bad error computing 1/z
peter at rulingia.com
Mon Jul 30 02:51:30 UTC 2012
[Excuse the long lines]
On 2012-Jul-29 19:51:22 -0500, Stephen Montgomery-Smith <stephen at missouri.edu> wrote:
>As I was debugging catanh, I noticed the following oddness.
>If z = cpack(x,y), where
>x = 1
>y = 0x1.25691d4068c910p+512, approximately 1.53672e+154
>then the real part of 1/z is wrong by about 4ULP. The real part of 1/z
>is approx 4.2346e-309. I know the reciprocal of this would cause an
I think I'm missing a bit of context here. If I calculate it using
doubles with gcc, I get:
z: +1.0000000000000000e+00 0x3ff0000000000000 +1.5367160172612364e+154 0x5ff25691d4068c91
1.0/z: +4.2346036163332497e-309 0x00030b8596824af0 -6.5073832039716640e-155 0x9febeb7fc100edd7
Doing the same thing with 40 digits precision in emacs-calc, I
calculate 1/z as:
That give me fractions of:
Therefore the perfectly rounded double result should be
ie, the real part calculates as 1ULP low when using doubles.
Note that since the real part is a denormal, it only has 50 bits of
precision. Are you sure your ULP calculations are taking into account
the number of fraction bits available when numbers are denormals?
-------------- next part --------------
A non-text attachment was scrubbed...
Name: not available
Size: 196 bytes
Desc: not available
Url : http://lists.freebsd.org/pipermail/freebsd-numerics/attachments/20120730/b4929679/attachment.pgp
More information about the freebsd-numerics