[Numpy-tickets] [NumPy] #924: problem with summing large array of float32

NumPy numpy-tickets@scipy....
Thu Oct 2 23:09:34 CDT 2008

```#924: problem with summing large array of float32
------------------------+---------------------------------------------------
Reporter:  emil        |       Owner:  somebody
Type:  defect      |      Status:  new
Priority:  high        |   Milestone:
Component:  numpy.core  |     Version:  1.0.1
Severity:  critical    |    Keywords:  sum, dot, float32
------------------------+---------------------------------------------------
In medical imaging large arrays of single floats are often used. In
summing over such an
array I discovered the following problem:
a=ones([840,2200,60],float32)
sum(a.flatten()) yields 16777216 (4096^2)   ???

However, sum(sum(sum(a,axis=0),axis=0),axis=0) yields 110880000 the
correct result.
The problem doesn't occur with float64
dot(a.flatten() , a.flatten() ) also gives 16777216

I don't think the problem lies only with numpy. I tried to circumvent the
problem by
writing my own summation in fortran (gfortran) and using f2py. The
following fortran snippet
gave me also incorrect results

Real*4 sum

sum = 0.

do i = 1,840
do j = 1,2200
do k = 1,60
sum = sum +1.0
end do
end do
end do

sum ends up being 16777216 again! the incorrect result!
on the other hand the following works:
Real*4 sum1, sum

sum = 0.

do i = 1,840
sum1 = 0.
do j = 1,2200
do k = 1,60
sum1 = sum1 +1.0
end do
end do
sum = sum+sum1
end do

The same problem seems to come up in matlab.

My computer has amd64's (opteron) and is running suse linux

--
Ticket URL: <http://scipy.org/scipy/numpy/ticket/924>
NumPy <http://projects.scipy.org/scipy/numpy>
The fundamental package needed for scientific computing with Python.
```