Implement numerically stable sum of floats

[cuviper]

From @vks on June 27, 2017 12:29

Here is an implementation based on the (hidden) one in the standard library:

use std::mem;

/// Numerically stable sum.
///
/// See http://www.cs.cmu.edu/~quake-papers/robust-arithmetic.ps for a proof of
/// correctness.
struct Sum {
    partials: Vec<f64>,
}

impl Sum {
    fn new() -> Sum {
        Sum {
            partials: vec![],
        }
    }

    fn add(&mut self, mut x: f64) {
        let mut j = 0;
        // This inner loop applies `hi`/`lo` summation to each
        // partial so that the list of partial sums remains exact.
        for i in 0..self.partials.len() {
            let mut y: f64 = self.partials[i];
            if x.abs() < y.abs() {
                mem::swap(&mut x, &mut y);
            }
            // Rounded `x+y` is stored in `hi` with round-off stored in
            // `lo`. Together `hi+lo` are exactly equal to `x+y`.
            let hi = x + y;
            let lo = y - (hi - x);
            if lo != 0.0 {
                self.partials[j] = lo;
                j += 1;
            }
            x = hi;
        }
        if j >= self.partials.len() {
            self.partials.push(x);
        } else {
            self.partials[j] = x;
            self.partials.truncate(j + 1);
        }
    }

    fn sum(&self) -> f64 {
        self.partials.iter().fold(0., |p, q| p + *q)
    }
}

impl std::iter::FromIterator<f64> for Sum {
    fn from_iter<T>(iter: T) -> Sum
        where T: IntoIterator<Item=f64>
    {
        let mut e = Sum::new();
        for i in iter {
            e.add(i);
        }
        e
    }
}

fn main() {
    const N: usize = 10_000;
    let mut v = Vec::with_capacity(4 * N);
    for _ in 0..N {
        v.push(1.);
        v.push(1e100);
        v.push(1.);
        v.push(-1e100);
    }
    let stable_sum: Sum = v.iter().map(|x| *x).collect();
    assert_eq!(stable_sum.sum(), (2 * N) as f64);
}

I'm not sure how the interface should ideally look like.

Copied from original issue: rust-num/num#309

[cuviper]

I see, this is from src/libtest/stats.rs.

I'm not sure about using an intermediate object like that. Can it just be an associated function?

trait Float {
...
   fn stable_sum<I: IntoIterator<Item = Self>(iter: I) -> Self { ... }
}

[cuviper]

From @vks on June 29, 2017 7:52

It could be an associated function, but then it would be less general:
incremental calculations would not be possible.

On Thu, Jun 29, 2017, 03:30 Josh Stone notifications@github.com wrote:

I see, this is from src/libtest/stats.rs
https://github.com/rust-lang/rust/blob/69c65d29615c391c958ebf75dd65258ec23e175c/src/libtest/stats.rs#L164-L196
.

I'm not sure about using an intermediate object like that. Can it just be
an associated function?

trait Float {
...
fn stable_sum<I: IntoIterator<Item = Self>(iter: I) -> Self { ... }
}

—
You are receiving this because you authored the thread.
Reply to this email directly, view it on GitHub
rust-num/num#309 (comment), or mute
the thread
https://github.com/notifications/unsubscribe-auth/AACCtIKZpUXsAvKzyGu7ZsffLK9fKU0kks5sIv4xgaJpZM4OGjrD
.

[true-dragon]

I believe the proposed implementation will give incorrectly-rounded results on some inputs. The function

    fn sum(&self) -> f64 {
        self.partials.iter().fold(0., |p, q| p + *q)
    }

sums the partials naively. Why this is a problem is explained here: https://code.activestate.com/recipes/393090/#c16

I've found a testcase that shows the difference: the vector [1e-16, 1.0, 1e16] should sum to 1.0000000000000002e+16, but instead sums to 1.0000000000000000e+16. This error is masked somehow by the 10,000x looping. Is there a reason the current testcase does this?

I'm not familiar enough with Rust to suggest a fix with confidence, but Python's fsum gives the correct result. Its source can be found here: https://hg.python.org/cpython/file/tip/Modules/mathmodule.c#l1112, with the interesting bit starting on line 1199 (The correct line number may change over time; I can't tell how to make a working permanent link.)

[true-dragon]

A note about the N times looping: for the [1e-16, 1, 1e16] case, N <= 2048 correctly causes the assertion to fail, but N > 2048 doesn't. It's a little late in the evening for me to pin down why things change at that power of two. I suspect a rounding error somewhere.