# Mean of integers overflows - bug or expected behaviour?
**URL:** https://discourse.julialang.org/t/mean-of-integers-overflows-bug-or-expected-behaviour/34767
**Category:** Statistics
**Tags:** integer-overflow
**Created:** [February 17, 2020, 4:54pm UTC](https://discourse.julialang.org/t/mean-of-integers-overflows-bug-or-expected-behaviour/34767 "2020-02-17T16:54:34Z")
**Posts on this page:** 20
**Page:** 3
### Author:  [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [February 18, 2020, 6:07pm UTC](https://discourse.julialang.org/t/mean-of-integers-overflows-bug-or-expected-behaviour/34767/41 "2020-02-18T18:07:39Z")
> [@Mikhail\_Kagalenko](#):
>
> The type-stable division for integers is div().
All methods of `/` I know of (in `Base` & standard libraries) are [type stable](https://docs.julialang.org/en/v1/manual/performance-tips/#Write-%22type-stable%22-functions-1). If you are aware of one that isn’t, that may be a bug.
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [February 18, 2020, 6:17pm UTC](https://discourse.julialang.org/t/mean-of-integers-overflows-bug-or-expected-behaviour/34767/42 "2020-02-18T18:17:58Z")
> [@Tamas\_Papp](#):
>
> as it converts an inherently associative summation problem to a tricky one with floating point.
The error you reported here is a relative error of one ulp, which is quite reasonable. Since `mean` is giving a floating-point output, the normal expectation would be to have some rounding error — the exact answer will generally not be representable.
(Furthermore, `Float64` addition is exact for integers up to 2^53, so you will only see rounding errors for integers so big that you will be risking `Int64` overflow anyway, and 32-bit machines are even worse.)
For algorithms that compute a floating-point result, the normal expectation is to avoid spurious overflow due to intermediate calculations if it is at all practical. For example, we give the correct answer for functions like `hypot` and `norm` for integer inputs:
```julia
julia> abs(typemax(Int) + typemax(Int)*im)
1.3043817825332783e19
julia> using LinearAlgebra
julia> norm(fill(typemax(Int), 1000))
2.916686334356758e20
```
Indeed, `norm` is a classic case where one accepts a performance price to avoid spurious overflow. On my machine we are paying a factor of 3 even for `Float64` inputs:
```julia
julia> unsafenorm(x) = sqrt(sum(abs2, x))
julia> x = rand(10^4);
julia> @btime norm($x);
2.621 μs (0 allocations: 0 bytes)
julia> @btime unsafenorm($x);
923.531 ns (0 allocations: 0 bytes)
```
We should certainly try to make a common function like `mean` as fast as possible, but not at the expense of correctness, _given that we are already returning a floating-point result_. Start with a correct implementation (e.g. [this one](https://github.com/JuliaLang/Statistics.jl/issues/22#issuecomment-587119323)), then make it faster it if we can.
(Honestly, I’m skeptical that many applications will even notice a constant-factor slowdown of `mean`, since I doubt it is the performance-critical step in almost any nontrivial problem.)
---
### Author:  [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [February 18, 2020, 6:36pm UTC](https://discourse.julialang.org/t/mean-of-integers-overflows-bug-or-expected-behaviour/34767/43 "2020-02-18T18:36:40Z")
> [@stevengj](#):
>
> not at the expense of correctness, _given that we are already returning a floating-point result_ .
I don’t mean to be pedantic here, but `mean` is “correct” in the sense of being equal to `sum(v) / length(v)`, where `sum` is “correct” in the sense of being equal to `reduce(+, v)`.
These of course fail for overflow because Julia chose [native integer arithmetic](https://docs.julialang.org/en/v1/manual/faq/#faq-integer-arithmetic-1). This was a decision with various trade-offs in mind, but it has pervasive implications for the language.
I think your suggestion is a reasonable alternative, but the current implementation is no more “incorrect” than `+(::Int, ::Int)` is. It is not at all obvious to me that mathematical correctness should always be the goal at the expense of speed in Julia.
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [February 18, 2020, 6:44pm UTC](https://discourse.julialang.org/t/mean-of-integers-overflows-bug-or-expected-behaviour/34767/44 "2020-02-18T18:44:01Z")
> [@Tamas\_Papp](#):
>
> but `mean` is “correct” in the sense of being equal to `sum(v) / length(v)`
This seems like a circular argument; that’s an algorithm, and you could _assert_ it as a definition, but I would say that it is not the _conventional_ definition. You could just as well “define” `norm(x) = sqrt(sum(abs2, x))` and say that `Inf` is therefore the “correct” result for `norm([1e300])`. Put another way, if we change the implementation of `sum` in a way that changes the roundoff error for floating-point inputs, are we changing the “definition” of `mean(::Vector{Float64})`?
**For a floating-point function** , it is the ordinary expectation to define “correctness” and “accuracy” with respect to the function as computed with _infinite_ precision (i.e. over ℝ). I think that expectation applies here, since `mean(::Vector{Int})` returns a `Float64`.
(Matters are different for arithmetic functions whose inputs _and_ outputs are a finite-precision integer type.)
---
### Author:  [@Palli](https://discourse.julialang.org/u/Palli)
#### Post date: [February 18, 2020, 7:27pm UTC](https://discourse.julialang.org/t/mean-of-integers-overflows-bug-or-expected-behaviour/34767/45 "2020-02-18T19:27:21Z")
> [@stillyslalom](#):
>
> The conversion operation ( `vcvtsi2sdq` on my machine) is expensive compared to `vaddsd` , which leads to a ~6x slowdown for `Int64` → `Float64` conversion
I’m fairly confident we can avoid slowdown, with tricks like Per’s. Just as with division, mean changes to floats, and we can and should define the operation to work without overflow. I could look into the performance issue if it turns out to be real.
Off-topic:
> [@Mikhail\_Kagalenko](#):
>
> Outputs of analog-to-digital converters are integer, and they are **the** most typical measurement data.
Yes, I think so, and for DAC, while I recall on the Acorn Archimedes, 8-bit sound was logarithmic, making it sound better (and volume control simpler, just subtraction, fast on the ARM with no DIV back then). Are logarithmic DAC/ADC outdated now?
---
### Author:  [@Elrod](https://discourse.julialang.org/u/Elrod)
#### Post date: [February 18, 2020, 7:42pm UTC](https://discourse.julialang.org/t/mean-of-integers-overflows-bug-or-expected-behaviour/34767/46 "2020-02-18T19:42:07Z")
> [@Palli](#):
>
> I’m fairly confident we can avoid slowdown, with trick’s like Per’s. Just as with division, mean changes to floats, and we can and should define the operation to work without overflow. I could look into the performance issue if it turns out to be real.
I’m sure it’s real, and doubt most of the slowdown is avoidable:
```julia
julia> function mysum1(x)
s = 0
for xᵢ ∈ x
s += xᵢ
end
s
end
mysum1 (generic function with 1 method)
julia> function mysum2(x)
s = zero(Int128)
for xᵢ ∈ x
s += xᵢ
end
s % Int64
end
mysum2 (generic function with 1 method)
julia> function mysum3(x)
s = 0.0
@fastmath for xᵢ ∈ x
s += xᵢ
end
s
end
mysum3 (generic function with 1 method)
julia> using BenchmarkTools
julia> x = rand(-100:100, 200);
julia> @btime mysum1($x)
17.941 ns (0 allocations: 0 bytes)
-115
julia> @btime mysum2($x)
131.092 ns (0 allocations: 0 bytes)
-115
julia> @btime mysum3($x)
158.761 ns (0 allocations: 0 bytes)
-115.0
julia> @btime mysum1($x)
17.935 ns (0 allocations: 0 bytes)
-115
```
This is much better with AVX512, where the conversion to `Float64` can be vectorized: (testing on AWS):
```julia
julia> function mysum1(x)
s = 0
for xᵢ ∈ x
s += xᵢ
end
s
end
mysum1 (generic function with 1 method)
julia> function mysum2(x)
s = zero(Int128)
for xᵢ ∈ x
s += xᵢ
end
s % Int64
end
mysum2 (generic function with 1 method)
julia> function mysum3(x)
s = 0.0
@fastmath for xᵢ ∈ x
s += xᵢ
end
s
end
mysum3 (generic function with 1 method)
julia> using BenchmarkTools
julia> x = rand(-100:100, 200);
julia> @btime mysum1($x)
14.732 ns (0 allocations: 0 bytes)
1008
julia> @btime mysum2($x)
227.857 ns (0 allocations: 0 bytes)
1008
julia> @btime mysum3($x)
38.964 ns (0 allocations: 0 bytes)
1008.0
julia> @btime mysum1($x)
14.715 ns (0 allocations: 0 bytes)
1008
```
That’s nearly a 9x penalty with AVX2, and about a 2.5x penalty with AVX512.
Whether more users are likely to complain about overflowing integer means, or a constant-factor performance hit is a valid question.
On the AVX2 computer, we can improve performance (over the 9x penalty) via:
```julia
julia> mysum4(x) = maximum(x)*Float64(length(x)) < typemax(Int) ? Float64(mysum1(x)) : mysum3(x)
mysum4 (generic function with 1 method)
julia> @btime mysum4($x)
69.165 ns (0 allocations: 0 bytes)
489.0
```
More elaborate techniques following this approach will probably be able to improve things further.
`mysum4` is just a smidge slower than `mysum3` on the AVX512 computer, so a slight imrpovement is all it would take to be faster than `mysum3` on both systems.
However, given longer vectors where the summation becomes memory bound, `mysum3` may approach 2x faster than `mysum4`.
EDIT:
I just realized it should be `maximum(abs, x)`, which is much slower:
```julia
julia> mysum4(x) = maximum(abs,x)*Float64(length(x)) < typemax(Int) ? Float64(mysum1(x)) : mysum3(x)
mysum4 (generic function with 1 method)
julia> @btime mysum4($x)
93.891 ns (0 allocations: 0 bytes)
489.0
```
---
### Author:  [@StefanKarpinski](https://discourse.julialang.org/u/StefanKarpinski)
#### Post date: [February 29, 2020, 4:57pm UTC](https://discourse.julialang.org/t/mean-of-integers-overflows-bug-or-expected-behaviour/34767/47 "2020-02-29T16:57:55Z")
Converting to float first is not strictly more correct. As has been shown there are cases where integers overflow, which leads to one kind error, but there as also cases where floats give the wrong answer because they lose precision and are not associative. It would be very surprising to me if computing `mean(v)` did not behave the same as `sum(v)/length(v)`. That means when the values are integers, I expect overflow as a possibility and do not expect lack of associativity.
---
### Author:  [@Mikhail\_Kagalenko](https://discourse.julialang.org/u/Mikhail_Kagalenko)
#### Post date: [February 29, 2020, 5:14pm UTC](https://discourse.julialang.org/t/mean-of-integers-overflows-bug-or-expected-behaviour/34767/48 "2020-02-29T17:14:25Z")
> Converting to float first is not strictly more correct. As has been shown there are cases where integers overflow, which leads to one kind error, but there as also cases where floats give the wrong answer because they lose precision and are not associative.
The former kind of error is catastrophic, that is, the answer has no valid digits. The later isn’t - the loss of accuracy in a floating-point sum is not unexpected.
> t would be very surprising to me if computing `mean(v)` did not behave the same as `sum(v)/length(v)` .
That looks like an argument in favour of making the result of mean() on integers to be an integer as well. That might be actually an improvement over the present situation.
---
### Author:  [@StefanKarpinski](https://discourse.julialang.org/u/StefanKarpinski)
#### Post date: [February 29, 2020, 5:27pm UTC](https://discourse.julialang.org/t/mean-of-integers-overflows-bug-or-expected-behaviour/34767/49 "2020-02-29T17:27:28Z")
It’s not because division of integers does not produce an integer.
---
### Author:  [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)
#### Post date: [February 29, 2020, 6:02pm UTC](https://discourse.julialang.org/t/mean-of-integers-overflows-bug-or-expected-behaviour/34767/50 "2020-02-29T18:02:31Z")
> [@Mikhail\_Kagalenko](#):
>
> the loss of accuracy in a floating-point sum is not unexpected
The point that some of us are trying to make here is that one may not expect the _sum of integers_ to be performed in floating point.
I fully understand why some users would prefer that `mean` protects them from “surprises” like that. However, the right question in Julia is usually not what you can improve for a very specific set of argument types, but what the API can/should promise for _generic_ types.
For example, `mean` converting an `AbstractVector{Int}` to `Float64` elementwise before doing the calculation is pretty useless unless we guarantee that something similar also happens for all subtypes of `Integer`, because unless we do this, it is just a special case that does not help generic code so it is hard to rely on.
---
### Author:  [@aplavin](https://discourse.julialang.org/u/aplavin)
#### Post date: [February 29, 2020, 6:08pm UTC](https://discourse.julialang.org/t/mean-of-integers-overflows-bug-or-expected-behaviour/34767/51 "2020-02-29T18:08:16Z")
A really nice general solution would be something like propagating @fastmath to all functions called. So the user can choose either the correct & slower computation, or almost correct & much faster.
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [February 29, 2020, 6:53pm UTC](https://discourse.julialang.org/t/mean-of-integers-overflows-bug-or-expected-behaviour/34767/52 "2020-02-29T18:53:53Z")
> [@Tamas\_Papp](#):
>
> For example, `mean` converting an `AbstractVector{Int}` to `Float64` elementwise before doing the calculation is pretty useless unless we guarantee that something similar also happens for all subtypes of `Integer` , because unless we do this, it is just a special case that does not help generic code so it is hard to rely on.
My suggested solution does just that in a type-generic way: [https://github.com/JuliaLang/Statistics.jl/issues/22#issuecomment-587119323](https://github.com/JuliaLang/Statistics.jl/issues/22#issuecomment-587119323)
I don’t think the principle here is very difficult to understand: if the outputs are floating point, then you want a computation that is numerically stable and reasonably accurate in the floating-point sense, even if the inputs are integers. (Example: the `norm` of an integer array.) The current `mean` algorithm, which can give a large error for integer inputs regardless of the floating-point precision, is numerically unstable according to the formal definition.
---
### Author:  [@Mikhail\_Kagalenko](https://discourse.julialang.org/u/Mikhail_Kagalenko)
#### Post date: [February 29, 2020, 6:53pm UTC](https://discourse.julialang.org/t/mean-of-integers-overflows-bug-or-expected-behaviour/34767/53 "2020-02-29T18:53:54Z")
> The point that some of us are trying to make here is that one may not expect the _sum of integers_ to be performed in floating point.
We are not discussing computation of the sum of integers, though.
---
### Author:  [@Mikhail\_Kagalenko](https://discourse.julialang.org/u/Mikhail_Kagalenko)
#### Post date: [February 29, 2020, 7:00pm UTC](https://discourse.julialang.org/t/mean-of-integers-overflows-bug-or-expected-behaviour/34767/54 "2020-02-29T19:00:17Z")
> It’s not because division of integers does not produce an integer.
Sometimes it does:
```julia
julia> typeof(div(1,2))
Int32
```
---
### Author:  [@StefanKarpinski](https://discourse.julialang.org/u/StefanKarpinski)
#### Post date: [February 29, 2020, 7:00pm UTC](https://discourse.julialang.org/t/mean-of-integers-overflows-bug-or-expected-behaviour/34767/55 "2020-02-29T19:00:56Z")
Sure I guess if you call a totally different function?
---
### Author:  [@Mikhail\_Kagalenko](https://discourse.julialang.org/u/Mikhail_Kagalenko)
#### Post date: [February 29, 2020, 7:02pm UTC](https://discourse.julialang.org/t/mean-of-integers-overflows-bug-or-expected-behaviour/34767/56 "2020-02-29T19:02:06Z")
It’s a “division of integers” function. Another one is `//`
---
### Author:  [@StefanKarpinski](https://discourse.julialang.org/u/StefanKarpinski)
#### Post date: [February 29, 2020, 7:14pm UTC](https://discourse.julialang.org/t/mean-of-integers-overflows-bug-or-expected-behaviour/34767/57 "2020-02-29T19:14:50Z")
Yes, there are many forms of integer division. I don’t really get what that shows. The only division operation that generally computes the mean in this case is the one that produces a float. `sum(v)/length(v)` is also the most obvious definition. Deviation from that, e.g. by converting to float first, could be opted into, but my point is that while you’re surprised that this overflows, other people (like myself) would be surprised that it suffers from the floating point issues given that it’s apparently an operation on integers. It’s not just the loss of precision, it’s the fact that for integers answers are perfectly reproducible, whereas for floats, they’re not—they depend on the architecture of the machine, the details of the summation algorithm and the number of threads.
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [February 29, 2020, 7:46pm UTC](https://discourse.julialang.org/t/mean-of-integers-overflows-bug-or-expected-behaviour/34767/58 "2020-02-29T19:46:14Z")
> [@StefanKarpinski](#):
>
> other people (like myself) would be surprised that it suffers from the floating point issues given that it’s apparently an operation on integers
It gives a floating-point answer, so I don’t understand why you’d be surprised that the algorithm might involve floating-point calculations, or that it would be designed to give an accurate answer in the floating-point sense (as opposed to the 2s-complement sense).
Are you unpleasantly surprised that `norm(x::Vector{Int})` does not give the same answer as `sqrt(sum(abs2, x))`? Would you rather have `norm([10^13])` give a `DomainError` (as the naive algorithm would produce)?
---
### Author:  [@StefanKarpinski](https://discourse.julialang.org/u/StefanKarpinski)
#### Post date: [February 29, 2020, 7:59pm UTC](https://discourse.julialang.org/t/mean-of-integers-overflows-bug-or-expected-behaviour/34767/59 "2020-02-29T19:59:35Z")
No, of course or would be good to give the numerically correct answer more of the time. My point is that what is a surprising result depends on your background and expectations. Doing integer summation followed by float division is obvious and fast in this case. If we can figure out a way to avoid overflow while preserving performance and accuracy, that would be ideal.
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [February 29, 2020, 8:03pm UTC](https://discourse.julialang.org/t/mean-of-integers-overflows-bug-or-expected-behaviour/34767/60 "2020-02-29T20:03:20Z")
The naive algorithm for `norm` is almost 20× faster than the algorithm we are currently using:
```julia
julia> x = rand(1:100, 1000);
julia> naivenorm(x) = sqrt(sum(abs2, x))
naivenorm (generic function with 1 method)
julia> @btime norm($x); @btime naivenorm($x);
3.397 μs (0 allocations: 0 bytes)
189.975 ns (0 allocations: 0 bytes)
```
We made the tradeoff here because it was more important to give an accurate answer that avoids spurious overflow catastrophies, and I think the same consideration applies to `mean`.
(Neither `norm` nor `mean` seems likely to be the performance-critical step in real applications, which is an additional reason to prefer accuracy and robustness to performance.)
[Previous page](https://discourse.julialang.org/t/mean-of-integers-overflows-bug-or-expected-behaviour/34767.md?page=2)
[Next page](https://discourse.julialang.org/t/mean-of-integers-overflows-bug-or-expected-behaviour/34767.md?page=4)