From mboxrd@z Thu Jan 1 00:00:00 1970 Path: main.gmane.org!not-for-mail From: Kevin Ryde Newsgroups: gmane.lisp.guile.devel Subject: scm_round certain wrong results Date: Sun, 18 Apr 2004 08:56:32 +1000 Sender: guile-devel-bounces+guile-devel=m.gmane.org@gnu.org Message-ID: <878ygurz6n.fsf@zip.com.au> NNTP-Posting-Host: deer.gmane.org Mime-Version: 1.0 Content-Type: multipart/mixed; boundary="=-=-=" X-Trace: sea.gmane.org 1082242937 607 80.91.224.253 (17 Apr 2004 23:02:17 GMT) X-Complaints-To: usenet@sea.gmane.org NNTP-Posting-Date: Sat, 17 Apr 2004 23:02:17 +0000 (UTC) Original-X-From: guile-devel-bounces+guile-devel=m.gmane.org@gnu.org Sun Apr 18 01:02:09 2004 Return-path: Original-Received: from monty-python.gnu.org ([199.232.76.173]) by deer.gmane.org with esmtp (Exim 3.35 #1 (Debian)) id 1BEyp6-0001K1-00 for ; Sun, 18 Apr 2004 01:02:09 +0200 Original-Received: from localhost ([127.0.0.1] helo=monty-python.gnu.org) by monty-python.gnu.org with esmtp (Exim 4.30) id 1BEykn-00053c-SF for guile-devel@m.gmane.org; Sat, 17 Apr 2004 18:57:41 -0400 Original-Received: from list by monty-python.gnu.org with tmda-scanned (Exim 4.30) id 1BEykP-00053M-Jy for guile-devel@gnu.org; Sat, 17 Apr 2004 18:57:17 -0400 Original-Received: from mail by monty-python.gnu.org with spam-scanned (Exim 4.30) id 1BEyjr-0004jo-Vk for guile-devel@gnu.org; Sat, 17 Apr 2004 18:57:15 -0400 Original-Received: from [61.8.0.85] (helo=mailout2.pacific.net.au) by monty-python.gnu.org with esmtp (Exim 4.30) id 1BEyjq-0004iL-VQ for guile-devel@gnu.org; Sat, 17 Apr 2004 18:56:43 -0400 Original-Received: from mailproxy2.pacific.net.au (mailproxy2.pacific.net.au [61.8.0.87]) by mailout2.pacific.net.au (8.12.3/8.12.3/Debian-6.6) with ESMTP id i3HMue5v021069 for ; Sun, 18 Apr 2004 08:56:40 +1000 Original-Received: from localhost (ppp25C7.dyn.pacific.net.au [61.8.37.199]) by mailproxy2.pacific.net.au (8.12.3/8.12.3/Debian-6.6) with ESMTP id i3HMucHV008046 for ; Sun, 18 Apr 2004 08:56:38 +1000 Original-Received: from gg by localhost with local (Exim 3.36 #1 (Debian)) id 1BEyjg-0005MU-00; Sun, 18 Apr 2004 08:56:32 +1000 Original-To: guile-devel@gnu.org Mail-Copies-To: never User-Agent: Gnus/5.110002 (No Gnus v0.2) Emacs/21.3 (gnu/linux) X-BeenThere: guile-devel@gnu.org X-Mailman-Version: 2.1.4 Precedence: list List-Id: "Developers list for Guile, the GNU extensibility library" List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , Errors-To: guile-devel-bounces+guile-devel=m.gmane.org@gnu.org Xref: main.gmane.org gmane.lisp.guile.devel:3606 X-Report-Spam: http://spam.gmane.org/gmane.lisp.guile.devel:3606 --=-=-= * numbers.c (scm_round): Test for x already an integer, to avoid bad rounding in x+0.5 when x is a big value already an integer. In various modes x+0.5 could give an adjacent integer, leading to that as the result, when we really just wanted x itself. * standalone/test-round.c: New file, exercising scm_round. * standalone/Makefile.am: Add it. * configure.in (AC_CHECK_FUNCS): Add fesetround. (AC_CHECK_HEADERS): Add fenv.h. The test program exercises the bad cases. In the usual default nearest-even hardware rounding I think x = +/- (2^53-1) is the only bad case. But if an application has gone to for instance round-upwards mode then lots of other big values come out wrong too. I believe it's worth scm_round ensuring it's not influenced by the hardware rounding mode, even if that's not something guile normally fiddles with. This would be for the 1.6 branch too. scm_round is used for `round' there, you can see a problem for instance with (define x #i#x1FFFFFFFFFFFFF) (integer? x) => #t (= x (round x)) => #f If you're doing this on an i386 or m68k you need numbers.c compiled unoptimized, since optimized will hold values in the coprocessor, giving enough extra precision to avoid this case (but not others). On other cpus I think it will happen always. --=-=-= Content-Disposition: inline; filename=numbers.c.round.diff --- numbers.c.~1.236.~ 2004-04-18 07:55:31.000000000 +1000 +++ numbers.c 2004-04-18 08:29:42.000000000 +1000 @@ -4884,11 +4884,42 @@ #endif } +/* scm_round is done using floor(x+0.5) to round to nearest with half-way + (ie. when x is an integer plus 0.5) cases going upwards. Then half-way + cases are identified and adjusted down if the round-upwards didn't give + the desired even integer. + + "plus_half == result" identifies a half-way case. If plus_half, which is + x + 0.5, is an integer then x must be an integer plus 0.5. + + An odd "result" value is identified with result/2 != floor(result/2). + This is done with plus_half, since that value is ready for use sooner in + a pipelined cpu, and we're already requiring plus_half == result. + + Note however that we need to be careful when x is big and already an + integer. In that case "x+0.5" may round to an adjacent integer, causing + us to return such a value, incorrectly. For instance if the hardware is + in the usual default nearest-even rounding, then for x = 0x1FFFFFFFFFFFFF + (ie. 53 one bits) we will have x+0.5 = 0x20000000000000 and that value + returned. Or if the hardware is in round-upwards mode, then other bigger + values like say x == 2^128 will see x+0.5 rounding up to the next higher + representable value, 2^128+2^76 (or whatever), again incorrect. + + These bad roundings of x+0.5 are avoided by testing at the start whether + x is already an integer. If it is then clearly that's the desired result + already. And if it's not then the exponent must be small enough to allow + an 0.5 to be represented, and hence added without a bad rounding. */ + double scm_round (double x) { - double plus_half = x + 0.5; - double result = floor (plus_half); + double plus_half, result; + + if (x == floor (x)) + return x; + + plus_half = x + 0.5; + result = floor (plus_half); /* Adjust so that the scm_round is towards even. */ return ((plus_half == result && plus_half / 2 != floor (plus_half / 2)) ? result - 1 --=-=-= Content-Type: text/x-csrc Content-Disposition: inline; filename=test-round.c /* Copyright (C) 2004 Free Software Foundation, Inc. * * This library is free software; you can redistribute it and/or * modify it under the terms of the GNU Lesser General Public * License as published by the Free Software Foundation; either * version 2.1 of the License, or (at your option) any later version. * * This library is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU * Lesser General Public License for more details. * * You should have received a copy of the GNU Lesser General Public * License along with this library; if not, write to the Free Software * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA */ #include "config.h" #include #include #include #if HAVE_FENV_H #include #endif #include "libguile.h" #define numberof(x) (sizeof (x) / sizeof ((x)[0])) static void test_scm_round () { /* FE constants are defined only where supported, in particular for instance some ARM systems have been seen with only a couple of modes */ static const int modes[] = { 0, #ifdef FE_TONEAREST FE_TONEAREST, #endif #ifdef FE_UPWARD FE_UPWARD, #endif #ifdef FE_DOWNWARD FE_DOWNWARD, #endif #ifdef FE_TOWARDZERO FE_TOWARDZERO, #endif }; double x, want; int i; for (i = 0; i < numberof (modes); i++) { /* First iteration is the default rounding mode, ie. no call to fesetround. Subsequent iterations are the FE modes from the table. */ if (i != 0) { #if HAVE_FESETROUND fesetround (modes[i]); #endif } assert (scm_round (0.0) == 0.0); assert (scm_round (1.0) == 1.0); assert (scm_round (-1.0) == -1.0); assert (scm_round (0.5) == 0.0); assert (scm_round (1.5) == 2.0); assert (scm_round (-1.5) == -2.0); assert (scm_round (2.5) == 2.0); assert (scm_round (-2.5) == -2.0); assert (scm_round (3.5) == 4.0); assert (scm_round (-3.5) == -4.0); /* 2^(DBL_MANT_DIG-1)-1+0.5 */ x = ldexp (1.0, DBL_MANT_DIG - 1) - 1.0 + 0.5; want = ldexp (1.0, DBL_MANT_DIG - 1); assert (scm_round (x) == want); /* -(2^(DBL_MANT_DIG-1)-1+0.5) */ x = - (ldexp (1.0, DBL_MANT_DIG - 1) - 1.0 + 0.5); want = - ldexp (1.0, DBL_MANT_DIG - 1); assert (scm_round (x) == want); /* 2^DBL_MANT_DIG-1 In the past scm_round had incorrectly incremented this value, due to the way that x+0.5 would round upwards (in the usual default nearest-even mode on most systems). */ x = ldexp (1.0, DBL_MANT_DIG) - 1.0; assert (x == floor (x)); /* should be an integer already */ assert (scm_round (x) == x); /* scm_round should return it unchanged */ /* -(2^DBL_MANT_DIG-1) */ x = - (ldexp (1.0, DBL_MANT_DIG) - 1.0); assert (x == floor (x)); /* should be an integer already */ assert (scm_round (x) == x); /* scm_round should return it unchanged */ /* 2^64 */ x = ldexp (1.0, 64); assert (scm_round (x) == x); /* -2^64 In the past scm_round had incorrectely incremented this value in any mode except FE_NEAREST, due to x+0.5 round either up or down to the next representable value (an integer). */ x = - ldexp (1.0, 64); assert (scm_round (x) == x); } } int main (int argc, char *argv[]) { scm_init_guile(); test_scm_round (); return 0; } --=-=-= Content-Type: text/plain; charset="us-ascii" MIME-Version: 1.0 Content-Transfer-Encoding: 7bit Content-Disposition: inline _______________________________________________ Guile-devel mailing list Guile-devel@gnu.org http://mail.gnu.org/mailman/listinfo/guile-devel --=-=-=--