From owner-cvs-src@FreeBSD.ORG Fri Mar 7 00:18:27 2008 Return-Path: Delivered-To: cvs-src@freebsd.org Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2001:4f8:fff6::34]) by hub.freebsd.org (Postfix) with ESMTP id A4E411065671 for ; Fri, 7 Mar 2008 00:18:27 +0000 (UTC) (envelope-from cperciva@freebsd.org) Received: from pd3mo1so.prod.shaw.ca (idcmail-mo1so.shaw.ca [24.71.223.10]) by mx1.freebsd.org (Postfix) with ESMTP id 799038FC13 for ; Fri, 7 Mar 2008 00:18:27 +0000 (UTC) (envelope-from cperciva@freebsd.org) Received: from pd3mr6so.prod.shaw.ca (pd3mr6so-qfe3.prod.shaw.ca [10.0.141.21]) by l-daemon (Sun ONE Messaging Server 6.0 HotFix 1.01 (built Mar 15 2004)) with ESMTP id <0JXC00M393HL8ME0@l-daemon> for cvs-src@freebsd.org; Thu, 06 Mar 2008 17:17:45 -0700 (MST) Received: from pn2ml7so.prod.shaw.ca ([10.0.121.151]) by pd3mr6so.prod.shaw.ca (Sun Java System Messaging Server 6.2-7.05 (built Sep 5 2006)) with ESMTP id <0JXC006J53HJVN70@pd3mr6so.prod.shaw.ca> for cvs-src@freebsd.org; Thu, 06 Mar 2008 17:17:45 -0700 (MST) Received: from hexahedron.daemonology.net ([24.80.10.198]) by l-daemon (Sun ONE Messaging Server 6.0 HotFix 1.01 (built Mar 15 2004)) with SMTP id <0JXC00M4X3HHDT60@l-daemon> for cvs-src@freebsd.org; Thu, 06 Mar 2008 17:17:43 -0700 (MST) Received: (qmail 1975 invoked from network); Fri, 07 Mar 2008 00:17:40 +0000 Received: from unknown (HELO hexahedron.daemonology.net) (127.0.0.1) by localhost with SMTP; Fri, 07 Mar 2008 00:17:40 +0000 Date: Thu, 06 Mar 2008 16:17:40 -0800 From: Colin Percival In-reply-to: <20080307085031.P10724@delplex.bde.org> To: Bruce Evans Message-id: <47D089A4.4060208@freebsd.org> MIME-version: 1.0 Content-type: text/plain; charset=windows-1252 Content-transfer-encoding: 8BIT X-Enigmail-Version: 0.95.5 References: <200803051121.m25BLE03035426@repoman.freebsd.org> <20080305182531.GS68971@server.vk2pj.dyndns.org> <20080306021222.GA46783@zim.MIT.EDU> <47CF5D19.3090100@freebsd.org> <20080306033246.GA47280@zim.MIT.EDU> <47CF7EBF.6000009@freebsd.org> <20080306063452.GB48339@zim.MIT.EDU> <47CF9586.70707@freebsd.org> <20080307085031.P10724@delplex.bde.org> User-Agent: Thunderbird 2.0.0.9 (X11/20071117) Cc: cvs-src@freebsd.org, src-committers@freebsd.org, cvs-all@freebsd.org Subject: Re: cvs commit: src/sys/i386/include _types.h X-BeenThere: cvs-src@freebsd.org X-Mailman-Version: 2.1.5 Precedence: list List-Id: CVS commit messages for the src tree List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Fri, 07 Mar 2008 00:18:27 -0000 Bruce Evans wrote: > On Wed, 5 Mar 2008, Colin Percival wrote: >> David Schultz wrote: >>> On Wed, Mar 05, 2008, Colin Percival wrote: >>>> Yes and no. Double rounding isn't allowed; >>> >>> Yes, it is. >> >> No, it isn't. The following code > > Yes it is :-). Extra precision for intermediate values implies double > rounding for the final results, and extra precision for intermediate > values is a fundamental part of C (it used to be only that float > expressions were (de-facto standard required to be) evaluated in double > precision, but now it affects double expressions too, and C99 has vast > complications to support i387'-ish extra precision. Let me clarify: Double rounding isn't allowed by IEEE754. To quote from the October 5th draft of IEEE 754R (which I happen to have in front of me right now): 5.1 Overview All conforming implementations of this standard shall provide the operations listed in this clause for all supported floating-point formats available for arithmetic. Each of the computational operations that return a numeric result specified by this standard shall be performed as if it first produced an intermediate result correct to infinite precision and with unbounded range, ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ and then rounded that intermediate result to fit in the destination’s format ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ This is clarified further in clause 10: 10.1 Expression evaluation rules Clause 5 of this standard specifies the result of a single arithmetic operation going to an explicit destination. Every operation has an explicit or implicit destination. One rounding occurs to fit the exact result into a ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ destination format. That result is reproducible in that the same operation ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ applied to the same operands under the same attributes produces the same ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ result on all conforming implementations in all languages. ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ Programming language standards might define syntax for expressions that combine one or more operations of this standard, producing a result to fit an explicit or implicit final destination. When a variable with a declared format ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ is a final destination, as in format conversion to a variable, that declared ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ format of that variable governs its rounding. The format of an implicit ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ destination, or of an explicit destination without a declared format, is defined by language standard expression evaluation rules. In other words, in the code double x, y, z; z = x * y + 1.0; it is perfectly legal for (x * y) to be computed with extra precision; but in the code double x, y, z; z = x * y; z = z - 1.0; both the multiplication and the subtraction must be rounded *once* to double precision. > In C, with FLT_EVAL_METHOD = 2 (evaluate all FP expressions in the range > and precision of the long double type) [...] As someone wrote earlier, the authors of C99 were a bit confused. I can only assume that they intended to be IEEE754-compliant and the rules concerning FLT_EVAL_METHOD apply to *implicit* destinations only. >>> It's impossible to implement efficient library functions that >>> produce correct long double results >> >> You could have stopped this sentence here -- for all practical purposes, >> correctly rounded trigonometric functions are not feasible. > > Nah, it's quite feasible, especially for transcendental functions since > transcendental numbers in general can be approximated more closely by > rationals than algebraic numbers. In fact, one of the common ways to prove that a real number is transcendental is to show that it can be approximated more closely by rationals than any algebraic function can be. > The evaluation just needs to be > done in extra precision How much extra precision do you need? If the exact value of a transcendental function is extremely close to halfway between two consecutive floating-point values, how do you decide which way to round it? Computing transcendental functions to within (0.5+x)ulp for small positive values of x is certainly feasible, but computing provably correctly rounded transcendental functions is very very hard. >> Fortunately, library functions aren't required to have any particular >> error bounds. > > We require error bounds of < 1 ulp and aim for < 0.5[0]1 ulps. What standard states that those bounds are required (or can be relied upon when proving that code is correct)? Intel FPUs over the years have been variously documented to compute transcendental functions to within 1.0, 2.5, or 3.5 ulp -- are you saying that the 8087 doesn't comply with the standard which it inspired? Colin Percival