From owner-freebsd-numerics@FreeBSD.ORG Mon Jul 30 03:05:14 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 605E6106566B for ; Mon, 30 Jul 2012 03:05:14 +0000 (UTC) (envelope-from stephen@missouri.edu) Received: from wilberforce.math.missouri.edu (wilberforce.math.missouri.edu [128.206.184.213]) by mx1.freebsd.org (Postfix) with ESMTP id 274E98FC16 for ; Mon, 30 Jul 2012 03:05:14 +0000 (UTC) Received: from [127.0.0.1] (wilberforce.math.missouri.edu [128.206.184.213]) by wilberforce.math.missouri.edu (8.14.5/8.14.5) with ESMTP id q6U35D3X062352 for ; Sun, 29 Jul 2012 22:05:13 -0500 (CDT) (envelope-from stephen@missouri.edu) Message-ID: <5015F9E9.9070408@missouri.edu> Date: Sun, 29 Jul 2012 22:05:13 -0500 From: Stephen Montgomery-Smith User-Agent: Mozilla/5.0 (X11; Linux i686; rv:14.0) Gecko/20120714 Thunderbird/14.0 MIME-Version: 1.0 To: freebsd-numerics@freebsd.org References: <5015DA8A.10200@missouri.edu> <20120730025118.GB7958@server.rulingia.com> In-Reply-To: <20120730025118.GB7958@server.rulingia.com> Content-Type: text/plain; charset=ISO-8859-1; format=flowed Content-Transfer-Encoding: 7bit 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 03:05:14 -0000 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 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