dashu_float/
iter.rs

1//! Implementation of core::iter traits
2
3use crate::{
4    error::assert_finite,
5    fbig::FBig,
6    repr::{Context, Repr, Word},
7    round::Round,
8};
9use core::iter::{Product, Sum};
10
11/// Correctly-rounded summation of finite floats.
12///
13/// Every addend is accumulated exactly at the [`Repr`] level (no per-step rounding); the exact total
14/// is then rounded a single time to the target context. The target context is `Context::max` over all
15/// addend contexts (matching chained `+`), and the final round reuses [`Context::repr_round`] — this
16/// yields the same result as rounding the mathematically exact sum (MPFR `mpfr_sum` semantics).
17///
18/// Because the accumulator is exact, summing addends with widely differing exponents can grow the
19/// intermediate significand to span the full exponent range of the inputs.
20fn precise_sum<R: Round, const B: Word>(
21    mut iter: impl Iterator<Item = (Repr<B>, Context<R>)>,
22) -> FBig<R, B> {
23    let (mut acc, mut context, mut all_neg_zero) = match iter.next() {
24        Some((repr, ctx)) => {
25            assert_finite(&repr);
26            let all_neg_zero = repr.is_neg_zero();
27            (repr, ctx, all_neg_zero)
28        }
29        None => return FBig::ZERO, // empty iterator → additive identity
30    };
31    for (repr, ctx) in iter {
32        assert_finite(&repr);
33        all_neg_zero &= repr.is_neg_zero();
34        acc = acc + &repr;
35        context = Context::max(context, ctx);
36    }
37    // Exact cancellation can leave `acc` with a zero significand whose exponent coincides with the
38    // `-0` sentinel (-1) — `Repr::new` then mislabels it `-0`. Canonicalize the sign per IEEE 754
39    // §6.3 folded left, as the chained `+` does: the exact zero sum is `-0` under
40    // roundTowardNegative or when *every* addend is `-0` (a sum of like-signed zeros keeps the
41    // sign), and `+0` otherwise (cancellation of nonzero addends, mixed-sign zeros).
42    if acc.significand.is_zero() {
43        acc = if R::IS_ROUND_TOWARD_NEGATIVE || all_neg_zero {
44            Repr::neg_zero()
45        } else {
46            Repr::zero()
47        };
48    }
49    FBig::new(context.repr_round(acc).value(), context)
50}
51
52impl<R: Round, const B: Word> Sum for FBig<R, B> {
53    fn sum<I: Iterator<Item = Self>>(iter: I) -> Self {
54        precise_sum(iter.map(|v| (v.repr, v.context)))
55    }
56}
57
58impl<'a, R: Round, const B: Word> Sum<&'a FBig<R, B>> for FBig<R, B> {
59    fn sum<I: Iterator<Item = &'a FBig<R, B>>>(iter: I) -> Self {
60        precise_sum(iter.map(|v| (v.repr.clone(), v.context)))
61    }
62}
63
64impl<R: Round, const B: Word> Product for FBig<R, B> {
65    fn product<I: Iterator<Item = Self>>(iter: I) -> Self {
66        iter.fold(FBig::ONE, |acc, x| acc * x)
67    }
68}
69
70impl<'a, R: Round, const B: Word> Product<&'a FBig<R, B>> for FBig<R, B> {
71    fn product<I: Iterator<Item = &'a FBig<R, B>>>(iter: I) -> Self {
72        iter.fold(FBig::ONE, |acc, x| acc * x)
73    }
74}
75
76#[cfg(test)]
77mod tests {
78    use super::*;
79    use crate::round::mode::{HalfAway, HalfEven, Zero};
80    use alloc::vec::Vec;
81    use core::str::FromStr;
82    use dashu_int::IBig;
83
84    type F = FBig<HalfAway, 10>;
85
86    fn r<const B: Word>(sig: i128, exp: isize) -> Repr<B> {
87        Repr::new(IBig::from(sig), exp)
88    }
89
90    #[test]
91    fn sum_empty() {
92        let s: F = core::iter::empty::<F>().sum();
93        assert_eq!(s, F::ZERO);
94        assert!(s.repr().is_pos_zero());
95    }
96
97    #[test]
98    fn sum_single() {
99        let a = F::from_str("1.234").unwrap();
100        let s: F = core::iter::once(a.clone()).sum();
101        assert_eq!(s, a);
102    }
103
104    #[test]
105    fn sum_cancellation_is_positive_zero() {
106        // 0.5 + (-0.5) must yield +0 under round-to-nearest (IEEE 754 §6.3), matching chained `+` —
107        // the exact-cancellation exponent (-1) must not be mistaken for the `-0` sentinel.
108        let half = F::from_str("0.5").unwrap();
109        let s: F = [half.clone(), -half].into_iter().sum();
110        assert!(s.repr().is_pos_zero());
111        assert!(!s.repr().is_neg_zero());
112    }
113
114    // The sign of a zero sum follows IEEE 754 §6.3 folded left: like-signed zeros keep their
115    // sign, a single addend is itself, and mixed-sign zeros cancel to +0 (nearest).
116    #[test]
117    fn sum_keeps_sign_of_all_neg_zero_addends() {
118        let nz = -F::ZERO; // -0
119        assert!(nz.repr().is_neg_zero());
120        let s: F = core::iter::once(nz.clone()).sum();
121        assert!(s.repr().is_neg_zero(), "sum of one -0 is -0");
122        let s: F = [nz.clone(), nz.clone()].into_iter().sum();
123        assert!(s.repr().is_neg_zero(), "sum of two -0 is -0");
124        let s: F = [nz, F::ZERO].into_iter().sum();
125        assert!(s.repr().is_pos_zero(), "-0 + (+0) is +0 (nearest)");
126    }
127
128    #[test]
129    fn sum_zero_is_neg_zero_under_down() {
130        use crate::round::mode::Down;
131        type FD = FBig<Down, 10>;
132        // exact cancellation of nonzero addends → -0 under roundTowardNegative
133        let one = FD::from_str("1.0").unwrap();
134        let s: FD = [one.clone(), -one].into_iter().sum();
135        assert!(s.repr().is_neg_zero());
136    }
137
138    #[test]
139    fn sum_ref_iter() {
140        let a = F::from_str("1.23").unwrap();
141        let b = F::from_str("4.56").unwrap();
142        let vals = [a, b];
143        let s: F = vals.iter().sum();
144        assert_eq!(s, F::from_str("5.79").unwrap());
145    }
146
147    #[test]
148    fn sum_exact_cancellation() {
149        let a = F::from_str("1.00").unwrap();
150        let b = F::from_str("-1.00").unwrap();
151        let s: F = [a, b].into_iter().sum();
152        assert!(s.repr().is_pos_zero()); // exact zero → +0
153    }
154
155    #[test]
156    fn sum_uses_max_precision() {
157        let a = F::from_str("1.234").unwrap(); // precision 4
158        let b = F::from_str("5.6").unwrap(); // precision 2
159        let s: F = [a, b].into_iter().sum();
160        assert_eq!(s, F::from_str("6.834").unwrap());
161        assert_eq!(s.precision(), 4); // max of the two operand precisions
162    }
163
164    // The headline regression: per-step truncation loses almost all the mass, the precise sum does
165    // not. 11 copies of 0.9 at precision 1, truncation: exact sum 9.9 → 9, naive fold → 1.
166    #[test]
167    fn sum_precise_beats_naive_truncation() {
168        type Z = FBig<Zero, 10>;
169        let vals: Vec<Z> = (0..11)
170            .map(|_| Z::from_parts(9.into(), -1)) // 0.9, precision 1
171            .collect();
172
173        let precise: Z = vals.iter().cloned().sum();
174        assert_eq!(precise, Z::from_parts(9.into(), 0)); // 9.9 truncated to 1 digit = 9
175
176        let naive: Z = vals.iter().cloned().fold(Z::ZERO, |acc, v| acc + v);
177        assert_eq!(naive, Z::from_parts(1.into(), 0)); // 0.9+0.9=1.8→1, then stuck at 1
178        assert_ne!(precise, naive);
179    }
180
181    // Cross-check against an independent oracle: an exact fold at unlimited precision (where `+` is
182    // exact) rounded once must equal the precise sum.
183    #[test]
184    fn sum_matches_unlimited_oracle() {
185        fn check<const B: Word>(strs: &[&str]) {
186            let vals: Vec<FBig<HalfEven, B>> = strs
187                .iter()
188                .map(|s| FBig::<HalfEven, B>::from_str(s).unwrap())
189                .collect();
190            let target = vals.iter().map(|v| v.precision()).max().unwrap_or(0);
191
192            let precise: FBig<HalfEven, B> = vals.clone().into_iter().sum();
193
194            let exact: FBig<HalfEven, B> = vals
195                .iter()
196                .map(|v| v.clone().with_precision(0).value())
197                .fold(FBig::<HalfEven, B>::ZERO, |acc, v| acc + v);
198            let oracle = exact.with_precision(target).value();
199
200            assert_eq!(precise.precision(), oracle.precision());
201            assert_eq!(precise, oracle);
202        }
203        check::<10>(&["1.234", "5.6", "0.001"]);
204        check::<10>(&["9.99", "0.009", "0.0009"]);
205        check::<2>(&["1.1", "0.01", "0.001", "0.0001"]); // 1.5, 0.25, 0.125, 0.0625
206    }
207
208    #[test]
209    fn product_owned_and_ref() {
210        let a = F::from_str("1.2").unwrap();
211        let b = F::from_str("3.0").unwrap();
212        let expected = &a * &b;
213        let owned: F = [a.clone(), b.clone()].into_iter().product();
214        let by_ref: F = [a, b].iter().product();
215        assert_eq!(owned, expected);
216        assert_eq!(by_ref, expected);
217    }
218
219    #[test]
220    #[should_panic(expected = "arithmetic operations with the infinity are not allowed")]
221    fn sum_infinite_panics() {
222        let _: F = core::iter::once(F::INFINITY).sum();
223    }
224
225    // Exercises the zero short-circuit in `impl Add for &Repr`: a `-0` sentinel exponent must not
226    // corrupt the accumulator exponent.
227    #[test]
228    fn repr_add_neg_zero_is_identity() {
229        let nz = Repr::<10>::neg_zero();
230        let x = r::<10>(5, -2); // 0.05
231                                // -0 + x == x, x + -0 == x (exact, pre-rounding)
232        assert_eq!(&nz + &x, x);
233        assert_eq!(&x + &nz, x);
234        assert_eq!(&nz + &nz, Repr::<10>::neg_zero());
235        // 1 + -1 cancels to +0 (not -0)
236        assert_eq!(&r::<10>(1, 0) + &r::<10>(-1, 0), Repr::<10>::zero());
237    }
238}
dashu_float/iter.rs

dashu_float/
iter.rs
