From owner-cvs-all@FreeBSD.ORG Sat Dec 3 07:23:32 2005 Return-Path: X-Original-To: cvs-all@FreeBSD.org Delivered-To: cvs-all@FreeBSD.org Received: from mx1.FreeBSD.org (mx1.freebsd.org [216.136.204.125]) by hub.freebsd.org (Postfix) with ESMTP id ED84F16A41F; Sat, 3 Dec 2005 07:23:31 +0000 (GMT) (envelope-from bde@FreeBSD.org) Received: from repoman.freebsd.org (repoman.freebsd.org [216.136.204.115]) by mx1.FreeBSD.org (Postfix) with ESMTP id 37B6B43D5C; Sat, 3 Dec 2005 07:23:31 +0000 (GMT) (envelope-from bde@FreeBSD.org) Received: from repoman.freebsd.org (localhost [127.0.0.1]) by repoman.freebsd.org (8.13.1/8.13.1) with ESMTP id jB37NU6A020627; Sat, 3 Dec 2005 07:23:30 GMT (envelope-from bde@repoman.freebsd.org) Received: (from bde@localhost) by repoman.freebsd.org (8.13.1/8.13.1/Submit) id jB37NUIM020626; Sat, 3 Dec 2005 07:23:30 GMT (envelope-from bde) Message-Id: <200512030723.jB37NUIM020626@repoman.freebsd.org> From: Bruce Evans Date: Sat, 3 Dec 2005 07:23:30 +0000 (UTC) To: src-committers@FreeBSD.org, cvs-src@FreeBSD.org, cvs-all@FreeBSD.org X-FreeBSD-CVS-Branch: HEAD Cc: Subject: cvs commit: src/lib/msun/src s_rint.c X-BeenThere: cvs-all@freebsd.org X-Mailman-Version: 2.1.5 Precedence: list List-Id: CVS commit messages for the entire tree List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Sat, 03 Dec 2005 07:23:32 -0000 bde 2005-12-03 07:23:30 UTC FreeBSD src repository Modified files: lib/msun/src s_rint.c Log: Fixed rint(x) in the following cases: (1) In round-to-nearest mode, on all machines, fdlibm rint() never worked for |x| = n+0.75 where n is an even integer between 262144 and 524286 inclusive (2*131072 cases). To avoid double rounding on some machines, we begin by adjusting x to a value with the 0.25 bit not set, essentially by moving the 0.25 bit to a lower bit where it works well enough as a guard, but we botched the adjustment when log2(|x|) == 18 (2*2**52 cases) and ended up just clearing the 0.25 bit then. Most subcases still worked accidentally since another lower bit serves as a guard. The case of odd n worked accidentally because the rounding goes the right way then. However, for even n, after mangling n+0.75 to 0.5, rounding gives n but the correct result is n+1. (2) In round-towards-minus-infinity mode, on all machines, fdlibm rint() never for x = n+0.25 where n is any integer between -524287 and -262144 inclusive (262144 cases). In these cases, after mangling n+0.25 to n, rounding gives n but the correct result is n-1. (3) In round-towards-plus-infinity mode, on all machines, fdlibm rint() never for x = n+0.25 where n is any integer between 262144 and 524287 inclusive (262144 cases). In these cases, after mangling n+0.25 to n, rounding gives n but the correct result is n+1. A variant of this bug was fixed for the float case in rev.1.9 of s_rintf.c, but the analysis there is incomplete (it only mentions (1)) and the fix is buggy. Example of the problem with double rounding: rint(1.375) on a machine which evaluates double expressions with just 1 bit of extra precision and is in round-to-nearest mode. We evaluate the result using (double)(2**52 + 1.375) - 2**52. Evaluating 2**52 + 1.375 in (53+1) bit prcision gives 2**52 + 1.5 (first rounding). (Second) rounding of this to double gives 2**52 + 2.0. Subtracting 2**52 from this gives 2.0 but we want 1.0. Evaluating 2**52 + 1.375 in double precision would have given the desired intermediate result of 2**52 + 1.0. The double rounding problem is relatively rare, so the botched adjustment can be fixed for most machines by removing the entire adjustment. This would be a wrong fix (using it is 1 of the bugs in rev.1.9 of s_rintf.c) since fdlibm is supposed to be generic, but it works in the following cases: - on all machines that evaluate double expressions in double precision, provided either long double has the same precision as double (alpha, and i386's with precision forced to double) or my earlier fix to use a long double 2**52 is modified to avoid using long double precision. - on all machines that evaluate double expressions in many more than 11 bits of extra precision. The 1 bit of extra precision in the example is the worst case. With N bits of extra precision, it sufices to adjust the bit N bits below the 0.5 bit. For N >= about 52 there is no such bit so the adjustment is both impossible and unnecessary. The fix in rev.1.9 of s_rintf.c apparently depends on corresponding magic in float precision: on all supported machines N is either 0 or >= 24, so double rounding doesn't occur in practice. - on all machines that don't use fdlibm rint*() (i386's). So under FreeBSD, the double rounding problem only affects amd64 now, but should only affect i386 in future (when double expressions are evaluated in long double precision). Revision Changes Path 1.12 +9 -0 src/lib/msun/src/s_rint.c