John Gustafson’s crusade to replace floating point with something better

Re:John Gustafson’s crusade to replace floating point with s... by theamk on Tuesday July 09, 2019 @01:22AM (#20389207)

We have this in our GPU computations, this is very annoying. The errors are usually too small to matter, but they make comparing binary results / reproducible builds impossible. Here is an example I just made:

    $ python3
    >>> s = [10**-x for x in range(16)]
    >>> s
    [1, 0.1, 0.01, 0.001, 0.0001, 1e-05, 1e-06, 1e-07, 1e-08, 1e-09, 1e-10, 1e-11, 1e-12, 1e-13, 1e-14, 1e-15]
    >>> sum(s)
    1.1111111111111112
    >>> sum(reversed(s))
    1.111111111111111

(it uses CPU for simplicity, but GPU has exactly the same math)

[ Parent | on HN ]

Re:John Gustafson’s crusade to replace floating point with s... by FabHK on Tuesday July 09, 2019 @10:36AM (#20392155)

Summation of an array is surprisingly non-trivial. Below a discussion about it with many interesting links (including Kahan summation and pairwise summation) if you want to go down some unexpected rabbit holes.

Somewhere I've seen a little routine (by Stefan Karpinski of Julia fame) that starts with a fixed pre-determined array of numbers, and you give it a desired target number, and it permutes the numbers of the array such that when you sum them up (naively, from left to right) the result is your desired target number [2].

https://discourse.julialang.org/t/accurate-summation-algorit...

EDIT: [2] https://discourse.julialang.org/t/array-ordering-and-naive-s...

[ Parent | on HN ]

Re:John Gustafson’s crusade to replace floating point with s... by dragontamer on Tuesday July 09, 2019 @10:52AM (#20392335)

If you don't care about cancellation error and are only concerned with perfect bitwise / reproducible results, then sort your numbers.

Sorting the numbers is necessary to minimize cancellation error. But absolute minimum cancellation error can only be done sequentially: the parallel version may return different results.

    l = list(s)
    l.sort()
    sum(l)
    l = list(reversed(s))
    l.sort()
    sum(l)

Both return 1.111111111111111 now. Sorting is the "synchronization" step you need to get all your parallel tasks back into a defined order.

[ Parent | on HN ]