numpy.mean along multiple axis gives wrong result for large arrays

A mean over an array containing only ones should obviously return only ones. However,

>>> import numpy as np
>>> test_array = np.ones((10000000, 4, 15), dtype = np.float32)
>>> print(test_array.mean(axis=(0,1)))
[ 0.4194304  0.4194304  0.4194304  0.4194304  0.4194304  0.4194304
  0.4194304  0.4194304  0.4194304  0.4194304  0.4194304  0.4194304
  0.4194304  0.4194304  0.4194304]

This returns the correct result:

>>> print(test_array.mean(axis=0).mean(axis=0))
[ 1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.]

I guess the reason for this problem is some overflow, because it does not appear when I use a test array with dtype = np.float64. I would expect numpy to either give the correct result or to at least give a warning whenever such an overflow happens.

surprisingly mean along all axis gives the correct result again:

>>> print(test_array.mean(axis=(0,1,2)))
1.0

(I used numpy-1.12.1 with python 3.5 on Ubuntu 16)

The issue provides no file or test entry point. Start by reproducing the float32 mean over the supplied large array and compare tuple-axis reduction with sequential reductions; done means the tuple-axis result is correct or the overflow behavior is explicitly addressed by a regression test.