From owner-freebsd-standards@FreeBSD.ORG Mon Oct 1 09:56:38 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 9020D16A469 for ; Mon, 1 Oct 2007 09:56:38 +0000 (UTC) (envelope-from brde@optusnet.com.au) Received: from mail15.syd.optusnet.com.au (mail15.syd.optusnet.com.au [211.29.132.196]) by mx1.freebsd.org (Postfix) with ESMTP id 0AD2B13C45A for ; Mon, 1 Oct 2007 09:56:37 +0000 (UTC) (envelope-from brde@optusnet.com.au) Received: from besplex.bde.org (c220-239-235-248.carlnfd3.nsw.optusnet.com.au [220.239.235.248]) by mail15.syd.optusnet.com.au (8.13.1/8.13.1) with ESMTP id l919uS9W023309 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=NO); Mon, 1 Oct 2007 19:56:34 +1000 Date: Mon, 1 Oct 2007 19:56:28 +1000 (EST) From: Bruce Evans X-X-Sender: bde@besplex.bde.org To: Steve Kargl In-Reply-To: <20070928152227.GA39233@troutmask.apl.washington.edu> Message-ID: <20071001173736.U1985@besplex.bde.org> References: <20070928152227.GA39233@troutmask.apl.washington.edu> MIME-Version: 1.0 Content-Type: TEXT/PLAIN; charset=US-ASCII; format=flowed Cc: freebsd-standards@freebsd.org 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: Mon, 01 Oct 2007 09:56:38 -0000 On Fri, 28 Sep 2007, Steve Kargl wrote: > So, in my never ending crusade to implement missing > C99 math functions, I decided to tackle sinl, cosl, > and tanl. I have implementations that use argument > reduction from k_rem_pio2.c and Chebyshev interpolation > for the trig functions over the range [0,pi/4]. > > On amd64, these routines give <= 1/2 ULP for the various > values I've tested in the reduced range. Now, if I test I doubt this :-). Getting <= 1 ULP is hard; getting <= 1/2 + 0.001 ULPs is harder, and getting <= 1/2 ULPs without using an enormous amount of extra precision is a research problem. > these on i386, I see between 400 and 1200 ULP. I spent > a week or so trying to understand the descrepancy and "fix" > the problem using a double-double precision arithemetic. > > I finally found the problem! /usr/include/float.h on i386 > is lying about its numerical representation of long double. This problem is well known. I must have referred to it many times in old mail to you. FreeBSD uses my old workaround for compiler bugs in double precision, of using a rounding precision of double for everything. Unfortunately, this is still necessary because the compiler bugs are still there. This mostly breaks long double precision, and it doesn't fix float precision, but double precision is more important. libm, in particular k_rem_pio2f.c (which is no longer used) needs tricky fixes to work around the bugs for float precision. The fixes have mostly not been merged to the double precision functions -- we just depend on double precision mostly working. > In particular, at least these 3 values appear to be wrong > > #define LDBL_MANT_DIG 64 > #define LDBL_EPSILON 1.0842021724855044340E-19L > #define LDBL_DIG 18 These were broken in rev.1.9 in 2002. The "slight rounding issues with long doubles" referred to in the log message are that long doubles are rounded to double precision when they are operated on, unless you change the rounding precision to long double using fpsetprec() or fesetprec(). Even if you change the rounding precision, long doubles can only appear to work, but not actually work, unless you fix the compiler bugs or work around them in libraries and in your application. It isn't clear if the "appear to work" in the log messages is with or without the change to the rounding precision. As mentioned in the log message, printf() in 2002 didn't support printing long doubles at all, so it was hard to see loss of precision in printf output. Now printf supports long doubles and %a, so it is easy to see lost precision using it. E.g.: %%% #include #include int main(void) { printf(" sin(1) = %La\n", (long double)(double)sin(1)); printf("2 * sin(1) = %La\n", (long double)(double)(2 * sin(1))); return (0); } %%% Output on i386: %%% sin(1) = 0xd.76aa47848677021p-4 2 * sin(1) = 0xd.76aa47848677p-3 %%% This exhibits a library bug, several compiler bugs, and the lost precision from my rounding precision hack: o library bug: sin() returns extra precision, despite its prototype saying that it returns double. Well, I used to think that this is a bug. I only leared a few months ago that it is a feature, since it matches the corresponding feature for a function implemented in C. On i386, sin() is implemented in asm. It just uses the hardware to evaluate sin() in long double precision, and returns the long double result of that. The corresponding hardware feature is that in a C function that does similar calculations is permitted to return extra precision: double x, y; long double foo() { return x+y; } This is permitted to do the following: - evaluate x+y in extra precision - return the extra precision, since unlike casts and assignments, the return operation is _not_ required to discard any extra precision. For gcc on i386, this is what it would do without my rounding precision hack, except possibly when compiled with -O0 and/or if the result of x+y is spilled to memory (-O0 makes spilling quite likely). The return statement must be written as "return (double)(x+y);" to ensure removing any extra precision _before_ returning. Now I wonder if the implicit casts for function args with prototyped functions are permitted to be as magic as the non-casts for return statments. Keeping extra precision in intermediate calculations is good, and you don't really want to lose it unconditionally at function call boundaries, and it's confusing to lose it on calls but not on returns. I'm using sin(1) mainly to generate a long double with nonzero in its 11 least significant mantissa bits. Note that initializing a long double to apparently-the-same value at compile time would not work: long double x = 0xd.76aa47848677021p-4L; /* loses extra precision */ This is because gcc on FreeBSD supports my precision hack to a fault, and has bugs in the fault. gcc now tries to do compile-time arithmetic in the same precision as the runtime arithmetic, i.e., in only double precision on FreeBSD. This is wrong if there is no arithmetic involved, as in the above initializer. The runtime rounding precision only applies to operations -- it doesn't prevent loading, storing or non-arithmetic (function) operations on long doubles with full precision. The above initializer corresponds to a store of a value that has somehow been calculated in full precision. Note that this bug makes it more difficult to implement functions in long double precision -- you cannot use tables of long double constants. OTOH, loads of long double constants are so slow on i386 that such constants should almost never be used. Instead represent long doubles as the sum of 2 doubles if possible. Loading 2 doubles and adding them is faster than loading 1 long double, at least if the 2 extra operations can be pipelined effectively. o compiler bug (?). Casts are required to remove any extra precision. I cast sin(1) to double to do this. gcc elides this cast since it thinks that sin(1) has no extra precision since it has type double, and it knows that sin(1) is in a long double register, so it thinks that nothing needs to be done for the combined cast (long double)(double). However, having type double doesn't imply no extra precision in general, and after a return statement is one class of cases where it never does. A compiler can only tell if there is extra precision by looking inside the function near its return statement, or perhaps using an attribute. o compiler bug. The casts in "(long double)(double)(2 * sin(1)))" are elided similarly. This is not visible in the printf output, since my rounding precision hack actually makes a difference here: ... o ... lost precision. The hack results in the multiplication not producing any extra precision, so the printf output is correct despite the compiler bug. Variations on the above: - with 1*sin(1) instead of 2*sin(1), optimizing away the multiplication is valid and done with -O, so the extra precision is retained. - with sin(1) == sin(1) instead of 2*sin(1), and %d for the format specifier, the precision bugs combined with spilling bugs result in the value 0. Similarly, sin(x) != sin(x) for almost all x. This behaviour is probably permitted, but it is surprising, and it isn't intentional. C99 apparently permits sin(x) to be evaluated _and_ returned in extra precision, and I don't know of any requirement that the return value always has the same amount of extra precision, and there are no casts in sight to clip the precision. Different amounts always result in practice, accidentally as follows: gcc doesn't really understand exztra precision, especially for spilling intermediate values. It just spills one of the sin(1)'s and thus loses the extra precision for only one of them. Loss of precision from spilling is normal near function calls because of the ABI, but uncommon (but much harder to predict) otherwise. > #include > > int > main (void) > { > int i; > long double a; > a = 1.L; > for (i = 1; i < 64; i++) > { > a /= 2; > if (1.L - a == 1.L) break; > } > printf("%d %Le\n", i, a); > return 0; > } > > amd64:sgk[206] ./z > 64 1.084202e-19 > > i386:kargl[205] ./z > 54 5.551115e-17 > > The above results for i386 are consistent with DBL_MANT_DIG > and DBL_EPSILON. Is this intentional, and should float.h be > fixed? LBDL_MANT_DIG, LDBL_EPSILON and LDBL_EPSILON were broken in 2002. Note that the compiler bugs make use of the epsilons fragile. The delicate casts or assignments needed to give the lower-precision epsilons a chance of working don't actually work unless you use hacks like -ffloat store and volatile variables. LDBL_MIN_* and LDBL_MAX_* were fixed in 2002. My precision hack only affects precision (only of +-*/ operations). Long doubles still have extra exponent range, so LDBL_MAX != DBL_MAX, etc. This was wrong before 2002. Bruce From owner-freebsd-standards@FreeBSD.ORG Mon Oct 1 11:08:45 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 12CF316A50A for ; Mon, 1 Oct 2007 11:08:45 +0000 (UTC) (envelope-from owner-bugmaster@FreeBSD.org) Received: from freefall.freebsd.org (freefall.freebsd.org [IPv6:2001:4f8:fff6::28]) by mx1.freebsd.org (Postfix) with ESMTP id EA4C813C4BE for ; Mon, 1 Oct 2007 11:08:44 +0000 (UTC) (envelope-from owner-bugmaster@FreeBSD.org) Received: from freefall.freebsd.org (gnats@localhost [127.0.0.1]) by freefall.freebsd.org (8.14.1/8.14.1) with ESMTP id l91B8iS4064631 for ; Mon, 1 Oct 2007 11:08:44 GMT (envelope-from owner-bugmaster@FreeBSD.org) Received: (from gnats@localhost) by freefall.freebsd.org (8.14.1/8.14.1/Submit) id l91B8hjI064627 for freebsd-standards@FreeBSD.org; Mon, 1 Oct 2007 11:08:43 GMT (envelope-from owner-bugmaster@FreeBSD.org) Date: Mon, 1 Oct 2007 11:08:43 GMT Message-Id: <200710011108.l91B8hjI064627@freefall.freebsd.org> X-Authentication-Warning: freefall.freebsd.org: gnats set sender to owner-bugmaster@FreeBSD.org using -f From: FreeBSD bugmaster To: freebsd-standards@FreeBSD.org Cc: Subject: Current problem reports assigned to you 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: Mon, 01 Oct 2007 11:08:45 -0000 Current FreeBSD problem reports Critical problems Serious problems S Tracker Resp. Description -------------------------------------------------------------------------------- o bin/25542 standards /bin/sh: null char in quoted string o kern/46239 standards posix semaphore implementation errors o stand/54410 standards one-true-awk not POSIX compliant (no extended REs) o stand/82654 standards C99 long double math functions are missing o stand/94729 standards [libc] fcntl() throws undocumented ENOTTY f stand/114910 standards getaddrinfo() fails to set ai_canonname 6 problems total. Non-critical problems S Tracker Resp. Description -------------------------------------------------------------------------------- o bin/21519 standards sys/dir.h should be deprecated some more o bin/24390 standards Replacing old dir-symlinks when using /bin/ln s stand/24590 standards timezone function not compatible witn Single Unix Spec s kern/28260 standards UIO_MAXIOV needs to be made public s stand/36076 standards Implementation of POSIX fuser command o stand/39256 standards snprintf/vsnprintf aren't POSIX-conformant for strings p stand/41576 standards POSIX compliance of ln(1) o stand/44425 standards getcwd() succeeds even if current dir has perm 000. o stand/46119 standards Priority problems for SCHED_OTHER using pthreads o stand/54833 standards [pcvt] more pcvt deficits o stand/54839 standards [pcvt] pcvt deficits p stand/55112 standards glob.h, glob_t's gl_pathc should be "size_t", not "int o stand/56476 standards cd9660 unicode support simple hack o stand/58676 standards grantpt(3) alters storage used by ptsname(3) s stand/62858 standards malloc(0) not C99 compliant s kern/64875 standards [libc] [patch] [feature request] add a system call: fd o stand/66357 standards make POSIX conformance problem ('sh -e' & '+' command- o stand/66531 standards _gettemp uses a far smaller set of filenames than docu o stand/70813 standards [PATCH] ls(1) not Posix compliant o stand/72006 standards floating point formating in non-C locales o stand/79056 standards regex(3) regression tests a stand/80293 standards sysconf() does not support well-defined unistd values o stand/81287 standards [PATCH]: fingerd(8) might send a line not ending in CR o stand/83845 standards [libm] [patch] add log2() and log2f() support for libm o stand/85080 standards output of long double subnormals (with printf) is wron o stand/92360 standards [headers] [patch] Missing TAB3 in kernel headers o stand/92362 standards [headers] [patch] Missing SIGPOLL in kernel headers o kern/93705 standards [headers] [patch] ENODATA and EGREGIOUS (for glibc com o stand/96016 standards [headers] clock_getres et al should be in o stand/96236 standards [PATCH] [POSIX] sed.1 incorrectly describes a function p stand/99517 standards Missing SIGRTMIN and SIGRTMAX signals o stand/99960 standards [Patch] make(1): Add -p flag o stand/100017 standards [Patch] Add fuser(1) functionality to fstat(1) o stand/104743 standards [headers] [patch] Wrong values for _POSIX_ minimal lim o stand/104841 standards [libm] [patch] C99 long double square root. o stand/107561 standards [patch] Missing SUS function tcgetsid() o kern/114578 standards [libc] wide character printing using swprintf(dst, n, o stand/114633 standards /etc/rc.subr: line 511: omits a quotation mark: "force o stand/116081 standards make does not work with the directive sinclude o stand/116221 standards SUS issue -- FreeBSD has not flag WNOWAIT for wait*() 40 problems total. From owner-freebsd-standards@FreeBSD.ORG Tue Oct 2 00:14:59 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 019BD16A47D for ; Tue, 2 Oct 2007 00:14:59 +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 EE77E13C4A3 for ; Tue, 2 Oct 2007 00:14:58 +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 l920C5Lg004084; Mon, 1 Oct 2007 17:12:05 -0700 (PDT) (envelope-from sgk@troutmask.apl.washington.edu) Received: (from sgk@localhost) by troutmask.apl.washington.edu (8.14.1/8.14.1/Submit) id l920Bsu1004083; Mon, 1 Oct 2007 17:11:54 -0700 (PDT) (envelope-from sgk) Date: Mon, 1 Oct 2007 17:11:54 -0700 From: Steve Kargl To: Bruce Evans Message-ID: <20071002001154.GA3782@troutmask.apl.washington.edu> References: <20070928152227.GA39233@troutmask.apl.washington.edu> <20071001173736.U1985@besplex.bde.org> Mime-Version: 1.0 Content-Type: text/plain; charset=us-ascii Content-Disposition: inline In-Reply-To: <20071001173736.U1985@besplex.bde.org> User-Agent: Mutt/1.4.2.3i Cc: freebsd-standards@freebsd.org 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: Tue, 02 Oct 2007 00:14:59 -0000 On Mon, Oct 01, 2007 at 07:56:28PM +1000, Bruce Evans wrote: > On Fri, 28 Sep 2007, Steve Kargl wrote: > > >So, in my never ending crusade to implement missing > >C99 math functions, I decided to tackle sinl, cosl, > >and tanl. I have implementations that use argument > >reduction from k_rem_pio2.c and Chebyshev interpolation > >for the trig functions over the range [0,pi/4]. > > > >On amd64, these routines give <= 1/2 ULP for the various > >values I've tested in the reduced range. Now, if I test > > I doubt this :-). Getting <= 1 ULP is hard; getting <= 1/2 + 0.001 > ULPs is harder, and getting <= 1/2 ULPs without using an enormous > amount of extra precision is a research problem. It appears that I was not clear in what I had intended to write. Let's try again. I've implemented the necessary argument reduction using k_rem_pio2.c for a 64-bit long double. I'm making no claims about its accuracy and the ULP. The comments in k_rem_pio2.c do not discuss the accuracy. It is noted that giving the reduced argument to one of my algorithms does not sudden recover any lost precision. However, given a reduced value in the range [0,pi/4], the algorithms that I have for sinl(), cosl(), and tanl() over this range give <= 1/2 ULP "for the various values I've tested." Note, I did not claim to test *all* values in [0,pi/4]. I've only tested the algorithms with several 10 of millions of randomly generated values. It should also be noted that I'm using Goldberg's definition of ULP where the "exact" result is computed via MPFR with 512 bits of precision and then the mpfr value is converted to long double. /* * y is the approximate value. z is the long double result * of converting the 512 bit MPFR number in round-to-nearest * mode. */ long double ulpl(long double y, long double z) { int n; long double f; f = frexpl(y, &n); f = fabsl(f - ldexpl(z, -n)); f = ldexpl(f, LDBL_MANT_DIG - 1); return (f); } > >these on i386, I see between 400 and 1200 ULP. I spent > >a week or so trying to understand the descrepancy and "fix" > >the problem using a double-double precision arithemetic. > > > >I finally found the problem! /usr/include/float.h on i386 > >is lying about its numerical representation of long double. > > This problem is well known. I must have referred to it many times > in old mail to you. You've told me many things about floating point in FreeBSD. Unfortunately, I have a small brain with only so much capacity. > FreeBSD uses my old workaround for compiler bugs in double > precision, of using a rounding precision of double for everything. I remembered the workaround this weekend, and in fact using fpsetprec to set 64-bit precision on my i386 system generated results consistent with my amd64 system. > Unfortunately, this is still necessary because the > compiler bugs are still there. This mostly breaks long double > precision, and it doesn't fix float precision, but double precision > is more important. libm, in particular k_rem_pio2f.c (which is no > longer used) needs tricky fixes to work around the bugs for float > precision. The fixes have mostly not been merged to the double > precision functions -- we just depend on double precision mostly > working. > > >In particular, at least these 3 values appear to be wrong > > > >#define LDBL_MANT_DIG 64 > >#define LDBL_EPSILON 1.0842021724855044340E-19L > >#define LDBL_DIG 18 > > These were broken in rev.1.9 in 2002. The "slight rounding issues with > long doubles" referred to in the log message are that long doubles are > rounded to double precision when they are operated on, unless you change > the rounding precision to long double using fpsetprec() or fesetprec(). > Even if you change the rounding precision, long doubles can only appear > to work, but not actually work, unless you fix the compiler bugs or > work around them in libraries and in your application. It isn't clear > if the "appear to work" in the log messages is with or without the change > to the rounding precision. As mentioned in the log message, printf() in > 2002 didn't support printing long doubles at all, so it was hard to see > loss of precision in printf output. Now printf supports long doubles and > %a, so it is easy to see lost precision using it. E.g.: > > %%% > #include > #include > > int > main(void) > { > printf(" sin(1) = %La\n", (long double)(double)sin(1)); > printf("2 * sin(1) = %La\n", (long double)(double)(2 * sin(1))); > return (0); > } > %%% > > Output on i386: > > %%% > sin(1) = 0xd.76aa47848677021p-4 > 2 * sin(1) = 0xd.76aa47848677p-3 > %%% > > This exhibits a library bug, several compiler bugs, and the lost precision > from my rounding precision hack: > > o library bug: sin() returns extra precision, despite its prototype saying > that it returns double. Well, I used to think that this is a bug. I > only leared a few months ago that it is a feature, since it matches > the corresponding feature for a function implemented in C. On i386, > sin() is implemented in asm. It just uses the hardware to evaluate > sin() in long double precision, and returns the long double result of > that. The corresponding hardware feature is that in a C function that > does similar calculations is permitted to return extra precision: > > double x, y; > long double foo() { return x+y; } > > This is permitted to do the following: > - evaluate x+y in extra precision > - return the extra precision, since unlike casts and assignments, the > return operation is _not_ required to discard any extra precision. > For gcc on i386, this is what it would do without my rounding precision > hack, except possibly when compiled with -O0 and/or if the result of > x+y is spilled to memory (-O0 makes spilling quite likely). The return > statement must be written as "return (double)(x+y);" to ensure removing > any extra precision _before_ returning. > > Now I wonder if the implicit casts for function args with prototyped > functions are permitted to be as magic as the non-casts for return > statments. Keeping extra precision in intermediate calculations is > good, and you don't really want to lose it unconditionally at function > call boundaries, and it's confusing to lose it on calls but not on > returns. Sorry about the lack of trimming the above passage, but I was afraid I might lose context if I deleted text. Anyway, the above looks like a catch-22 for i386. I can wrap my implementation in fpsetprec pairs to get the desired precision, but someone using libm may unknowingly lose precision. The value of LDBL_MANT_DIG in float.h actually caused my initial post because I knew I lost at least 11 bits of precision (until I remembered your workaround). I thought that my polynomial approximations were in error, so I tried many variations on the Chebyshev interpolation method, Horner's method, summation with compensation, etc. > I'm using sin(1) mainly to generate a long double with nonzero in its > 11 least significant mantissa bits. Note that initializing a long > double to apparently-the-same value at compile time would not work: > > long double x = 0xd.76aa47848677021p-4L; /* loses extra precision */ > > This is because gcc on FreeBSD supports my precision hack to a fault, > and has bugs in the fault. gcc now tries to do compile-time arithmetic > in the same precision as the runtime arithmetic, i.e., in only double > precision on FreeBSD. This is wrong if there is no arithmetic involved, > as in the above initializer. The runtime rounding precision only applies > to operations -- it doesn't prevent loading, storing or non-arithmetic > (function) operations on long doubles with full precision. The above > initializer corresponds to a store of a value that has somehow been > calculated in full precision. > > Note that this bug makes it more difficult to implement functions in > long double precision -- you cannot use tables of long double constants. Ah ha! This explains why my Chebyshev polynomial approximations (ie nearly mimimax approximations) can have between 10 and 100 ULP of errors. > OTOH, loads of long double constants are so slow on i386 that such > constants should almost never be used. Instead represent long doubles > as the sum of 2 doubles if possible. Loading 2 doubles and adding them > is faster than loading 1 long double, at least if the 2 extra operations > can be pipelined effectively. Looks like I spend some time this week splitting the coefficients into high and low parts. Once again, thanks for a detailed explanation of some of the finer points of the 53 vs 64 bit issue. Like many of you other responses, I'll save this one and hopefully retain some of the info. -- Steve From owner-freebsd-standards@FreeBSD.ORG Tue Oct 2 02:25:06 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 74D5B16A418 for ; Tue, 2 Oct 2007 02:25:06 +0000 (UTC) (envelope-from brde@optusnet.com.au) Received: from mail12.syd.optusnet.com.au (mail12.syd.optusnet.com.au [211.29.132.193]) by mx1.freebsd.org (Postfix) with ESMTP id 2102213C467 for ; Tue, 2 Oct 2007 02:25:05 +0000 (UTC) (envelope-from brde@optusnet.com.au) Received: from besplex.bde.org (c220-239-235-248.carlnfd3.nsw.optusnet.com.au [220.239.235.248]) by mail12.syd.optusnet.com.au (8.13.1/8.13.1) with ESMTP id l922P1cR000332 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=NO); Tue, 2 Oct 2007 12:25:03 +1000 Date: Tue, 2 Oct 2007 12:25:01 +1000 (EST) From: Bruce Evans X-X-Sender: bde@besplex.bde.org To: Steve Kargl In-Reply-To: <20071002001154.GA3782@troutmask.apl.washington.edu> Message-ID: <20071002105817.V7990@besplex.bde.org> References: <20070928152227.GA39233@troutmask.apl.washington.edu> <20071001173736.U1985@besplex.bde.org> <20071002001154.GA3782@troutmask.apl.washington.edu> MIME-Version: 1.0 Content-Type: TEXT/PLAIN; charset=US-ASCII; format=flowed Cc: freebsd-standards@freebsd.org 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: Tue, 02 Oct 2007 02:25:06 -0000 On Mon, 1 Oct 2007, Steve Kargl wrote: > On Mon, Oct 01, 2007 at 07:56:28PM +1000, Bruce Evans wrote: >> On Fri, 28 Sep 2007, Steve Kargl wrote: >>> On amd64, these routines give <= 1/2 ULP for the various >>> values I've tested in the reduced range. Now, if I test >> >> I doubt this :-). Getting <= 1 ULP is hard; getting <= 1/2 + 0.001 >> ULPs is harder, and getting <= 1/2 ULPs without using an enormous >> amount of extra precision is a research problem. > > It appears that I was not clear in what I had intended to write. > Let's try again. I've implemented the necessary argument > reduction using k_rem_pio2.c for a 64-bit long double. I'm > making no claims about its accuracy and the ULP. The comments > in k_rem_pio2.c do not discuss the accuracy. It is noted > that giving the reduced argument to one of my algorithms > does not sudden recover any lost precision. I think it is accurate to within one bit provided it is set up correctly. Setup includes passing it a large table of bits of 2/pi. I don't know if the existing tables with 1584 bits are large enough for long double precision. If long doubles have 64 bits, then the recipe in k_rem_pio_2 returns 106 mantissa bits. I would expect a max relative error of 2^-105. > However, given a reduced value in the range [0,pi/4], the > algorithms that I have for sinl(), cosl(), and tanl() over > this range give <= 1/2 ULP "for the various values I've > tested." Note, I did not claim to test *all* values in > [0,pi/4]. I've only tested the algorithms with several 10 > of millions of randomly generated values. If you start with a relative error of 2^-105, then it is not very hard to get a relative error of 2^-104 in the result, but it is not very easy to get this either without having very slow, large and/or tricky code. k_cos.c only wants a relative error of about 2^-60, so it only needs a small amount of tricky code to be fast (it uses essentially cos(x+y) ~= cos(x) - x*y; getting a smaller relative error would require many more terms; everything should really be evaluated in double-double precision, but k_cos.c uses tricks so ordinary double precision error can be used to get a relative error of strictly less than 2^-54 (half an ulp) if the final addition were in infinite precision. Rounding then gives a final relative error of strictly less than 1 ulp. If you can limit the relative error to < 2^-104 before rounding to long double precision, then (with i386 long doubles but of course not with sparc64 long doubles) the final result has a very large chance of being perfectly rounded, but since the result space is huge there is also a very large chance of wrong rounding in a huge number of cases. (Chance of wrong rounding ~= 2^-40 in an arg space of size 2^64 gives about 2^20 cases incorrectly rounded. The arg space is so large that you will probably never find any incorrectly rounded cases by blind searching.) > It should also be noted that I'm using Goldberg's definition > of ULP where the "exact" result is computed via MPFR with > 512 bits of precision and then the mpfr value is converted > to long double. If you can limit the relative error to < 2^-512 before the final rounding step, then you probably have perfect rounding in all cases but no one can prove it because the arg space is too large. (Chance of wrong rounding ~= 2^-448 in an arg space of size 2^64 gives about 2^-384 cases incorrectly rounded). I think doing any calculations in 512-bit precision would be too slow to be useful for long double precision, except for verifying results. > /* > * y is the approximate value. z is the long double result > * of converting the 512 bit MPFR number in round-to-nearest > * mode. > */ > long double > ulpl(long double y, long double z) > { > int n; > long double f; > f = frexpl(y, &n); > f = fabsl(f - ldexpl(z, -n)); > f = ldexpl(f, LDBL_MANT_DIG - 1); > return (f); > } You should try to use all the bits in the MFPR number, to verify that when the rounded results differ, it is because the result is near a half-way case (as indicated by lots of 1's bits in the MFPR number after the rounding point). >> [rounding precision hack] > Sorry about the lack of trimming the above passage, but I was afraid > I might lose context if I deleted text. All gone now :-). > Anyway, the above looks like a catch-22 for i386. I can wrap In 1988 I thought it would be fixed by now. Now It looks like i386 will be obsolete first. >> Note that this bug makes it more difficult to implement functions in >> long double precision -- you cannot use tables of long double constants. > > Ah ha! This explains why my Chebyshev polynomial approximations > (ie nearly mimimax approximations) can have between 10 and 100 > ULP of errors. And not all 11 bits (2048 ulps) are normally lost since the constants are usually exact for leading terms. >> OTOH, loads of long double constants are so slow on i386 that such >> constants should almost never be used. Instead represent long doubles >> as the sum of 2 doubles if possible. Loading 2 doubles and adding them >> is faster than loading 1 long double, at least if the 2 extra operations >> can be pipelined effectively. > > Looks like I spend some time this week splitting the coefficients > into high and low parts. I use pari scripts to generate coeffiencts, and they handle the splitting with minor adjustements. Bruce From owner-freebsd-standards@FreeBSD.ORG Tue Oct 2 08:20:04 2007 Return-Path: Delivered-To: freebsd-standards@hub.freebsd.org Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2001:4f8:fff6::34]) by hub.freebsd.org (Postfix) with ESMTP id 9E0F016A41B; Tue, 2 Oct 2007 08:20:04 +0000 (UTC) (envelope-from jinmei@FreeBSD.org) Received: from freefall.freebsd.org (freefall.freebsd.org [IPv6:2001:4f8:fff6::28]) by mx1.freebsd.org (Postfix) with ESMTP id 7107613C469; Tue, 2 Oct 2007 08:20:04 +0000 (UTC) (envelope-from jinmei@FreeBSD.org) Received: from freefall.freebsd.org (jinmei@localhost [127.0.0.1]) by freefall.freebsd.org (8.14.1/8.14.1) with ESMTP id l928K4Dt040687; Tue, 2 Oct 2007 08:20:04 GMT (envelope-from jinmei@freefall.freebsd.org) Received: (from jinmei@localhost) by freefall.freebsd.org (8.14.1/8.14.1/Submit) id l928K3qX040675; Tue, 2 Oct 2007 08:20:03 GMT (envelope-from jinmei) Date: Tue, 2 Oct 2007 08:20:03 GMT Message-Id: <200710020820.l928K3qX040675@freefall.freebsd.org> To: marka@isc.org, jinmei@FreeBSD.org, freebsd-standards@FreeBSD.org From: jinmei@FreeBSD.org Cc: Subject: Re: standards/114910: getaddrinfo() fails to set ai_canonname 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: Tue, 02 Oct 2007 08:20:04 -0000 Synopsis: getaddrinfo() fails to set ai_canonname State-Changed-From-To: feedback->closed State-Changed-By: jinmei State-Changed-When: Tue Oct 2 08:19:32 UTC 2007 State-Changed-Why: no objection to the committed fix. http://www.freebsd.org/cgi/query-pr.cgi?pr=114910 From owner-freebsd-standards@FreeBSD.ORG Tue Oct 2 13:40:03 2007 Return-Path: Delivered-To: freebsd-standards@hub.freebsd.org Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2001:4f8:fff6::34]) by hub.freebsd.org (Postfix) with ESMTP id DD01416A468 for ; Tue, 2 Oct 2007 13:40:03 +0000 (UTC) (envelope-from gnats@FreeBSD.org) Received: from freefall.freebsd.org (freefall.freebsd.org [IPv6:2001:4f8:fff6::28]) by mx1.freebsd.org (Postfix) with ESMTP id A9E4513C480 for ; Tue, 2 Oct 2007 13:40:03 +0000 (UTC) (envelope-from gnats@FreeBSD.org) Received: from freefall.freebsd.org (gnats@localhost [127.0.0.1]) by freefall.freebsd.org (8.14.1/8.14.1) with ESMTP id l92De3IL061457 for ; Tue, 2 Oct 2007 13:40:03 GMT (envelope-from gnats@freefall.freebsd.org) Received: (from gnats@localhost) by freefall.freebsd.org (8.14.1/8.14.1/Submit) id l92De3hE061456; Tue, 2 Oct 2007 13:40:03 GMT (envelope-from gnats) Resent-Date: Tue, 2 Oct 2007 13:40:03 GMT Resent-Message-Id: <200710021340.l92De3hE061456@freefall.freebsd.org> Resent-From: FreeBSD-gnats-submit@FreeBSD.org (GNATS Filer) Resent-To: freebsd-standards@FreeBSD.org Resent-Reply-To: FreeBSD-gnats-submit@FreeBSD.org, Roy Marples Received: from mx1.freebsd.org (mx1.freebsd.org [IPv6:2001:4f8:fff6::34]) by hub.freebsd.org (Postfix) with ESMTP id 6787216A41A for ; Tue, 2 Oct 2007 13:30:56 +0000 (UTC) (envelope-from nobody@FreeBSD.org) Received: from www.freebsd.org (www.freebsd.org [IPv6:2001:4f8:fff6::21]) by mx1.freebsd.org (Postfix) with ESMTP id 2FF5A13C447 for ; Tue, 2 Oct 2007 13:30:56 +0000 (UTC) (envelope-from nobody@FreeBSD.org) Received: from www.freebsd.org (localhost [127.0.0.1]) by www.freebsd.org (8.14.1/8.14.1) with ESMTP id l92DUu22001703 for ; Tue, 2 Oct 2007 13:30:56 GMT (envelope-from nobody@www.freebsd.org) Received: (from nobody@localhost) by www.freebsd.org (8.14.1/8.14.1/Submit) id l92DUt2Y001701; Tue, 2 Oct 2007 13:30:55 GMT (envelope-from nobody) Message-Id: <200710021330.l92DUt2Y001701@www.freebsd.org> Date: Tue, 2 Oct 2007 13:30:55 GMT From: Roy Marples To: freebsd-gnats-submit@FreeBSD.org X-Send-Pr-Version: www-3.1 Cc: Subject: standards/116826: [PATCH] sh support for POSIX character classes 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: Tue, 02 Oct 2007 13:40:03 -0000 >Number: 116826 >Category: standards >Synopsis: [PATCH] sh support for POSIX character classes >Confidential: no >Severity: non-critical >Priority: low >Responsible: freebsd-standards >State: open >Quarter: >Keywords: >Date-Required: >Class: sw-bug >Submitter-Id: current-users >Arrival-Date: Tue Oct 02 13:40:03 GMT 2007 >Closed-Date: >Last-Modified: >Originator: Roy Marples >Release: FreeBSD-6.2 >Organization: Gentoo >Environment: FreeBSD uberlaptop 6.2-RELEASE FreeBSD Gentoo 6.2-r3 #0: Thu Sep 13 17:44:37 BST 2007 root@uberlaptop:/usr/src/sys-6.2-r3/i386/compile/UBERKERNEL i386 >Description: sh does not support POSIX character classes http://www.opengroup.org/onlinepubs/009695399/utilities/xcu_chap02.html#tag_02_13 http://www.opengroup.org/onlinepubs/009695399/basedefs/xbd_chap09.html#tag_09_03_05 >How-To-Repeat: $ sh -c 'case foo in [[:alpha:]]*) echo "alpha";; *) echo "not alpha";; esac' not alpha I would expect it to echo alpha. >Fix: Patch attached with submission follows: diff -u a/sh/expand.c b/sh/expand.c --- a/sh/expand.c 2005-11-06 20:39:47 +0000 +++ b/sh/expand.c 2007-10-02 13:46:28 +0100 @@ -1320,6 +1320,42 @@ } +STATIC int ccmatch(char *p, int chr, char **r) +{ + static const struct class { + char name[10]; + int (*fn)(int); + } classes[] = { + { .name = ":alnum:]", .fn = isalnum }, + { .name = ":cntrl:]", .fn = iscntrl }, + { .name = ":lower:]", .fn = islower }, + { .name = ":space:]", .fn = isspace }, + { .name = ":alpha:]", .fn = isalpha }, + { .name = ":digit:]", .fn = isdigit }, + { .name = ":print:]", .fn = isprint }, + { .name = ":upper:]", .fn = isupper }, + { .name = ":blank:]", .fn = isblank }, + { .name = ":graph:]", .fn = isgraph }, + { .name = ":punct:]", .fn = ispunct }, + { .name = ":xdigit:]", .fn = isxdigit }, + }; + const struct class *class, *end; + char *q; + + end = classes + sizeof(classes) / sizeof(classes[0]); + for (class = classes; class < end; class++) { + q = prefix(class->name, p); + if (!q) + continue; + *r = q; + return class->fn(chr); + } + + *r = 0; + return 0; +} + + STATIC int pmatch(char *pattern, char *string, int squoted) { @@ -1405,6 +1441,15 @@ continue; if (c == CTLESC) c = *p++; + else if (c == '[') { + char *r; + + found |= ccmatch(p, chr, &r); + if (r) { + p = r; + continue; + } + } if (*p == '-' && p[1] != ']') { p++; while (*p == CTLQUOTEMARK) diff -u a/sh/mystring.c b/sh/mystring.c --- a/sh/mystring.c 2004-04-06 21:06:51 +0100 +++ b/sh/mystring.c 2007-10-02 13:45:31 +0100 @@ -88,14 +88,14 @@ * prefix -- see if pfx is a prefix of string. */ -int +char * prefix(const char *pfx, const char *string) { while (*pfx) { if (*pfx++ != *string++) return 0; } - return 1; + return (char *)string; } diff -u a/sh/mystring.h b/sh/mystring.h --- a/sh/mystring.h 2004-04-06 21:06:51 +0100 +++ b/sh/mystring.h 2007-10-02 13:45:35 +0100 @@ -36,7 +36,7 @@ #include void scopyn(const char *, char *, int); -int prefix(const char *, const char *); +char *prefix(const char *, const char *); int number(const char *); int is_number(const char *); >Release-Note: >Audit-Trail: >Unformatted: From owner-freebsd-standards@FreeBSD.ORG Tue Oct 2 17:35:37 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 DEA5316A417 for ; Tue, 2 Oct 2007 17:35:37 +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 9773D13C461 for ; Tue, 2 Oct 2007 17:35:37 +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 l92HWbWO012763; Tue, 2 Oct 2007 10:32:37 -0700 (PDT) (envelope-from sgk@troutmask.apl.washington.edu) Received: (from sgk@localhost) by troutmask.apl.washington.edu (8.14.1/8.14.1/Submit) id l92HWbLx012762; Tue, 2 Oct 2007 10:32:37 -0700 (PDT) (envelope-from sgk) Date: Tue, 2 Oct 2007 10:32:37 -0700 From: Steve Kargl To: Bruce Evans , freebsd-standards@FreeBSD.ORG Message-ID: <20071002173237.GA12586@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> Mime-Version: 1.0 Content-Type: text/plain; charset=us-ascii Content-Disposition: inline In-Reply-To: <20071002172317.GA95181@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: Tue, 02 Oct 2007 17:35:38 -0000 On Tue, Oct 02, 2007 at 01:23:18PM -0400, David Schultz wrote: > > Although it would be nice to get all this stuff right the first > time, very few people are going to care if our trig functions are > accurate to within 1 ulp for huge inputs; many competing math > libraries don't guarantee that anyway. A programmer who asks for > sinl(1000000000*PI + 0.01) is going to be disappointed regardless, > because you can't represent the input accurately using IEEE-754 > floating point. Most people just care about taking a program that > uses sinl() and getting it to compile and run on FreeBSD, and most > of those programs don't call sinl() with huge arguments. > > Anyway, my point is that if you have something that works > reasonably well and doesn't have egregious errors, my suggestion > is to just commit it and not kill yourself over whether the > argument reduction is correct in the last ulp. Hi David, There are a few problems: 1) I don't have commit access. 2) I need to do style(9) clean-up pass. 3) I've only tested these on i386 and amd64, and I know these fail for sparc64. 4) Most importantly, I don't know how the project wants to handle the 53 vs 64 bit default fpu setting on i386. PS: There is also the minor possibility of giving bde either a heartache or death via laughter when he finally sees what I have done. :) -- Steve From owner-freebsd-standards@FreeBSD.ORG Tue Oct 2 17:59:55 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 F050B16A418 for ; Tue, 2 Oct 2007 17:59:55 +0000 (UTC) (envelope-from das@FreeBSD.ORG) Received: from VARK.MIT.EDU (VARK.MIT.EDU [18.95.3.179]) by mx1.freebsd.org (Postfix) with ESMTP id 9546413C45A for ; Tue, 2 Oct 2007 17:59:55 +0000 (UTC) (envelope-from das@FreeBSD.ORG) Received: from VARK.MIT.EDU (localhost [127.0.0.1]) by VARK.MIT.EDU (8.14.1/8.14.1) with ESMTP id l92HNIpL095269; Tue, 2 Oct 2007 13:23:18 -0400 (EDT) (envelope-from das@FreeBSD.ORG) Received: (from das@localhost) by VARK.MIT.EDU (8.14.1/8.14.1/Submit) id l92HNIPq095268; Tue, 2 Oct 2007 13:23:18 -0400 (EDT) (envelope-from das@FreeBSD.ORG) Date: Tue, 2 Oct 2007 13:23:18 -0400 From: David Schultz To: Steve Kargl Message-ID: <20071002172317.GA95181@VARK.MIT.EDU> Mail-Followup-To: Steve Kargl , Bruce Evans , freebsd-standards@FreeBSD.ORG References: <20070928152227.GA39233@troutmask.apl.washington.edu> <20071001173736.U1985@besplex.bde.org> <20071002001154.GA3782@troutmask.apl.washington.edu> Mime-Version: 1.0 Content-Type: text/plain; charset=us-ascii Content-Disposition: inline In-Reply-To: <20071002001154.GA3782@troutmask.apl.washington.edu> Cc: freebsd-standards@FreeBSD.ORG 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: Tue, 02 Oct 2007 17:59:56 -0000 Just a quick note... Although it would be nice to get all this stuff right the first time, very few people are going to care if our trig functions are accurate to within 1 ulp for huge inputs; many competing math libraries don't guarantee that anyway. A programmer who asks for sinl(1000000000*PI + 0.01) is going to be disappointed regardless, because you can't represent the input accurately using IEEE-754 floating point. Most people just care about taking a program that uses sinl() and getting it to compile and run on FreeBSD, and most of those programs don't call sinl() with huge arguments. Anyway, my point is that if you have something that works reasonably well and doesn't have egregious errors, my suggestion is to just commit it and not kill yourself over whether the argument reduction is correct in the last ulp. From owner-freebsd-standards@FreeBSD.ORG Tue Oct 2 21:15:07 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 E799316A420 for ; Tue, 2 Oct 2007 21:15:07 +0000 (UTC) (envelope-from das@FreeBSD.ORG) Received: from VARK.MIT.EDU (VARK.MIT.EDU [18.95.3.179]) by mx1.freebsd.org (Postfix) with ESMTP id 0C38E13C4A3 for ; Tue, 2 Oct 2007 21:15:06 +0000 (UTC) (envelope-from das@FreeBSD.ORG) Received: from VARK.MIT.EDU (localhost [127.0.0.1]) by VARK.MIT.EDU (8.14.1/8.14.1) with ESMTP id l92LEvkp096756; Tue, 2 Oct 2007 17:14:57 -0400 (EDT) (envelope-from das@FreeBSD.ORG) Received: (from das@localhost) by VARK.MIT.EDU (8.14.1/8.14.1/Submit) id l92LEvW0096755; Tue, 2 Oct 2007 17:14:57 -0400 (EDT) (envelope-from das@FreeBSD.ORG) Date: Tue, 2 Oct 2007 17:14:57 -0400 From: David Schultz To: Steve Kargl Message-ID: <20071002211457.GA96707@VARK.MIT.EDU> Mail-Followup-To: Steve Kargl , Bruce Evans , freebsd-standards@FreeBSD.ORG 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> Mime-Version: 1.0 Content-Type: text/plain; charset=us-ascii Content-Disposition: inline In-Reply-To: <20071002173237.GA12586@troutmask.apl.washington.edu> Cc: freebsd-standards@FreeBSD.ORG 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: Tue, 02 Oct 2007 21:15:08 -0000 On Tue, Oct 02, 2007, Steve Kargl wrote: > There are a few problems: 1) I don't have commit access. 2) I > need to do style(9) clean-up pass. 3) I've only tested these > on i386 and amd64, and I know these fail for sparc64. 4) Most > importantly, I don't know how the project wants to handle the > 53 vs 64 bit default fpu setting on i386. > > PS: There is also the minor possibility of giving bde either > a heartache or death via laughter when he finally sees what I > have done. :) I have some other stuff in my tree that I need to commit at some point, so I can throw it your changes too. It needs to wait until after SOSP though. We should see about getting you commit access some time. I don't have time to be a very good mentor at the moment, but maybe someone else could. For the 128-bit format used on sparc64, you probably need a separate implementation (or at least different constants), but hopefully that's not too hard once you've done it once. Dealing with our reduced precision on x86 is easy, since in that case you can just fall back on the ordinary sin() routine. From owner-freebsd-standards@FreeBSD.ORG Wed Oct 3 01:13:44 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 E65A716A419 for ; Wed, 3 Oct 2007 01:13:43 +0000 (UTC) (envelope-from brde@optusnet.com.au) Received: from mail16.syd.optusnet.com.au (mail16.syd.optusnet.com.au [211.29.132.197]) by mx1.freebsd.org (Postfix) with ESMTP id 8E8B413C448 for ; Wed, 3 Oct 2007 01:13:43 +0000 (UTC) (envelope-from brde@optusnet.com.au) Received: from c220-239-235-248.carlnfd3.nsw.optusnet.com.au (c220-239-235-248.carlnfd3.nsw.optusnet.com.au [220.239.235.248]) by mail16.syd.optusnet.com.au (8.13.1/8.13.1) with ESMTP id l931Dc02020774 (version=TLSv1/SSLv3 cipher=DHE-RSA-AES256-SHA bits=256 verify=NO); Wed, 3 Oct 2007 11:13:40 +1000 Date: Wed, 3 Oct 2007 11:13:38 +1000 (EST) From: Bruce Evans X-X-Sender: bde@delplex.bde.org To: Steve Kargl In-Reply-To: <20071002173237.GA12586@troutmask.apl.washington.edu> Message-ID: <20071003103519.X14175@delplex.bde.org> 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> MIME-Version: 1.0 Content-Type: TEXT/PLAIN; charset=US-ASCII; format=flowed Cc: freebsd-standards@FreeBSD.org 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, 03 Oct 2007 01:13:44 -0000 On Tue, 2 Oct 2007, Steve Kargl wrote: > On Tue, Oct 02, 2007 at 01:23:18PM -0400, David Schultz wrote: >> ... >> Anyway, my point is that if you have something that works >> reasonably well and doesn't have egregious errors, my suggestion >> is to just commit it and not kill yourself over whether the >> argument reduction is correct in the last ulp. > There are a few problems: 1) I don't have commit access. 2) I > need to do style(9) clean-up pass. 3) I've only tested these > on i386 and amd64, and I know these fail for sparc64. 4) Most > importantly, I don't know how the project wants to handle the > 53 vs 64 bit default fpu setting on i386. (3) is most important. The C versions will never even be used for i386 or amd64, since these arches have sin and cos in hardware and (I believe) it is impossible to beat the hardware on these arches. The asm versions have always been trivial to implement (change 2 bytes in s_sin.S (compare e_sqrt.S with e_sqrtf.S; the long double changes go the other way). amd64 actually uses the C versions for double precision becase the C versions are competitive with the hardware for speed and beat it for accuracy on large args and on small args aear a multiple of pi/2. (4): These functions are useless unless the default is changed. Users have always been able to get equivalent _in_accuracy on i386 by just using the double precision versions. Hacks like #define sinl sin almost work for just getting things to compile. In fact, on i386, due to the bugfeatures outlined in my previous mail, users can even get full long double accuracy using hacks like that, even if the rounding precision is 53 bits: %%% #include long double x: main() { x = (long double)sin(1); } %%% On i386, this stores the full long double result for sinl(1) in x, due to the bugfeatures -- sin(1) returns extra precision in a register, and gcc doesn't understand the extra precision so it just stores all bits of the register in x. The behaviour would be similar if sin(1) were used in an expression. On amd64, the ABI prevents this from working (doubles must be returned in a 64-bit XMM register). > PS: There is also the minor possibility of giving bde either > a heartache or death via laughter when he finally sees what I > have done. :) I would be unhappy with any versions that aren't either: - as identical as possible in algorithm and style to the double precision versions, in the same way that the float versions are as identical as possible, or - much better in algorithm and style and efficiency and accuracy than the double precision versions. I don't want to maintain different versions. I tried to turn your original efforts for logl() into the latter for logf(), log() and logl(), but so far they are only a little better. In particular they are not good enough to commit due to point (3) -- the long double version is only implemented for machines with 80-bit long doubles, and all versions only beat the hardware for efficiency on some machines (amd64, but not i386 with P[2-4]). I think the long double version beats the hardware for accuracy but haven't finished testing this. The hardware is much easier to beat for efficiency for the log family than for the sin family. The hardware is much easier to beat for accuracy for logl() than for anything else. The hardware is accurate enough for logf() and logl(), so is unbeatable for accuracy of these functions unless a slow version is used to get more accuracy. Bruce