From owner-freebsd-numerics@FreeBSD.ORG  Mon Jul 30 00:51:23 2012
Return-Path: <owner-freebsd-numerics@FreeBSD.ORG>
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 50CC51065670
	for <freebsd-numerics@freebsd.org>;
	Mon, 30 Jul 2012 00:51:23 +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 010658FC08
	for <freebsd-numerics@freebsd.org>;
	Mon, 30 Jul 2012 00:51:22 +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
	q6U0pMVb053462 for <freebsd-numerics@freebsd.org>;
	Sun, 29 Jul 2012 19:51:22 -0500 (CDT)
	(envelope-from stephen@missouri.edu)
Message-ID: <5015DA8A.10200@missouri.edu>
Date: Sun, 29 Jul 2012 19:51:22 -0500
From: Stephen Montgomery-Smith <stephen@missouri.edu>
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
Content-Type: text/plain; charset=ISO-8859-1; format=flowed
Content-Transfer-Encoding: 7bit
Subject: 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."
	<freebsd-numerics.freebsd.org>
List-Unsubscribe: <http://lists.freebsd.org/mailman/listinfo/freebsd-numerics>, 
	<mailto:freebsd-numerics-request@freebsd.org?subject=unsubscribe>
List-Archive: <http://lists.freebsd.org/pipermail/freebsd-numerics>
List-Post: <mailto:freebsd-numerics@freebsd.org>
List-Help: <mailto:freebsd-numerics-request@freebsd.org?subject=help>
List-Subscribe: <http://lists.freebsd.org/mailman/listinfo/freebsd-numerics>, 
	<mailto:freebsd-numerics-request@freebsd.org?subject=subscribe>
X-List-Received-Date: Mon, 30 Jul 2012 00:51:23 -0000

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 tried writing my own reciprocal algorithm, and it had exactly the same 
error.

Does anyone have an explanation?

I thought it might be that the number was somehow "denormalized" and it 
couldn't store enough digits in the mantissa.  Or maybe the error is 
when mpfr converts the correct answer to double.  I don't know.