I mean we can add an sum function in std that minimizes the loss of precision. Just like math.fsum in python. In fact, there are already some implementations in library/test/stats.rs(line 141). If you do some more improvements, you can directly add it to f64/f32. Below is my implement, it may not be perfect, which adds INF, NAN and overflow handling based on stats.rs, and made some optimizations at the end with reference to the implementation in Cpython(line 1424). I want to know if it has the value and feasibility of adding to std.

```
fn main() {
let v = vec![1e15, 0.05, 0.05];
assert_eq!(f64::stable_sum(v), 1000000000000000.1);
let v = vec![1e-16, 1., 1e16];
assert_eq!(f64::stable_sum(v), 10000000000000002.0);
// infinity
let v = vec![0., f64::INFINITY];
assert_eq!(f64::stable_sum(v), f64::INFINITY);
// NaN
let v = vec![f64::NEG_INFINITY, f64::INFINITY];
assert!(f64::stable_sum(v).is_nan());
// overflow
let v = vec![f64::MAX, 1e292];
assert_eq!(f64::stable_sum(v), f64::INFINITY);
}
trait StableSum: Sized {
fn stable_sum<I: IntoIterator<Item = f64>>(into_iter: I) -> f64;
}
impl StableSum for f64 {
fn stable_sum<I: IntoIterator<Item = f64>>(into_iter: I) -> f64 {
let mut iter = into_iter.into_iter();
let mut partials = vec![];
let mut special_sum = 0.;
let mut inf_sum = 0.;
while let Some(mut x) = iter.next() {
let xsave = x;
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..partials.len() {
let mut y: f64 = partials[i];
if x.abs() < y.abs() {
let t = x;
x = y;
y = t;
}
// 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 {
partials[j] = lo;
j += 1;
}
x = hi;
}
if (!x.is_finite()) {
/* a nonfinite x could arise either as
a result of intermediate overflow, or
as a result of a nan or inf in the
summands */
if (xsave.is_finite()) {
//intermediate overflow
return x;
}
if (xsave.is_infinite()) {
inf_sum += xsave;
}
special_sum += xsave;
// reset partials
partials.clear();
}
if j >= partials.len() {
partials.push(x);
} else {
partials[j] = x;
partials.truncate(j + 1);
}
}
if (special_sum != 0.0) {
if (inf_sum.is_nan()) {
return f64::NAN;
} else {
return special_sum;
}
}
let mut len = partials.len();
let mut hi = partials[len - 1];
let mut lo = 0.;
len -= 1;
while len > 0 {
let x = hi;
let y = partials[len - 1];
hi = x + y;
lo = y - (hi - x);
if lo != 0. {
break;
}
}
if len > 1 && ((lo < 0. && partials[len - 2] < 0.) || (lo > 0. && partials[len - 2] > 0.)) {
let y = lo + lo;
let x = hi + y;
if y == x - hi {
hi = x;
}
}
hi
}
}
```