Bad error computing 1/z

Stephen Montgomery-Smith stephen at
Mon Jul 30 03:05:14 UTC 2012

On 07/29/2012 09:51 PM, Peter Jeremy wrote:
> [Excuse the long lines]
> On 2012-Jul-29 19:51:22 -0500, Stephen Montgomery-Smith <stephen at> 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
>> overflow.
> 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:
> (16#C.2E165A092BC21FF1FBEB5F45FAB8DA33*16.^-257, -16#D.F5BFE08076EBB470CE55E8BF55A1B0B5*16.^-129)
> (4.234603616333252354574574265966349487777e-309, -6.507383203971664334288137901223264158669e-155)
> That give me fractions of:
> 0x0.30B8596824AF087FC7EFAD7D17EAE368D  0x1.BEB7FC100EDD768E19CABD17EAB43617
> Therefore the perfectly rounded double result should be
>    0x00030b8596824af1  0x9febeb7fc100edd7
> 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?

You answered my question.  The real part is denormal, and so my ULP 
calculation were misleading.

Thanks, Stephen

More information about the freebsd-numerics mailing list