[Numpy-discussion] Rounding errors in add.reduce?

Sue Giller sag at hydrosphere.com
Thu Dec 6 12:30:19 EST 2001

import MA
I have been using add.reduce on some arrays with single precision 
data in them.  At some point, the reduction seems to be producing 
incorrect values, caused I presume by floating point rounding 
errors.  THe values are correct if the same array is created as 
double precision.

The odd things is that I can do a straight summing of those values 
into long, single or double variable and get correct answers.  Is this 
a bug or what?  If the difference is due to rounding error, I would 
expect the same errors to show up in the cases of summing the 
individual values.

The output from the following code shows the following.  Note that 
the add.reduce from the single precision array is different from all 
the others.  It doesn't matter the rank or size of array, just so sum of 
values gets to certain size.

accumulated sums
	float	75151440.0 	long   75151440
accumulated from raveled array
	float	75151440.0 	double 75151440.0
	single	75150288.0 	double 75151440.0

--- code ---
import MA
# this causes same problem if array is 1d of [amax*bmax*cmax] in 
amax = 4
bmax = 31
cmax = 12
farr = MA.zeros((amax, bmax, cmax), 'f')           # single float
darr = MA.zeros((amax, bmax, cmax), 'd')           # double float
sum = 0.0
lsum = 0
value = 50505			# reducing this can cause all values to agree
for a in range (0, amax):
    for b in range (0, bmax):        
        for c in range (0, cmax):
            farr[a, b, c] = value
            darr[a, b, c] = value            
            sum = sum + value
            lsum = lsum + value

fflat = MA.ravel(farr)
dflat = MA.ravel(darr)
fsum = dsum = 0.0
for value in fflat:
    fsum = fsum + value
for value in dflat:
    dsum = dsum + value
freduce = MA.add.reduce(fflat)
dreduce = MA.add.reduce(dflat)
print "accumulated sums"
print "\tfloat\t", sum, "\tlong  ", lsum
print "accumulated from raveled array"
print "\tfloat\t", fsum, "\tdouble", dsum
print "add.reduce"
print "\tsingle\t", freduce, "\tdouble", dreduce

More information about the NumPy-Discussion mailing list