From owner-freebsd-hackers@FreeBSD.ORG Wed Jul 14 00:54:59 2004 Return-Path: Delivered-To: freebsd-hackers@freebsd.org Received: from mx1.FreeBSD.org (mx1.freebsd.org [216.136.204.125]) by hub.freebsd.org (Postfix) with ESMTP id 73B0016A4CE for ; Wed, 14 Jul 2004 00:54:59 +0000 (GMT) Received: from alcanet.com.au (mail2.alcanet.com.au [203.62.196.17]) by mx1.FreeBSD.org (Postfix) with ESMTP id 594BF43D1F for ; Wed, 14 Jul 2004 00:54:56 +0000 (GMT) (envelope-from peter.jeremy@alcatel.com.au) Received: from gsmx07.alcatel.com.au (IDENT:root@localhost.localdomain [127.0.0.1])i6E0sp3l031533; Wed, 14 Jul 2004 10:54:51 +1000 Received: from gsmx07.alcatel.com.au (localhost [127.0.0.1]) i6E0spks036007; Wed, 14 Jul 2004 10:54:51 +1000 (EST) (envelope-from peter.jeremy@alcatel.com.au) Received: (from pjeremy@localhost) by gsmx07.alcatel.com.au (8.12.9p2/8.12.9/Submit) id i6E0sjHc036003; Wed, 14 Jul 2004 10:54:45 +1000 (EST) (envelope-from peter.jeremy@alcatel.com.au) Date: Wed, 14 Jul 2004 10:54:44 +1000 From: Peter Jeremy To: Dmitry Morozovsky Message-ID: <20040714005444.GA35946@gsmx07.alcatel.com.au> Mail-Followup-To: Dmitry Morozovsky , freebsd-hackers@freebsd.org References: <20040711210219.J84500@woozle.rinet.ru> Mime-Version: 1.0 Content-Type: text/plain; charset=us-ascii Content-Disposition: inline In-Reply-To: <20040711210219.J84500@woozle.rinet.ru> User-Agent: Mutt/1.4.2i X-Mailman-Approved-At: Wed, 14 Jul 2004 12:06:33 +0000 cc: freebsd-hackers@freebsd.org Subject: Re: gcc strangeness X-BeenThere: freebsd-hackers@freebsd.org X-Mailman-Version: 2.1.1 Precedence: list List-Id: Technical Discussions relating to FreeBSD List-Unsubscribe: , List-Archive: List-Post: List-Help: List-Subscribe: , X-List-Received-Date: Wed, 14 Jul 2004 00:54:59 -0000 On 2004-Jul-11 21:06:32 +0400, Dmitry Morozovsky wrote: >marck@woozle:/tmp/tsostik> cat x.c >#include >int main () >{ > float a; > for(a=0.01;a<=0.1; a+=0.01) > printf("%f %.3f %d\n", a*100, a*100, (int)(a*100)); >return 0; >} >marck@woozle:/tmp/tsostik> cc x.c >marck@woozle:/tmp/tsostik> ./a.out >1.000000 1.000 0 >2.000000 2.000 1 >3.000000 3.000 2 >4.000000 4.000 3 >5.000000 5.000 5 >6.000000 6.000 6 >7.000000 7.000 7 >8.000000 8.000 7 >9.000000 9.000 8 >9.999999 10.000 9 As various other people have pointed out, this is an artifact of binary floating point arithmetic, rather than gcc. Note that single-precision floating point only gives you 24 binary (or 7.22 decimal) digits of precision. Expecting to accurately display 7 digits may be somewhat hopeful. The correct solution is to use rint(3), rather than casting to int (which is truncation). ceil(3) will solve this specific issue but is just as susceptable to unexpected rounding. Switching to double precision will help but is still not a complete solution. I looked into a similar problem this a few years ago following a complaint that int(0.29 * 100) == 28 in double precision. The following is my analysis of that problem. The underlying problem is that 1/100 is not representable in a finite number of bits. Specifically, 1/100 = 0d0.01 = 0b0.0000[00101000111101011100] = 0x0.0[28F5C] (where the brackets enclose the repeating digits). As an example, 0.29 = 0x0.4[A3D70]. Expressed as a rounded IEEE 754 double precision number, this is 0x3fd28f5c28f5c28f, which underestimates the correct value by ~2e-17 (0.36ULP[1]). When this number is multiplied by 100, the underestimate becomes 0.56ULP, causing the result to round down (ie the result is 1ULP = ~3.4e-15 low). Since conversion-to-integer is defined to round towards zero, this gives an integer value of 28. A good (though not hyper-linked) bibliography on floating point can be found at http://www.cse.nd.edu/~milind/references/floatingpoint. Prof William Kahan's writings (http://www.cs.berkeley.edu/~wkahan/) also make interesting reading. [1] Units in the Last Place. Peter