Skip site navigation (1)Skip section navigation (2)
Date:      Mon, 30 Jul 2012 12:51:18 +1000
From:      Peter Jeremy <peter@rulingia.com>
To:        freebsd-numerics@freebsd.org
Subject:   Re: Bad error computing 1/z
Message-ID:  <20120730025118.GB7958@server.rulingia.com>
In-Reply-To: <5015DA8A.10200@missouri.edu>
References:  <5015DA8A.10200@missouri.edu>

next in thread | previous in thread | raw e-mail | index | archive | help

--CE+1k2dSO48ffgeK
Content-Type: text/plain; charset=us-ascii
Content-Disposition: inline
Content-Transfer-Encoding: quoted-printable

[Excuse the long lines]

On 2012-Jul-29 19:51:22 -0500, Stephen Montgomery-Smith <stephen@missouri.e=
du> wrote:
>As I was debugging catanh, I noticed the following oddness.
>
>If z =3D cpack(x,y), where
>x =3D 1
>y =3D 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=20
>is approx 4.2346e-309.  I know the reciprocal of this would cause an=20
>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+1=
54 0x5ff25691d4068c91
1.0/z: +4.2346036163332497e-309 0x00030b8596824af0   -6.5073832039716640e-1=
55 0x9febeb7fc100edd7

Doing the same thing with 40 digits precision in emacs-calc, I
calculate 1/z as:
(16#C.2E165A092BC21FF1FBEB5F45FAB8DA33*16.^-257, -16#D.F5BFE08076EBB470CE55=
E8BF55A1B0B5*16.^-129)
(4.234603616333252354574574265966349487777e-309, -6.50738320397166433428813=
7901223264158669e-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?

--=20
Peter Jeremy

--CE+1k2dSO48ffgeK
Content-Type: application/pgp-signature

-----BEGIN PGP SIGNATURE-----
Version: GnuPG v2.0.19 (FreeBSD)

iEYEARECAAYFAlAV9qYACgkQ/opHv/APuIcvPgCgqBX4lYN/wwnAlkEk1lzgUBbh
eJ8AoMD3ABR4zgi9TC2ew2WSLl3Gz4Pq
=nz4n
-----END PGP SIGNATURE-----

--CE+1k2dSO48ffgeK--



Want to link to this message? Use this URL: <https://mail-archive.FreeBSD.org/cgi/mid.cgi?20120730025118.GB7958>