From owner-freebsd-standards@FreeBSD.ORG Wed Dec 5 20:36:26 2007 Return-Path: Delivered-To: freebsd-standards@FreeBSD.ORG Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2001:4f8:fff6::34]) by hub.freebsd.org (Postfix) with ESMTP id 4444916A419 for ; Wed, 5 Dec 2007 20:36:26 +0000 (UTC) (envelope-from sgk@troutmask.apl.washington.edu) Received: from troutmask.apl.washington.edu (troutmask.apl.washington.edu [128.208.78.105]) by mx1.freebsd.org (Postfix) with ESMTP id 1B76713C465 for ; Wed, 5 Dec 2007 20:36:26 +0000 (UTC) (envelope-from sgk@troutmask.apl.washington.edu) Received: from troutmask.apl.washington.edu (localhost.apl.washington.edu [127.0.0.1]) by troutmask.apl.washington.edu (8.14.1/8.14.1) with ESMTP id lB5KWZRK044448; Wed, 5 Dec 2007 12:32:35 -0800 (PST) (envelope-from sgk@troutmask.apl.washington.edu) Received: (from sgk@localhost) by troutmask.apl.washington.edu (8.14.1/8.14.1/Submit) id lB5KWZ08044447; Wed, 5 Dec 2007 12:32:35 -0800 (PST) (envelope-from sgk) Date: Wed, 5 Dec 2007 12:32:35 -0800 From: Steve Kargl To: Bruce Evans , freebsd-standards@FreeBSD.ORG Message-ID: <20071205203235.GA40539@troutmask.apl.washington.edu> References: <20070928152227.GA39233@troutmask.apl.washington.edu> <20071001173736.U1985@besplex.bde.org> <20071002001154.GA3782@troutmask.apl.washington.edu> <20071002172317.GA95181@VARK.MIT.EDU> <20071002173237.GA12586@troutmask.apl.washington.edu> <20071003103519.X14175@delplex.bde.org> <20071010204249.GA7446@troutmask.apl.washington.edu> <20071203074407.GA10989@VARK.MIT.EDU> <20071203145105.GA16203@troutmask.apl.washington.edu> <20071205185132.GA91591@VARK.MIT.EDU> Mime-Version: 1.0 Content-Type: text/plain; charset=us-ascii Content-Disposition: inline In-Reply-To: <20071205185132.GA91591@VARK.MIT.EDU> User-Agent: Mutt/1.4.2.3i Cc: Subject: Re: long double broken on i386? X-BeenThere: freebsd-standards@freebsd.org X-Mailman-Version: 2.1.5 Precedence: list List-Id: Standards compliance List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Wed, 05 Dec 2007 20:36:26 -0000 On Wed, Dec 05, 2007 at 01:51:32PM -0500, David Schultz wrote: > On Mon, Dec 03, 2007, Steve Kargl wrote: > > On Mon, Dec 03, 2007 at 02:44:08AM -0500, David Schultz wrote: > > > Is the latest version of the patch the one you sent to the list? > > > If cosl and friends produce accurate answers within a reasonable > > > domain, it's worth committing; the whole business about how > > > accurate cosl(1000000000) is can be dealt with later. > > > > No, I've rewritten most of the code (numerous times in fact :-). > > > > First, I followed Bruce's recommendation to compute the > > ULP in the higher precision available in GMP/MPFR instead > > of the long double type. The old method appears to > > have been using a chop rounding because I would get > > only 0.5, 1.0, 1.5, ... ULPs (ie, I never saw anything like > > 0.6312 ULP). With the new method, I see ULPs in the 0.6 > > to 0.75 range. I haven't had time to implement a full > > blown Remes algorithm to see if I can reduce the ULP. > > > > Second, I've implemented the algorithms found in > > s_cos.c, s_sin.c, and s_tan.c that use the extra > > precision from the argument reduction. If I use > > these algorithms, then the ULPs go up! For example, > > the ULP for sinl() goes up to about 1.5. Bruce > > had suggested that the 396 hex digits in the table > > for 2/pi may need to be checked. I found a paper > > by K.C. Ng at David Hough's website that (if I read it > > correctly) suggests that we need a much larger > > table. So, I need to check how well argument reduction > > is working. > > Having a version that works for machines with the 128-bit floating > point format is pretty important. (These would probably be two > separate files; I've been thinking we should add a msun/ieee96 > directory and a msun/ieee128 directory or something.) This is probably a good idea if freebsd wants to maintain the same algorithm across, say, sinf, sin, and sinl. I can produce code for the ieee128 case, but I have no way to test it. An alternative may be to use table-driven implementations where the intervals are defined by exactly representable values in ieee96 (i.e., 64-bit significand), which are then also exactly representable in ieee128. I haven't investigated how many intervals one would need nor the best interpolation scheme. > But honestly, I've tried to wrestle with the argument reduction > stuff before, and my advice is to not kill yourself over it. (large arg discussion elided) > That's not to say that we need no argument reduction at all. For > instance, cosl(5*pi/3) should still give an accurate answer. But > when the input is so huge that the exponent is bigger than the > number of significant bits in the mantissa, you lose so much > precision in the input that it's not as important anymore. That's > probably why Intel decided to make its trig instructions only work > up to 2^63 before requiring explicit argument reduction. This is what I need to check. For reasonable args, does the arg reduction provide extra precision and is it useful? It seems that you may have already checked this. I hope to revisit some of this on Saturday. -- Steve