Last active
December 31, 2015 07:39
-
-
Save jfager/7955901 to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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]); | |
}) | |
} | |
} |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment