msum.rs

extern mod std;
extern mod extra;

use std::{num, util};


fn msum_10927(vals: &[f64]) -> f64 {
    let mut partials : ~[f64] = ~[];
    // `partials` is already borrowed when looping
    // use `modifications` to update `partials`
    let mut modifications: ~[f64] = ~[];

    for &mut x in vals.iter() {
        let mut i = 0;
        modifications.truncate(0);
        // This inner loop applies hi/lo summation to each
        // partial so that the list of partial sums remains exact.
        for &mut y in partials.iter() {
            if num::abs(x) < num::abs(y) {
                util::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 {
                modifications.push(lo);
                i += 1;
            }
            x = hi;
        }
        for (f1, partial) in modifications.iter().
            zip(partials.mut_iter()) {
            *partial = *f1;
        }
        if partials.len() > i as uint {
            partials[i] = x;
        } else {
            partials.push(x);
        }
    }
    partials.iter().fold(0.0, |p, q| p + *q)
}

fn msum_10927_swap(vals: &[f64]) -> f64 {
    let mut partials : ~[f64] = ~[];
    // `partials` is already borrowed when looping
    // use `modifications` to update `partials`
    let mut modifications: ~[f64] = ~[];

    for &mut x in vals.iter() {
        modifications.truncate(0);
        // This inner loop applies hi/lo summation to each
        // partial so that the list of partial sums remains exact.
        for &mut y in partials.iter() {
            if num::abs(x) < num::abs(y) {
                util::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 {
                modifications.push(lo);
            }
            x = hi;
        }
        if x != 0.0 {
            modifications.push(x);
        }
        util::swap(&mut modifications, &mut partials);
    }
    partials.iter().fold(0.0, |p, q| p + *q)
}

fn msum_mine(vals: &[f64]) -> f64 {
    let mut partials: ~[f64] = ~[];
    for &mut x in vals.iter() {
        let mut i = 0;
        let mut j = 0;
        while i < partials.len() {
            let mut y = partials[i];
            if num::abs(x) < num::abs(y) {
                util::swap(&mut x, &mut y);
            }
            let hi = x + y;
            let lo = y - (hi - x);
            if lo != 0f64 {
                partials[j] = lo;
                j += 1;
            }
            x = hi;
            i += 1;
        }
        if j >= partials.len() {
            partials.push(x);
        } else {
            partials[j] = x;
            partials.truncate(j+1);
        }
    }
    //TODO: can't just sum partials, need to fix
    partials.iter().fold(0.0, |a, &b| a+b)
}


#[cfg(test)]
mod tests {

    use super::{msum_10927, msum_10927_swap, msum_mine};

    // test cases from http://bugs.python.org/file10357/msum4.py
    fn cases() -> ~[(~[f64], f64)] {
        ~[
            (~[1e30, 1.2, -1e30], 1.2),
            (~[0.5, 3.2321, 1.5678], 5.2999),
            (~[1.0, 1.0e100, 1.0, -1.0e100], 2.0),
            (~[], 0.0),
            (~[0.0], 0.0),
            (~[1e100, 1.0, -1e100, 1e-100, 1e50, -1.0, -1e50], 1e-100),
            // (~[1e308, 1e308, -1e308], 1e308),
            // (~[-1e308, 1e308, 1e308], 1e308),
            // (~[1e308, -1e308, 1e308], 1e308),
            // (~[2.0**1023, 2.0**1023, -2.0**1000], 1.7976930277114552e+308),
            // (~[twopow, twopow, twopow, twopow, -twopow, -twopow, -twopow], 8.9884656743115795e+307),
            // (~[2.0**53, -0.5, -2.0**-54], 2.0**53-1.0),
            // (~[2.0**53, 1.0, 2.0**-100], 2.0**53+2.0),
            // (~[2.0**53+10.0, 1.0, 2.0**-100], 2.0**53+12.0),
            // (~[2.0**53-4.0, 0.5, 2.0**-54], 2.0**53-3.0),
            // (~[2.0**1023-2.0**970, -1.0, 2.0**1023], 1.7976931348623157e+308),
            // (~[float_info.max, float_info.max*2.**-54], float_info.max),
            // (~[float_info.max, float_info.max*2.**-53], OverflowError),
            // (~[1./n for n in range(1, 1001)], 7.4854708605503451),
            // (~[(-1.)**n/n for n in range(1, 1001)], -0.69264743055982025),
            // (~[1.7**(i+1)-1.7**i for i in range(1000)] + [-1.7**1000], -1.0),
            // (~[inf, -inf, nan], nan),
            // (~[nan, inf, -inf], nan),
            // (~[inf, nan, inf], nan),
            // (~[inf, inf], inf),
            // (~[inf, -inf], ValueError),
            // (~[-inf, 1e308, 1e308, -inf], -inf),
            // (~[2.0**1023-2.0**970, 0.0, 2.0**1023], OverflowError),
            // (~[2.0**1023-2.0**970, 1.0, 2.0**1023], OverflowError),
            // (~[2.0**1023, 2.0**1023], OverflowError),
            // (~[2.0**1023, 2.0**1023, -1.0], OverflowError),
            // (~[twopow, twopow, twopow, twopow, -twopow, -twopow], OverflowError),
            // (~[twopow, twopow, twopow, twopow, -twopow, twopow], OverflowError),
            // (~[-twopow, -twopow, -twopow, -twopow], OverflowError),
            // (~[2.**1023, 2.**1023, -2.**971], float_info.max),
            // (~[2.**1023, 2.**1023, -2.**970], OverflowError),
            // (~[-2.**970, 2.**1023, 2.**1023, -2.**-1074], float_info.max),
            // (~[2.**1023, 2.**1023, -2.**970, 2.**-1074], OverflowError),
            // (~[-2.**1023, 2.**971, -2.**1023], -float_info.max),
            // (~[-2.**1023, -2.**1023, 2.**970], OverflowError),
            // (~[-2.**1023, -2.**1023, 2.**970, 2.**-1074], -float_info.max),
            // (~[-2.**-1074, -2.**1023, -2.**1023, 2.**970], OverflowError),
            // (~[2.**930, -2.**980, 2.**1023, 2.**1023, twopow, -twopow], 1.7976931348622137e+308),
            // (~[2.**1023, 2.**1023, -1e307], 1.6976931348623159e+308),
            // (~[1e16, 1., 1e-16], 10000000000000002.0)
        ]
    }

    #[test]
    fn test_msum_10927() {
        let cs = cases();
        for (arr, s) in cs.move_iter() {
            assert_eq!(msum_10927(arr), s);
        }
    }

    #[test]
    fn test_msum_10927_swap() {
        let cs = cases();
        for (arr, s) in cs.move_iter() {
            assert_eq!(msum_10927_swap(arr), s);
        }
    }

    #[test]
    fn test_msum_mine() {
        let cs = cases();
        for (arr, s) in cs.move_iter() {
            assert_eq!(msum_mine(arr), s);
        }
    }
}

#[cfg(test)]
mod bench {
    use extra::test::BenchHarness;
    use super::{msum_10927, msum_10927_swap, msum_mine};


    #[bench]
    fn msum_10927_single_f64(bh: &mut BenchHarness) {
        bh.iter(|| {
            msum_10927([0.0]);
        })
    }
    #[bench]
    fn msum_10927_three_items(bh: &mut BenchHarness) {
        bh.iter(|| {
            msum_10927([1e20, 1.5, -1e20]);
        })
    }

    #[bench]
    fn msum_10927_ten_items(bh: &mut BenchHarness) {
        bh.iter(|| {
            msum_10927([1e20, 1.5, -1e20, 1e-20, -1e-20, 0.1, 1.0, -1.0, -0.1, 10.0]);
        })
    }

    #[bench]
    fn msum_10927_swap_single_f64(bh: &mut BenchHarness) {
        bh.iter(|| {
            msum_10927_swap([0.0]);
        })
    }
    #[bench]
    fn msum_10927_swap_three_items(bh: &mut BenchHarness) {
        bh.iter(|| {
            msum_10927_swap([1e20, 1.5, -1e20]);
        })
    }

    #[bench]
    fn msum_10927_swap_ten_items(bh: &mut BenchHarness) {
        bh.iter(|| {
            msum_10927_swap([1e20, 1.5, -1e20, 1e-20, -1e-20, 0.1, 1.0, -1.0, -0.1, 10.0]);
        })
    }

    #[bench]
    fn msum_mine_single_f64(bh: &mut BenchHarness) {
        bh.iter(|| {
            msum_mine([0.0]);
        })
    }
    #[bench]
    fn msum_mine_three_items(bh: &mut BenchHarness) {
        bh.iter(|| {
            msum_mine([1e20, 1.5, -1e20]);
        })
    }
    #[bench]
    fn msum_mine_ten_items(bh: &mut BenchHarness) {
        bh.iter(|| {
            msum_mine([1e20, 1.5, -1e20, 1e-20, -1e-20, 0.1, 1.0, -1.0, -0.1, 10.0]);
        })
    }

}