spurious overflow in `std`

[stevengj]

julia> std([1e300, -1e300])
Inf

whereas the correct answer is ≈ 1.4142e300, which is representable.

This is the same issue that arises in computing the norm of a vector, for which the naive algorithm overflows/underflows. This solution is simply to scale the array by dividing out maximum(abs, x) (this can be skipped if the maximum absolute value isn't too big or too small).

[stevengj]

numpy.std also gives the wrong answer here:

>>> numpy.std([1e300, -1e300])
/Users/stevenj/.julia/conda/3/lib/python3.10/site-packages/numpy/core/_methods.py:176: RuntimeWarning: overflow encountered in multiply
  x = um.multiply(x, x, out=x)
inf

but for that matter Numpy's norm function is also broken (numpy/numpy#19097):

>>> numpy.linalg.norm([1e300, -1e300])
inf
>>> numpy.linalg.norm(numpy.asarray([1e300, -1e300]))
inf

whereas Julia's is correct:

julia> norm([1e300, -1e300])
1.4142135623730952e300

and Scipy's is norm correct:

>>> scipy.linalg.norm([1e300, -1e300])
1.4142135623730952e+300

but Scipy's tstd is wrong:

>>> scipy.stats.tstd([1e300, -1e300])
/Users/stevenj/.julia/conda/3/lib/python3.10/site-packages/numpy/lib/nanfunctions.py:1746: RuntimeWarning: overflow encountered in multiply
  sqr = np.multiply(arr, arr, out=arr, where=where)
inf

and GNU Octave's std function is also wrong:

octave:1> std([1e300, -1e300])
ans = Inf

(Would be interesting to check Matlab and R.)

[nalimilan]

R is also wrong:

> sd(c(1e300, -1e300))
[1] Inf

So it seems virtually nobody gets this right.