round() wrong in Python 2.4?

[Magnus Lycka]

Nils Grimsmo wrote:
> (Is this due to the different GCC used?)

Yes, but there are probably other nasty values with the
other CGG. Basically, what the code does, for a positive
number, is to calculate floor(0.0225*1000.0+0.5)/1000.0.

As others have said: Don't trust this. If you use Python 2.4,
you can take advantage of the new decimal module, where no
floating point calculations are involved.

I tried with Python 2.2.3 on RH EL3, and for half a million
tested values that ended with 5 in the 4th decimal and used
ndigits=3, I got almost 3000 where it rounded towards zero,
and not towards infinity as the docs say. I.e. almost 0.6%
wrong.

Here is a more direct description of the problem in my system:
 >>> math.floor(4.0925*1.0*10.0*10.0*10.0+0.5)
4093.0
 >>> math.floor(4.0935*1.0*10.0*10.0*10.0+0.5)
4093.0
 >>> math.floor(4.0945*1.0*10.0*10.0*10.0+0.5)
4095.0

Your 2.4 system is still strange though. I tried the
program below on a range of systems: RH Linux, HP-UX,
AIX, Solaris, on Sparc, PowerPC, PA-RISC, Intel Pentium
and AMD 64, and they always gave the same results with
Python 2.2.3 or Python 2.3.1.

Program:

for N in (1000,2000,3000,5000,10000,100000,1000000):
     buggy=0
     for i in range(1,N,2):
         f=i/2000.0
         r=round(f,3)
         if r<f:
             buggy+=1
     print "%7i %10f %5i %f%%" % (N/2,f,buggy,buggy*200./N)

Consistent output:

     500   0.499500     0 0.000000%
    1000   0.999500    12 1.200000%
    1500   1.499500    12 0.800000%
    2500   2.499500    24 0.960000%
    5000   4.999500    47 0.940000%
   50000  49.999500   369 0.738000%
  500000 499.999500  2950 0.590000%

So, while N*1000.0 + 0.5 might sometimes be a little less than
an integer, even though N is an odd integer divided with 2000.0,
it seems that machines handling IEEE floating point numbers
agree about which numbers are affected, and 0.0225 should not
be a problem number:
 >>> round(0.0225,3)
0.023

There have been problems with GCC's float() before though...
http://lists.debian.org/debian-gcc/2002/04/msg00056.html

> How do you correctly output floating-point numbers in 2.4?

There is no change here.
0.023 => 0.023 and 0.022 => 0.021999999999999999
in different Python versions. Use str() or %s etc.

BTW, the C source code looks like this:

static PyObject *
builtin_round(PyObject *self, PyObject *args)
{
	double x;
	double f;
	int ndigits = 0;
	int i;

	if (!PyArg_ParseTuple(args, "d|i:round", &x, &ndigits))
			return NULL;
	f = 1.0;
	i = abs(ndigits);
	while  (--i >= 0)
		f = f*10.0;
	if (ndigits < 0)
		x /= f;
	else
		x *= f;
	if (x >= 0.0)
		x = floor(x + 0.5);
	else
		x = ceil(x - 0.5);
	if (ndigits < 0)
		x *= f;
	else
		x /= f;
	return PyFloat_FromDouble(x);
}

Perhaps one could argue that the code should be changed to
	if (x >= 0.0)
		x = floor(x + d + 0.5);
	else
		x = ceil(x - d - 0.5);
where d is a fairly small number, but this doesn't help in
the long run... For large enough floating point numbers,
the resolution of the floating point system gets bigger than
1! It might well be possible to make a round() function that
works just right in e.g. business accounting applications, where
money ranges between perhaps 0.01 and 1,000,000,000,000.00, but
it's much more difficult to make such a thing work for the
standard library, where we might want to use the whole range
available to floats.