# Bad error computing 1/z

Stephen Montgomery-Smith stephen at missouri.edu
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 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
>> 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:
> 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