[celrod]

Here is a fascinating post by Stefan Karpinski, one of the creators of Julia: https://discourse.julialang.org/t/array-ordering-and-naive-s...

He shows off a function named "sumsto" which takes a single positive double precision floating point number "x" as an argument.

"sumsto" always returns the same vector of 2046 floating point numbers -- just in an order so that naive left-to-right summation returns "x". Just changing the order of that vector lets you sum to almost any positive double precision number.

If you want to run the code, I didn't see where "realmax" was defined, but this works:

realmax() = prevfloat(typemax(Float64))

[StefanKarpinski]

`realmax` got renamed to `floatmax` in 1.0; similarly, `realmin` got renamed to `floatmin` (they're both very much specific to floating point types).

[joe_the_user]

You think that's bad? Theoretical math has "conditionally convergent series'" [1] Sum in order and it can convergence to a given value. Rearrange the terms and any real number is possible.

I have a hunch there could be a bit of relation between these things.

[1] https://en.wikipedia.org/wiki/Conditional_convergence

[StefanKarpinski]

It's an interesting connection. Any finite collection of real values (represented in floating-point or otherwise) has a definite correct sum. If you're constrained to the floating-point to represent that value, you should ideally round the true value to the closest float and return that. Computing this true sum of a sequence of floating-point values is surprisingly hard though. The state of the art seems to be Radford Neal's paper on superaccumulators from 2015:

https://arxiv.org/abs/1505.05571

Fortunately pairwise summation (Julia's default) is fast and fairly hard to trip up, although it certainly is possible to trick it.