Add compensated_add for floats

[clarfonthey]

This method implements Kahan summation, which allows pretty good, mostly numerically stable summing of floats without any allocation or sorting, just tracking a single error term.

A "perfect" summation algorithm would likely be much less performant and be implemented as a method of [f32] and [f64], but for cases where people want [0.1; 10].iter().sum() == 1.0, this is adequate.

[TimNN]

@rust-highfive seems to be down, so randomly assigning to a libs team member:

r? @BurntSushi

[rust-highfive]

The job x86_64-gnu-llvm-3.9 of your PR failed on Travis (raw log). Through arcane magic we have determined that the following fragments from the build log may contain information about the problem.

Click to expand the log.

[00:54:13] ........................test [run-pass] run-pass/mir_heavy_promoted.rs has been running for over 60 seconds
[00:54:33] ............................................................................
[00:54:53] ......................................................................................ii............
[00:55:44] ..................................................i.................................................
[00:55:55] ...i.ii..........................test [run-pass] run-pass/saturating-float-casts.rs has been running for over 60 seconds
[00:56:39] ...........iiiiiii..................................................................................
[00:56:59] ....................................................................................................
[00:57:16] ....................................................................................................
[00:57:35] ....................................................................................
---
[01:12:26] ....................................................................................................
[01:12:39] ....................................................................................................
[01:12:56] ....................................................................................................
[01:13:11] ........................i...........................................................................
[01:13:22] ............................F.F.............F..F....................................................
[01:13:46] ....................................................................................................
[01:13:59] ....................................................................................................
[01:14:11] ....................................................................................................
[01:14:23] ....................................................................................................
---
[01:16:43] ...................................i................................................................
[01:16:44] .
[01:16:44] failures:
[01:16:44] 
[01:16:44] ---- num/f32.rs - f32::f32::compensated_add (line 432) stdout ----
[01:16:44]  error[E0599]: no method named `compensated_add` found for type `{float}` in the current scope
[01:16:44]  --> num/f32.rs:438:27
[01:16:44]   |
[01:16:44] 9 |     let (s, c) = comp_sum.compensated_add(x, comp);
[01:16:44] 
[01:16:44] 
[01:16:44] thread 'num/f32.rs - f32::f32::compensated_add (line 432)' panicked at 'couldn't compile the test', librustdoc/test.rs:325:17
[01:16:44] 
[01:16:44] 
[01:16:44] ---- num/f32.rs - f32::f32::compensated_add (line 449) stdout ----
[01:16:44]  error[E0599]: no method named `compensated_add` found for type `{float}` in the current scope
[01:16:44]  --> num/f32.rs:455:23
[01:16:44]   |
[01:16:44] 9 |         |(s, c), x| s.compensated_add(x, c)
[01:16:44] 
[01:16:44] 
[01:16:44] thread 'num/f32.rs - f32::f32::compensated_add (line 449)' panicked at 'couldn't compile the test', librustdoc/test.rs:325:17
[01:16:44] 
[01:16:44] ---- num/f64.rs - f64::f64::compensated_add (line 443) stdout ----
[01:16:44]  error[E0599]: no method named `compensated_add` found for type `{float}` in the current scope
[01:16:44]  --> num/f64.rs:449:27
[01:16:44]   |
[01:16:44] 9 |     let (s, c) = comp_sum.compensated_add(x, comp);
[01:16:44] 
[01:16:44] 
[01:16:44] thread 'num/f64.rs - f64::f64::compensated_add (line 443)' panicked at 'couldn't compile the test', librustdoc/test.rs:325:17
[01:16:44] 
[01:16:44] ---- num/f64.rs - f64::f64::compensated_add (line 460) stdout ----
[01:16:44]  error[E0599]: no method named `compensated_add` found for type `{float}` in the current scope
[01:16:44]  --> num/f64.rs:466:23
[01:16:44]   |
[01:16:44] 9 |         |(s, c), x| s.compensated_add(x, c)
[01:16:44] 
[01:16:44] 
[01:16:44] thread 'num/f64.rs - f64::f64::compensated_add (line 460)' panicked at 'couldn't compile the test', librustdoc/test.rs:325:17
[01:16:44] 
[01:16:44] failures:
[01:16:44] failures:
[01:16:44]     num/f32.rs - f32::f32::compensated_add (line 432)
[01:16:44]     num/f32.rs - f32::f32::compensated_add (line 449)
[01:16:44]     num/f64.rs - f64::f64::compensated_add (line 443)
[01:16:44]     num/f64.rs - f64::f64::compensated_add (line 460)
[01:16:44] test result: FAILED. 1995 passed; 4 failed; 2 ignored; 0 measured; 0 filtered out
[01:16:44] 
[01:16:44] error: test failed, to rerun pass '--doc'
[01:16:44] 
[01:16:44] 
[01:16:44] 
[01:16:44] command did not execute successfully: "/checkout/obj/build/x86_64-unknown-linux-gnu/stage0/bin/cargo" "test" "--target" "x86_64-unknown-linux-gnu" "--release" "--locked" "--color" "always" "--features" "panic-unwind jemalloc backtrace" "--manifest-path" "/checkout/src/libstd/Cargo.toml" "-p" "core" "--" "--quiet"
[01:16:44] 
[01:16:44] 
[01:16:44] failed to run: /checkout/obj/build/bootstrap/debug/bootstrap test
[01:16:44] Build completed unsuccessfully in 0:32:25
[01:16:44] Build completed unsuccessfully in 0:32:25
[01:16:44] Makefile:58: recipe for target 'check' failed
[01:16:44] make: *** [check] Error 1

The command "stamp sh -x -c "$RUN_SCRIPT"" exited with 2.
travis_time:start:0b33c6ff
$ date && (curl -fs --head https://google.com | grep ^Date: | sed 's/Date: //g' || true)

I'm a bot! I can only do what humans tell me to, so if this was not helpful or you have suggestions for improvements, please ping or otherwise contact @TimNN. (Feature Requests)

[scottmcm]

Compensated summation is a great algorithm, but does this really need to be in std? It feels like once more than basic .sum() is needed it makes sense to go to a crate.

[clarfonthey]

@scottmcm normally I'd agree, although based upon the criteria of:

1. It gets rid of a lot of the weirdness of float arithmetic
2. It's nontrivial to remember and implement
3. It can be used to build more complicated, more stable sum algorithms

I'd say that it's worthy of inclusion. I'd draw parallels to this and the Euclidean division/modulo functions which were recently added to libstd: it's nontrivial to implement, gets rid of a "weirdness" of the usual division/modulo, and helps act as a building block for more complicated things.

I'd definitely be willing to do an RFC instead if you'd prefer once before such a thing is merged, though.

[scottmcm]

Hmm, I'll try and tease out why I have a different reaction to this one compared to some others.

• I was convinced that mod_euc (Euclidean modulo rfcs#2169 (comment)) is basically always what you want when the LHS is signed, without needing other things atop it for normal cases, and that there isn't an obvious domain more than the operator itself that makes sense to turn into a full crate
• I like mul_with_carry (widening_mul rfcs#2417) even though it needs other things built atop it since std can plausibly do it more directly using target information or LLVM tricks

Whereas I feel like I'd still need to get out Numerical Recipes to be confident an algorithm is stable even if I was using compensated_add, and the domain of "numerically-stable routines to calculate various things" seems like a coherent crate, including, for example, variance as well.

I'm not on libs, though, so it's not up to me, and I don't know if an RFC would be desired.

[pietroalbini]

Ping from triage @BurntSushi! This PR needs your review.

[BurntSushi]

@pietroalbini I took myself off the PR review rotation because I just don't have the bandwidth to do it unfortunately. :-(

r? @alexcrichton

[alexcrichton]

Does a method like this have precedence for inclusion in the standard libraries of other languages?

[clarfonthey]

@alexcrichton Python offers fsum which performs a much more accurate, but much less performant stable sum. It uses Shewchuck's algorithm, which is basically a combination of this, reordering of arguments depending on absolute value, sorting, and tracking multiple compensation terms.

I went back to look at the original justification for adding fsum to Python, and it appears that people were in favour of having a simple, accurate way of summing floats in the standard library. In that case, simply changing sum and math.fsum was a quick way of making a lot of programs more accurate.

My main reasoning for offering the simple primitive through Kahan summing was that most people are not going to want the performance cut of a complicated algorithm like fsum's, and Kahan summing is usually good enough to get rid of minor weirdness for numbers of the same magnitude. As stated, these kinds of methods can go in other crates.

Offering simply a compensated_add function allows people to carry the compensated term throughout multiple sums, whereas a simple iter.stable_sum() would discard the compensation term after the iteration is finished. This could mean that a lot of people could end up nesting stable_sums like regular addition, which would ultimately lead to huge performance losses for no gain in accuracy. Considering how simple the usage is, I figured that this would be the best solution if it was seen as worthwhile to have a more-numerically-stable summing algorithm in Rust by default.

[alexcrichton]

Hm so given the history with Python, is there a strong reason as to why this should be included in the standard library? It seems like with two possible options the way we'd prefer to go is to have both evolve on crates.io and if one gains enough traction we pull it into libstd?

[pietroalbini]

So, what should we do with this PR? Close it?

[alexcrichton]

Sounds like it!