From owner-freebsd-numerics@FreeBSD.ORG Mon Jul 30 02:51:30 2012 Return-Path: Delivered-To: freebsd-numerics@freebsd.org Received: from mx1.freebsd.org (mx1.freebsd.org [69.147.83.52]) by hub.freebsd.org (Postfix) with ESMTP id 956991065670 for ; Mon, 30 Jul 2012 02:51:30 +0000 (UTC) (envelope-from peter@rulingia.com) Received: from vps.rulingia.com (host-122-100-2-194.octopus.com.au [122.100.2.194]) by mx1.freebsd.org (Postfix) with ESMTP id 0B3CD8FC0A for ; Mon, 30 Jul 2012 02:51:29 +0000 (UTC) Received: from server.rulingia.com (c220-239-247-45.belrs5.nsw.optusnet.com.au [220.239.247.45]) by vps.rulingia.com (8.14.5/8.14.5) with ESMTP id q6U2pQ7Y081126 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=OK) for ; Mon, 30 Jul 2012 12:51:27 +1000 (EST) (envelope-from peter@rulingia.com) X-Bogosity: Ham, spamicity=0.000000 Received: from server.rulingia.com (localhost.rulingia.com [127.0.0.1]) by server.rulingia.com (8.14.5/8.14.5) with ESMTP id q6U2pJqC082836 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=NO) for ; Mon, 30 Jul 2012 12:51:19 +1000 (EST) (envelope-from peter@server.rulingia.com) Received: (from peter@localhost) by server.rulingia.com (8.14.5/8.14.5/Submit) id q6U2pI2D082835 for freebsd-numerics@freebsd.org; Mon, 30 Jul 2012 12:51:18 +1000 (EST) (envelope-from peter) Date: Mon, 30 Jul 2012 12:51:18 +1000 From: Peter Jeremy To: freebsd-numerics@freebsd.org Message-ID: <20120730025118.GB7958@server.rulingia.com> References: <5015DA8A.10200@missouri.edu> MIME-Version: 1.0 Content-Type: multipart/signed; micalg=pgp-sha1; protocol="application/pgp-signature"; boundary="CE+1k2dSO48ffgeK" Content-Disposition: inline In-Reply-To: <5015DA8A.10200@missouri.edu> X-PGP-Key: http://www.rulingia.com/keys/peter.pgp User-Agent: Mutt/1.5.21 (2010-09-15) Subject: Re: Bad error computing 1/z X-BeenThere: freebsd-numerics@freebsd.org X-Mailman-Version: 2.1.5 Precedence: list List-Id: "Discussions of high quality implementation of libm functions." List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Mon, 30 Jul 2012 02:51:30 -0000 --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 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--