# Julia equivalent of Python's "fsum" for floating point summation

**URL:** <https://discourse.julialang.org/t/julia-equivalent-of-pythons-fsum-for-floating-point-summation/17785>\
**Category:** General Usage\
**Tags:** python\
**Created:** [November 20, 2018, 5:16pm UTC](https://discourse.julialang.org/t/julia-equivalent-of-pythons-fsum-for-floating-point-summation/17785 "2018-11-20T17:16:28Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![anon32405968](https://avatars.discourse-cdn.com/v4/letter/a/57b2e6/32.png) [@anon32405968](https://discourse.julialang.org/u/anon32405968)\
**Post date:** [November 20, 2018, 5:16pm UTC](https://discourse.julialang.org/t/julia-equivalent-of-pythons-fsum-for-floating-point-summation/17785/1 "2018-11-20T17:16:28Z")

</div>

I have been looking for a Julia equivalent of Python’s `fsum` and have thus far been unable to find one. Is there an implementation already available? If not I will write my own and put it on Github for others to use.

**From the Python docs:**

`math.fsum` ( _iterable_ )

_Return an accurate floating point sum of values in the iterable. Avoids loss of precision by tracking multiple intermediate partial sums:_

_The algorithm’s accuracy depends on IEEE-754 arithmetic guarantees and the typical case where the rounding mode is half-even. On some non-Windows builds, the underlying C library uses extended precision addition and may occasionally double-round an intermediate sum causing it to be off in its least significant bit._

_For further discussion and two alternative approaches, see the [ASPN cookbook recipes for accurate floating point summation](https://code.activestate.com/recipes/393090/)._

---

<div class="post-metadata">

**Author:** ![stillyslalom](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stillyslalom/32/45687_2.png) [@stillyslalom](https://discourse.julialang.org/u/stillyslalom)\
**Post date:** [November 20, 2018, 5:26pm UTC](https://discourse.julialang.org/t/julia-equivalent-of-pythons-fsum-for-floating-point-summation/17785/2 "2018-11-20T17:26:17Z")

</div>

There was once a Kahan summation function in Base, but it got split off to a package prior to 1.0: [GitHub - JuliaMath/KahanSummation.jl: Sum and cumulative sum using the Kahan-Babuska-Neumaier algorithm](https://github.com/JuliaMath/KahanSummation.jl)

---

<div class="post-metadata">

**Author:** ![Tero\_Frondelius](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tero_frondelius/32/7629_2.png) [@Tero\_Frondelius](https://discourse.julialang.org/u/Tero_Frondelius)\
**Post date:** [November 20, 2018, 5:26pm UTC](https://discourse.julialang.org/t/julia-equivalent-of-pythons-fsum-for-floating-point-summation/17785/3 "2018-11-20T17:26:24Z")

</div>

Maybe this is relevant: [https://github.com/dpsanders/IntervalArithmetic.jl](https://github.com/dpsanders/IntervalArithmetic.jl)

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [November 20, 2018, 6:15pm UTC](https://discourse.julialang.org/t/julia-equivalent-of-pythons-fsum-for-floating-point-summation/17785/4 "2018-11-20T18:15:55Z")

</div>

My understanding is that Python’s `fsum` function doesn’t use Kahan summation, it uses an [algorithm by Shewchuk](http://code.activestate.com/recipes/393090/) that guarantees an exactly rounded result (i.e. the result as if you did the sum in infinite precision and then rounded it to the closest floating-point value).

Kahan summation, in contrast, doesn’t guarantee an exactly rounded result, but its error (relative to the L₁ norm of the input) is bounded independent of the number of terms. The default `sum` function in Julia is also quite accurate because it uses [pairwise summation](https://en.wikipedia.org/wiki/Pairwise_summation), with an error bound that grows only as the log of the number of terms.

(My understanding is that Shewchuk’s algorithm is especially useful if you want a small relative error even for ill-conditioned sums: sums of both positive and negative values where the result is nearly zero due to cancellations. I’ve never used the algorithm myself, though.)

---

<div class="post-metadata">

**Author:** ![StefanKarpinski](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stefankarpinski/32/24_2.png) [@StefanKarpinski](https://discourse.julialang.org/u/StefanKarpinski)\
**Post date:** [November 20, 2018, 8:19pm UTC](https://discourse.julialang.org/t/julia-equivalent-of-pythons-fsum-for-floating-point-summation/17785/5 "2018-11-20T20:19:15Z")

</div>

See also Radford Neal’s paper on exact summation with superaccumulators:

> **[Fast exact summation using small and large superaccumulators](https://arxiv.org/abs/1505.05571)**
>
> I present two new methods for exactly summing a set of floating-point numbers, and then correctly rounding to the nearest floating-point number. Higher accuracy than simple summation (rounding after each addition) is important in many applications,...

If someone wants to implement exact summation, it might be worth comparing different approaches.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [November 20, 2018, 9:30pm UTC](https://discourse.julialang.org/t/julia-equivalent-of-pythons-fsum-for-floating-point-summation/17785/6 "2018-11-20T21:30:40Z")

</div>

Someone may want to contact the author to ask for the C code under an MIT or BSD license, as otherwise translating the code from the paper may be problematic.

---

<div class="post-metadata">

**Author:** ![StefanKarpinski](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stefankarpinski/32/24_2.png) [@StefanKarpinski](https://discourse.julialang.org/u/StefanKarpinski)\
**Post date:** [November 20, 2018, 10:19pm UTC](https://discourse.julialang.org/t/julia-equivalent-of-pythons-fsum-for-floating-point-summation/17785/7 "2018-11-20T22:19:18Z")

</div>

Or implement the algorithms without copying the code.

---

<div class="post-metadata">

**Author:** ![JeffreySarnoff](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jeffreysarnoff/32/1980_2.png) [@JeffreySarnoff](https://discourse.julialang.org/u/JeffreySarnoff)\
**Post date:** [November 21, 2018, 5:56am UTC](https://discourse.julialang.org/t/julia-equivalent-of-pythons-fsum-for-floating-point-summation/17785/8 "2018-11-21T05:56:02Z")

</div>

I wrote Radford, and asked about use of the C code under the MIT or BSD License.

---

<div class="post-metadata">

**Author:** ![Sukera](https://avatars.discourse-cdn.com/v4/letter/s/ce7236/32.png) [@Sukera](https://discourse.julialang.org/u/Sukera)\
**Post date:** [November 21, 2018, 6:32am UTC](https://discourse.julialang.org/t/julia-equivalent-of-pythons-fsum-for-floating-point-summation/17785/9 "2018-11-21T06:32:54Z")

</div>

There’s also the possibility of one person transpiling the code from the paper to pseudo code and someone else implementing the algorithm based solely on the pseudo code (without reading the original implementation!) - that way, no explicit copying is going on, which is the part that would be problematic.

I’ve not read the paper yet, so I’d be happy to do the second part sometime next week if someone can provide a pseudo code version of this.

---

<div class="post-metadata">

**Author:** ![StefanKarpinski](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stefankarpinski/32/24_2.png) [@StefanKarpinski](https://discourse.julialang.org/u/StefanKarpinski)\
**Post date:** [November 21, 2018, 2:52pm UTC](https://discourse.julialang.org/t/julia-equivalent-of-pythons-fsum-for-floating-point-summation/17785/10 "2018-11-21T14:52:19Z")

</div>

That’s the most paranoid version of GPL avoidance, but it’s not really legally necessary. As long as you don’t derive your code from GPL code, you can license your own code however you want. Of course the question is what counts as “deriving” your code from GPL code? Do you have to explicitly copy or translate it? Or can you be subconsciously influenced by GPL code you’ve seen but not consciously copied?

---

<div class="post-metadata">

**Author:** ![Sukera](https://avatars.discourse-cdn.com/v4/letter/s/ce7236/32.png) [@Sukera](https://discourse.julialang.org/u/Sukera)\
**Post date:** [November 21, 2018, 2:58pm UTC](https://discourse.julialang.org/t/julia-equivalent-of-pythons-fsum-for-floating-point-summation/17785/11 "2018-11-21T14:58:09Z")

</div>

Not my idea, I’ve only described the technique as told [here](https://lobste.rs/s/q94fty/what_does_gpl_say_about_translating_some#c_xwg1bs) - as far as I understand, any derivation of GPL’d code has to be GPL’d itself. By describing the algorithm with various sources without directly translating it to make it as easy as possible to reimplement it, it’s possible to (somewhat) circumvent this.

IANAL though, so this might just be bonkers 😃

---

<div class="post-metadata">

**Author:** ![anon32405968](https://avatars.discourse-cdn.com/v4/letter/a/57b2e6/32.png) [@anon32405968](https://discourse.julialang.org/u/anon32405968)\
**Post date:** [November 21, 2018, 3:49pm UTC](https://discourse.julialang.org/t/julia-equivalent-of-pythons-fsum-for-floating-point-summation/17785/12 "2018-11-21T15:49:46Z")

</div>

Thanks for the feedback everyone — looks like I have plenty of alternatives to choose from at the moment.

> The default `sum` function in Julia is also quite accurate because it uses [pairwise summation](https://en.wikipedia.org/wiki/Pairwise_summation)

I didn’t realize the regular `sum` function uses pairwise summation, this will probably be adequate for my present requirements. I need to hand in a PhD thesis in the 2 months time so it sounds like someone else will get to implementing Shewchuck’s (or Neal’s) algorithm before me. If not I will look at it again at that time.

**EDIT** : Out of curiosity, would this `fsum` equivalent belong in the stdlib or somewhere else?

---

<div class="post-metadata">

**Author:** ![StefanKarpinski](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stefankarpinski/32/24_2.png) [@StefanKarpinski](https://discourse.julialang.org/u/StefanKarpinski)\
**Post date:** [November 21, 2018, 4:36pm UTC](https://discourse.julialang.org/t/julia-equivalent-of-pythons-fsum-for-floating-point-summation/17785/13 "2018-11-21T16:36:18Z")

</div>

> [@Sukera](#):
>
> Not my idea, I’ve only described the technique as told [here](https://lobste.rs/s/q94fty/what_does_gpl_say_about_translating_some#c_xwg1bs) - as far as I understand, any derivation of GPL’d code has to be GPL’d itself

Please do not do that. Translating GPL code to pseudocode that “looks suspiciously like Python” and then having someone else translate that pseudocode into actual Python (or in this case Julia) is far more legally dicey than just reading the paper and implementing the algorithm based on the description while disregarding the GPL code that’s included. Translating twice does not break the “derived work link”—you can translate through as many different programming languages as you want and it’s still a derived work.

---

<div class="post-metadata">

**Author:** ![StefanKarpinski](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stefankarpinski/32/24_2.png) [@StefanKarpinski](https://discourse.julialang.org/u/StefanKarpinski)\
**Post date:** [November 21, 2018, 6:04pm UTC](https://discourse.julialang.org/t/julia-equivalent-of-pythons-fsum-for-floating-point-summation/17785/14 "2018-11-21T18:04:01Z")

</div>

> [@anon32405968](#):
>
> I need to hand in a PhD thesis in the 2 months time so it sounds like someone else will get to implementing Shewchuck’s (or Neal’s) algorithm before me.

What better time for procrastination? 😬

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [November 21, 2018, 10:29pm UTC](https://discourse.julialang.org/t/julia-equivalent-of-pythons-fsum-for-floating-point-summation/17785/15 "2018-11-21T22:29:48Z")

</div>

(Note that the [published code](https://www.cs.toronto.edu/~radford/xsum.software.html) by Radford is not GPL—in fact, it has no license statement at all, which means that you cannot redistribute it, or code derived from it, in any form. Clean-room re-implementation is not about “GPL avoidance,” it’s about working around copyright law’s definition of “derived work.” In this case, the paper itself contains code excerpts, so it is especially tricky to re-implement the algorithm in a clean way from a copyright perspective.)

The author emailed me to say he was planning on posting a cleaned-up version of the code under a free/open-source license (his initial impulse was GPL, but I’ve urged him to consider a more permissive license for a small library routine like this). Not sure about the timeframe.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [November 22, 2018, 3:37pm UTC](https://discourse.julialang.org/t/julia-equivalent-of-pythons-fsum-for-floating-point-summation/17785/16 "2018-11-22T15:37:53Z")

</div>

5 posts were split to a new topic: [Copyright issues for code excerpts](https://discourse.julialang.org/t/copyright-issues-for-code-excerpts/17868)

---

<div class="post-metadata">

**Author:** ![y4lu](https://avatars.discourse-cdn.com/v4/letter/y/47e85d/32.png) [@y4lu](https://discourse.julialang.org/u/y4lu)\
**Post date:** [November 22, 2018, 9:03am UTC](https://discourse.julialang.org/t/julia-equivalent-of-pythons-fsum-for-floating-point-summation/17785/18 "2018-11-22T09:03:50Z")

</div>

A Kahan / Pairwise polyalgorithm seems like it could still be pretty fast, it might be +15\*3 operations for summing 256 numbers depending on where you switch. Not sure how much extra precision you’d get

---

<div class="post-metadata">

**Author:** ![StefanKarpinski](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stefankarpinski/32/24_2.png) [@StefanKarpinski](https://discourse.julialang.org/u/StefanKarpinski)\
**Post date:** [November 22, 2018, 4:58pm UTC](https://discourse.julialang.org/t/julia-equivalent-of-pythons-fsum-for-floating-point-summation/17785/19 "2018-11-22T16:58:43Z")

</div>

Pairwise summation is the default in Julia; it is just as fast as naive summation and much more accurate. Kahan summation is even more accurate but significantly slower, which is why it is not the default.

---

<div class="post-metadata">

**Author:** ![y4lu](https://avatars.discourse-cdn.com/v4/letter/y/47e85d/32.png) [@y4lu](https://discourse.julialang.org/u/y4lu)\
**Post date:** [November 23, 2018, 7:58am UTC](https://discourse.julialang.org/t/julia-equivalent-of-pythons-fsum-for-floating-point-summation/17785/20 "2018-11-23T07:58:22Z")

</div>

It’s a pretty neat trick too, and definitely a hard default to argue against. I was just noticing that each pass of the pairwise sum makes a partial Kahan sum on the remaining values 50% cheaper

Another approach

```julia
function finiteDif_sum(arr, d= 0.0)
  x = 0.0; y=0.0;
  y += d;
  for i in arr
    x = x + x-y + d + i; ## +y-x -d ?
    y = y + i + y + i - x - d; ## and + x - y - i + d ?
    end;
  [x,y - d]
  end;

```

Add- might’ve gotten the x-y around the wrong way

---

<div class="post-metadata">

**Author:** ![JeffreySarnoff](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jeffreysarnoff/32/1980_2.png) [@JeffreySarnoff](https://discourse.julialang.org/u/JeffreySarnoff)\
**Post date:** [November 26, 2018, 1:14am UTC](https://discourse.julialang.org/t/julia-equivalent-of-pythons-fsum-for-floating-point-summation/17785/21 "2018-11-26T01:14:57Z")

</div>

I have posted a copy of Radford’s paper with all C code removed.  
[https://github.com/JeffreySarnoff/SuperaccumulatorPaperNoC/blob/master/superaccumulator.pdf](https://github.com/JeffreySarnoff/SuperaccumulatorPaperNoC/blob/master/superaccumulator.pdf)

[Next page](https://discourse.julialang.org/t/julia-equivalent-of-pythons-fsum-for-floating-point-summation/17785.md?page=2)