MathUtils 1.0.1

//! Precise floating point summation according to ECMAScript `Math.sumPrecise`.

#[inline]
fn two_sum(a: f64, b: f64) -> (f64, f64) {
    let s = a + b;
    let v = s - a;
    let e = (a - (s - v)) + (b - v);
    (s, e)
}

fn grow_expansion(expansion: &mut Vec<f64>, b: f64) {
    let mut q = b;
    for x in expansion.iter_mut() {
        let (sum, err) = two_sum(q, *x);
        *x = err;
        q = sum;
    }
    if q != 0.0 {
        expansion.push(q);
    }
}

/// Computes the exact sum of an iterator or slice of floating-point numbers without
/// intermediate rounding error, corresponding to ECMAScript `Math.sumPrecise(items)`.
///
/// # Special Values
///
/// - If any value is `NaN` or both `+Infinity` and `-Infinity` are present, returns `NaN`.
/// - If only `+Infinity` is present, returns `+Infinity`.
/// - If only `-Infinity` is present, returns `-Infinity`.
/// - If the collection is empty, returns `0.0`.
/// - If all values in the collection are `-0.0`, returns `-0.0`.
///
/// # Examples
///
/// ```
/// use MathUtils::Math;
///
/// // Demonstrating precision where standard float addition loses digits:
/// let numbers = [1e16, 1.0, -1e16];
/// assert_eq!(Math::sumPrecise(&numbers), 1.0);
/// assert_eq!(Math::sum_precise(&numbers), 1.0);
/// ```
pub fn sum_precise<'a, I>(items: I) -> f64
where
    I: IntoIterator<Item = &'a f64>,
{
    let mut has_nan = false;
    let mut has_pos_inf = false;
    let mut has_neg_inf = false;
    let mut all_neg_zero = true;
    let mut count = 0usize;

    let mut expansion: Vec<f64> = Vec::new();

    for &x in items {
        count += 1;
        if x.is_nan() {
            has_nan = true;
        } else if x.is_infinite() {
            if x.is_sign_positive() {
                has_pos_inf = true;
            } else {
                has_neg_inf = true;
            }
        } else {
            if !(x == 0.0 && x.is_sign_negative()) {
                all_neg_zero = false;
            }
            if x != 0.0 {
                grow_expansion(&mut expansion, x);
            }
        }
    }

    if has_nan || (has_pos_inf && has_neg_inf) {
        return f64::NAN;
    }
    if has_pos_inf {
        return f64::INFINITY;
    }
    if has_neg_inf {
        return f64::NEG_INFINITY;
    }
    if count == 0 {
        return 0.0;
    }
    if all_neg_zero {
        return -0.0;
    }
    if expansion.is_empty() {
        return 0.0;
    }

    // Filter out zeros in expansion and sum
    expansion.retain(|&x| x != 0.0);
    if expansion.is_empty() {
        return 0.0;
    }

    expansion.sort_by(|a, b| a.abs().partial_cmp(&b.abs()).unwrap_or(std::cmp::Ordering::Equal));

    let mut sum = 0.0;
    let mut c = 0.0;
    for x in expansion {
        let y = x - c;
        let t = sum + y;
        c = (t - sum) - y;
        sum = t;
    }
    sum
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn test_sum_precise_basic() {
        assert_eq!(sum_precise(&[]), 0.0);
        assert_eq!(sum_precise(&[1.0, 2.0, 3.0]), 6.0);
        assert_eq!(sum_precise(&[1e16, 1.0, -1e16]), 1.0);
        assert_eq!(sum_precise(&[1e-16, 1.0, -1.0]), 1e-16);
    }

    #[test]
    fn test_sum_precise_special_cases() {
        assert!(sum_precise(&[f64::NAN]).is_nan());
        assert!(sum_precise(&[1.0, f64::NAN]).is_nan());
        assert!(sum_precise(&[f64::INFINITY, f64::NEG_INFINITY]).is_nan());
        assert_eq!(sum_precise(&[f64::INFINITY, 100.0]), f64::INFINITY);
        assert_eq!(sum_precise(&[f64::NEG_INFINITY, -100.0]), f64::NEG_INFINITY);

        let neg_zeros = [-0.0, -0.0];
        let res = sum_precise(&neg_zeros);
        assert_eq!(res, 0.0);
        assert!(res.is_sign_negative());

        let mixed_zeros = [-0.0, 0.0];
        let res2 = sum_precise(&mixed_zeros);
        assert_eq!(res2, 0.0);
        assert!(res2.is_sign_positive());
    }
}