Solve the overflow in mean() on integers by promoting accumulator

src/Statistics.jl

@@ -167,7 +167,16 @@ julia> mean(A, dims=2)
 mean(A::AbstractArray; dims=:) = _mean(A, dims)
 
 _mean(A::AbstractArray{T}, region) where {T} = mean!(Base.reducedim_init(t -> t/2, +, A, region), A)
-_mean(A::AbstractArray, ::Colon) = sum(A) / length(A)
+function _mean(A::AbstractArray, ::Colon)
+    isempty(A) && return sum(A)/0
+    n = length(A)
+    x1 = first(A) / n
+    _prom(x::T, y::S) where {T,S} = begin
+        R = promote_type(T, S)
+        return convert(R, x)
+    end
+    return sum(x->_prom(x,x1), A) / n
+end
 
 function mean(r::AbstractRange{<:Real})
     isempty(r) && return oftype((first(r) + last(r)) / 2, NaN)

test/runtests.jl

@@ -130,6 +130,20 @@ end
         @test mean(identity, x) == mean(identity, g) == typemax(T)
         @test mean(x, dims=2) == [typemax(T)]'
     end
+    # Check that mean of integers does not cause catastrophic loss of accuracy
+    let x = fill(typemax(Int), 10)
+        @test (mean(x) == mean(x, dims=1)[] == mean(float, x)
+               ≈ float(typemax(Int)))  # avoid integer overflow (#22)
+    end
+    let x = rand(10000)  # mean should use sum's accurate pairwise algorithm
+        @test mean(x) == sum(x) / length(x)
+    end
+    @test mean(Number[1, 1.5, 2+3im]) === 1.5+1im # mixed-type array
+    @test isequal(mean(Float64[]), NaN)
+    @test isequal(mean(Int[]), NaN)
+    @inferred mean(Int[])
+    @inferred mean(Float32[])
+    @test isequal(typeof(mean(Float32[])), typeof(mean(Float32[1])))
 end
 
 @testset "mean/median for ranges" begin
@@ -710,7 +724,7 @@ end
     x = Any[1, 2, 4, 10]
     y = Any[1, 2, 4, 10//1]
     @test var(x) === 16.25
-    @test var(y) === 65//4
+    @test var(y) === 16.25
     @test std(x) === sqrt(16.25)
     @test quantile(x, 0.5)  === 3.0
     @test quantile(x, 1//2) === 3//1