Skip to main content

dashu_float/
add.rs

1use crate::{
2    error::assert_finite_operands,
3    fbig::FBig,
4    helper_macros,
5    repr::{Context, Repr, Word},
6    round::{Round, Rounded},
7    utils::{digit_len, shl_digits, shl_digits_in_place, split_digits, split_digits_ref},
8};
9use core::{
10    cmp::Ordering,
11    ops::{Add, AddAssign, Sub, SubAssign},
12};
13
14use dashu_base::Sign::{self, *};
15use dashu_int::{IBig, UBig};
16
17impl<R: Round, const B: Word> Add for FBig<R, B> {
18    type Output = Self;
19
20    #[inline]
21    fn add(self, rhs: Self) -> Self::Output {
22        add_val_val(self, rhs, Positive)
23    }
24}
25
26impl<R: Round, const B: Word> Add<&FBig<R, B>> for FBig<R, B> {
27    type Output = Self;
28
29    #[inline]
30    fn add(self, rhs: &FBig<R, B>) -> Self::Output {
31        add_val_ref(self, rhs, Positive)
32    }
33}
34
35impl<R: Round, const B: Word> Add<FBig<R, B>> for &FBig<R, B> {
36    type Output = FBig<R, B>;
37
38    #[inline]
39    fn add(self, rhs: FBig<R, B>) -> Self::Output {
40        add_ref_val(self, rhs, Positive)
41    }
42}
43
44impl<R: Round, const B: Word> Add<&FBig<R, B>> for &FBig<R, B> {
45    type Output = FBig<R, B>;
46
47    #[inline]
48    fn add(self, rhs: &FBig<R, B>) -> Self::Output {
49        add_ref_ref(self, rhs, Positive)
50    }
51}
52
53impl<R: Round, const B: Word> Sub for FBig<R, B> {
54    type Output = Self;
55
56    #[inline]
57    fn sub(self, rhs: Self) -> Self::Output {
58        add_val_val(self, rhs, Negative)
59    }
60}
61
62impl<R: Round, const B: Word> Sub<&FBig<R, B>> for FBig<R, B> {
63    type Output = Self;
64
65    #[inline]
66    fn sub(self, rhs: &FBig<R, B>) -> Self::Output {
67        add_val_ref(self, rhs, Negative)
68    }
69}
70
71impl<R: Round, const B: Word> Sub<FBig<R, B>> for &FBig<R, B> {
72    type Output = FBig<R, B>;
73
74    #[inline]
75    fn sub(self, rhs: FBig<R, B>) -> Self::Output {
76        add_ref_val(self, rhs, Negative)
77    }
78}
79
80impl<R: Round, const B: Word> Sub<&FBig<R, B>> for &FBig<R, B> {
81    type Output = FBig<R, B>;
82
83    #[inline]
84    fn sub(self, rhs: &FBig<R, B>) -> Self::Output {
85        add_ref_ref(self, rhs, Negative)
86    }
87}
88
89helper_macros::impl_binop_assign_by_taking!(impl AddAssign<Self>, add_assign, add);
90helper_macros::impl_binop_assign_by_taking!(impl SubAssign<Self>, sub_assign, sub);
91
92macro_rules! impl_add_sub_primitive_with_fbig {
93    ($($t:ty)*) => {$(
94        helper_macros::impl_binop_with_primitive!(impl Add<$t>, add);
95        helper_macros::impl_binop_assign_with_primitive!(impl AddAssign<$t>, add_assign);
96        helper_macros::impl_binop_with_primitive!(impl Sub<$t>, sub);
97        helper_macros::impl_binop_assign_with_primitive!(impl SubAssign<$t>, sub_assign);
98    )*};
99}
100impl_add_sub_primitive_with_fbig!(u8 u16 u32 u64 u128 usize UBig i8 i16 i32 i64 i128 isize IBig);
101
102fn add_val_val<R: Round, const B: Word>(
103    lhs: FBig<R, B>,
104    mut rhs: FBig<R, B>,
105    rhs_sign: Sign,
106) -> FBig<R, B> {
107    assert_finite_operands(&lhs.repr, &rhs.repr);
108
109    let context = Context::max(lhs.context, rhs.context);
110    rhs.repr.significand *= rhs_sign;
111    let sum = if lhs.repr.is_zero() {
112        rhs.repr
113    } else if rhs.repr.is_zero() {
114        lhs.repr
115    } else {
116        match lhs.repr.exponent.cmp(&rhs.repr.exponent) {
117            Ordering::Equal => context.repr_round(Repr::new(
118                lhs.repr.significand + rhs.repr.significand,
119                lhs.repr.exponent,
120            )),
121            Ordering::Greater => context.repr_add_large_small(lhs.repr, &rhs.repr, Positive),
122            Ordering::Less => context.repr_add_small_large(lhs.repr, &rhs.repr, Positive),
123        }
124        .value()
125    };
126    FBig::new(sum, context)
127}
128
129fn add_val_ref<R: Round, const B: Word>(
130    lhs: FBig<R, B>,
131    rhs: &FBig<R, B>,
132    rhs_sign: Sign,
133) -> FBig<R, B> {
134    assert_finite_operands(&lhs.repr, &rhs.repr);
135
136    let context = Context::max(lhs.context, rhs.context);
137    let sum = if lhs.repr.is_zero() {
138        let mut repr = rhs.repr.clone();
139        repr.significand *= rhs_sign;
140        repr
141    } else if rhs.repr.is_zero() {
142        lhs.repr
143    } else {
144        match lhs.repr.exponent.cmp(&rhs.repr.exponent) {
145            Ordering::Equal => {
146                let sum_signif = match rhs_sign {
147                    Positive => lhs.repr.significand + &rhs.repr.significand,
148                    Negative => lhs.repr.significand - &rhs.repr.significand,
149                };
150                context.repr_round(Repr::new(sum_signif, lhs.repr.exponent))
151            }
152            Ordering::Greater => context.repr_add_large_small(lhs.repr, &rhs.repr, rhs_sign),
153            Ordering::Less => context.repr_add_small_large(lhs.repr, &rhs.repr, rhs_sign),
154        }
155        .value()
156    };
157    FBig::new(sum, context)
158}
159
160fn add_ref_val<R: Round, const B: Word>(
161    lhs: &FBig<R, B>,
162    mut rhs: FBig<R, B>,
163    rhs_sign: Sign,
164) -> FBig<R, B> {
165    assert_finite_operands(&lhs.repr, &rhs.repr);
166
167    let context = Context::max(lhs.context, rhs.context);
168    rhs.repr.significand *= rhs_sign;
169    let sum = if lhs.repr.is_zero() {
170        rhs.repr
171    } else if rhs.repr.is_zero() {
172        lhs.repr.clone()
173    } else {
174        match lhs.repr.exponent.cmp(&rhs.repr.exponent) {
175            Ordering::Equal => context.repr_round(Repr::new(
176                &lhs.repr.significand + rhs.repr.significand,
177                lhs.repr.exponent,
178            )),
179            Ordering::Greater => context.repr_add_small_large(rhs.repr, &lhs.repr, Positive),
180            Ordering::Less => context.repr_add_large_small(rhs.repr, &lhs.repr, Positive),
181        }
182        .value()
183    };
184    FBig::new(sum, context)
185}
186
187fn add_ref_ref<R: Round, const B: Word>(
188    lhs: &FBig<R, B>,
189    rhs: &FBig<R, B>,
190    rhs_sign: Sign,
191) -> FBig<R, B> {
192    assert_finite_operands(&lhs.repr, &rhs.repr);
193
194    let context = Context::max(lhs.context, rhs.context);
195    let sum = if lhs.repr.is_zero() {
196        let mut repr = rhs.repr.clone();
197        repr.significand *= rhs_sign;
198        repr
199    } else if rhs.repr.is_zero() {
200        lhs.repr.clone()
201    } else {
202        match lhs.repr.exponent.cmp(&rhs.repr.exponent) {
203            Ordering::Equal => context.repr_round(Repr::new(
204                &lhs.repr.significand + rhs_sign * rhs.repr.significand.clone(),
205                lhs.repr.exponent,
206            )),
207            Ordering::Greater => {
208                context.repr_add_large_small(lhs.repr.clone(), &rhs.repr, rhs_sign)
209            }
210            Ordering::Less => context.repr_add_small_large(lhs.repr.clone(), &rhs.repr, rhs_sign),
211        }
212        .value()
213    };
214    FBig::new(sum, context)
215}
216
217impl<R: Round> Context<R> {
218    /// Round sum = `significand * B ^ exponent` with the low part (value, precision).
219    /// If the sum is actually from a subtraction and the low part is not zero, `is_sub` should be true.
220    fn repr_round_sum<const B: Word>(
221        &self,
222        mut significand: IBig,
223        mut exponent: isize,
224        mut low: (IBig, usize),
225        is_sub: bool,
226    ) -> Rounded<Repr<B>> {
227        if !self.is_limited() {
228            // short cut for unlimited precision
229            return Rounded::Exact(Repr::new(significand, exponent));
230        }
231
232        // use one extra digit to prevent cancellation in rounding
233        let rnd_precision = self.precision + is_sub as usize;
234
235        // align to precision again
236        let digits = digit_len::<B>(&significand);
237        match digits.cmp(&rnd_precision) {
238            Ordering::Equal => {}
239            Ordering::Greater => {
240                // Shrink if the result has more digits than desired precision
241                /*
242                 * lhs:         |=========0000|
243                 * rhs:              |========|xxxxx|
244                 * sum:        |==============|xxxxx|
245                 * precision:  |<----->|
246                 * shrink:     |=======|xxxxxxxxxxxx|
247                 */
248                let shift = digits - rnd_precision;
249                let (signif_hi, mut signif_lo) = split_digits::<B>(significand, shift);
250                significand = signif_hi;
251                exponent += shift as isize;
252                shl_digits_in_place::<B>(&mut signif_lo, low.1);
253                low.0 += signif_lo;
254                low.1 += shift;
255            }
256            Ordering::Less => {
257                // Expand to low parts if the result has less digits than desired precision.
258                /*
259                 * A possible case when lhs and rhs have different sign:
260                 * lhs:  |=========0000|
261                 * rhs:  |=============|xxxxx|
262                 * sum:          |=====|xxxxx|
263                 * precision+1:  |<------>|
264                 * shift:              |<>|
265                 * expanded:     |========|xx|
266                 */
267                if !low.0.is_zero() {
268                    let (low_val, low_prec) = low;
269                    let shift = low_prec.min(rnd_precision - digits);
270                    let (pad, low_val) = split_digits::<B>(low_val, low_prec - shift);
271                    shl_digits_in_place::<B>(&mut significand, shift);
272                    exponent -= shift as isize;
273                    significand += pad;
274                    low = (low_val, low_prec - shift);
275                }
276            }
277        };
278
279        // perform rounding
280        if low.0.is_zero() {
281            Rounded::Exact(Repr::new(significand, exponent))
282        } else {
283            // By now significand should have at least full precision. After adjustment, the digits length
284            // could be one more than the precision. We don't shrink the extra digit.
285            let adjust = R::round_fract::<B>(&significand, low.0, low.1);
286            Rounded::Inexact(Repr::new(significand + adjust, exponent), adjust)
287        }
288    }
289
290    // lhs + rhs_sign * rhs, assuming lhs.exponent >= rhs.exponent
291    fn repr_add_large_small<const B: Word>(
292        &self,
293        mut lhs: Repr<B>,
294        rhs: &Repr<B>,
295        rhs_sign: Sign,
296    ) -> Rounded<Repr<B>> {
297        debug_assert!(lhs.exponent >= rhs.exponent);
298
299        // use one extra digit when subtracting to prevent cancellation in rounding
300        let is_sub = lhs.significand.sign() != rhs_sign * rhs.significand.sign();
301        let rnd_precision = self.precision + is_sub as usize;
302
303        let ediff = (lhs.exponent - rhs.exponent) as usize;
304        let ldigits = lhs.digits();
305        let rdigits_est = rhs.digits_ub(); // overestimate
306
307        // align the exponent
308        let low: (IBig, usize); // (value of low part, precision of the low part)
309        let (significand, exponent) =
310            if self.is_limited() && is_sub && rdigits_est + self.precision >= ldigits + ediff {
311                // The smaller operand (`rhs`, lower exponent) reaches the larger
312                // operand's `precision`-digit window — its top digit is at or above the
313                // window edge (`rdigits + precision >= ldigits + ediff`, i.e. `rhs`'s top
314                // position `rdigits - ediff` is `>= ldigits - precision`). An effective
315                // subtraction can then cancel and lose leading digits, which the trimmed
316                // path cannot recover (its single re-expand in `repr_round_sum` collapses
317                // a genuinely small difference to the wrong value — e.g. `1.00 -
318                // 0.99999999` at precision 3 to `0` instead of `1e-8`, or `0.5 - 0.4375`
319                // at precision 1 to `0` instead of `0.0625`). So form the exact difference
320                // at full operand width and let the shared `repr_round_sum` round it once
321                // (with the same guard digit as the trimmed path, so no low tail is
322                // needed). The complement (`<`) is the trimmed/negligible region where
323                // `rhs` stays strictly below the window and no cancellation is possible.
324                shl_digits_in_place::<B>(&mut lhs.significand, ediff);
325                low = (IBig::ZERO, 0);
326                match rhs_sign {
327                    Positive => (lhs.significand + &rhs.significand, rhs.exponent),
328                    Negative => (lhs.significand - &rhs.significand, rhs.exponent),
329                }
330            } else if self.is_limited()
331                && rdigits_est + 1 < ediff
332                && rdigits_est + 1 + rnd_precision < ldigits + ediff
333            {
334                // rhs is entirely below lhs's rounding window, so only its sign
335                // contributes to the rounding; replace it with a unit sticky tail
336                // (`|low| = 1`).
337                //
338                // The sticky must be positioned at rhs's *real* magnitude, i.e. `ediff`
339                // digits below lhs's exponent — NOT at `precision - ldigits`. Positioning
340                // by `ediff` keeps the sticky genuinely sub-ULP (|1| << B^ediff, and the
341                // branch guard guarantees ediff >= 3), so it can never land on a rounding
342                // tie. Positioning by `precision - ldigits` instead let the re-expand drag
343                // the sticky up to the LSB, where for base 2 + round-to-nearest it equals
344                // exactly half (1 == B^0 == ½·B^1) and injected a spurious ULP — e.g.
345                // `1 + 2^-100` at precision 10 returned `513·2^-9` instead of `1`.
346                low = (rhs_sign * rhs.significand.signum(), ediff);
347                (lhs.significand, lhs.exponent)
348            } else if self.is_limited() && ldigits >= self.precision {
349                // if the lhs already exceeds the desired precision, just align rhs
350                /* Before:
351                 * lhs: |==============|
352                 * rhs:      |==============|
353                 *              ediff  |<-->|
354                 *    precision  |<--->|
355                 *
356                 * After:
357                 * lhs: |==============|
358                 * rhs:      |=========|xxxx|
359                 */
360                let (rhs_signif, r) = split_digits_ref::<B>(&rhs.significand, ediff);
361                low = (rhs_sign * r, ediff);
362                (lhs.significand + rhs_sign * rhs_signif, lhs.exponent)
363            } else if self.is_limited() && ediff + ldigits > self.precision {
364                // if the shifted lhs exceeds the desired precision, align lhs and rhs to precision
365                /* Before:
366                 * lhs: |=========|
367                 * rhs:      |==============|
368                 *                |< ediff >|
369                 *      |< precision >|
370                 *
371                 * After:
372                 * lhs: |=========0000|
373                 * rhs:      |========|xxxxx|
374                 *        lshift  |<->|
375                 *            rshift  |<--->|
376                 */
377                let lshift = self.precision - ldigits;
378                let rshift = ediff - lshift;
379                let (rhs_signif, r) = split_digits_ref::<B>(&rhs.significand, rshift);
380                shl_digits_in_place::<B>(&mut lhs.significand, lshift);
381
382                low = (rhs_sign * r, rshift);
383                (lhs.significand + rhs_sign * rhs_signif, lhs.exponent - lshift as isize)
384            } else {
385                // otherwise directly shift lhs to required position
386                /* Before:
387                 * lhs: |==========|
388                 * rhs:       |==============|
389                 *                 |< ediff >|
390                 *      |<------ precision ------>|
391                 *
392                 * After:
393                 * lhs: |==========0000000000|
394                 * rhs:       |==============|
395                 */
396                shl_digits_in_place::<B>(&mut lhs.significand, ediff);
397                low = (IBig::ZERO, 0);
398                match rhs_sign {
399                    Positive => (lhs.significand + &rhs.significand, rhs.exponent),
400                    Negative => (lhs.significand - &rhs.significand, rhs.exponent),
401                }
402            };
403
404        self.repr_round_sum(significand, exponent, low, is_sub)
405    }
406
407    // lhs + rhs_sign * rhs, assuming lhs.exponent <= rhs.exponent
408    fn repr_add_small_large<const B: Word>(
409        &self,
410        lhs: Repr<B>,
411        rhs: &Repr<B>,
412        rhs_sign: Sign,
413    ) -> Rounded<Repr<B>> {
414        debug_assert!(lhs.exponent <= rhs.exponent);
415
416        // the following implementation should be exactly the same as `repr_add_large_small`
417        // other than lhs and rhs are swapped. See `repr_add_large_small` for full documentation
418        let is_sub = lhs.significand.sign() != rhs_sign * rhs.significand.sign();
419        let rnd_precision = self.precision + is_sub as usize;
420
421        let ediff = (rhs.exponent - lhs.exponent) as usize;
422        let rdigits = rhs.digits();
423        let ldigits_est = lhs.digits_ub();
424
425        // align the exponent
426        let low: (IBig, usize);
427        let (significand, exponent) =
428            if self.is_limited() && is_sub && ldigits_est + self.precision >= rdigits + ediff {
429                // Symmetric counterpart of the guard in `repr_add_large_small` (see there
430                // for the rationale); here the lower-exponent operand is `lhs`. Form the
431                // exact difference at full operand width and let the shared
432                // `repr_round_sum` round it once.
433                let rhs_signif = shl_digits::<B>(&rhs.significand, ediff);
434                low = (IBig::ZERO, 0);
435                (rhs_sign * rhs_signif + lhs.significand, lhs.exponent)
436            } else if self.is_limited()
437                && ldigits_est + 1 < ediff
438                && ldigits_est + 1 + rnd_precision < rdigits + ediff
439            {
440                // lhs is entirely below rhs's rounding window, so only its sign
441                // contributes; replace it with a unit sticky tail positioned at lhs's
442                // real magnitude (`ediff` digits below rhs's exponent). See
443                // `repr_add_large_small` for why the position must be `ediff` and not
444                // `precision - rdigits`.
445                low = (lhs.significand.signum(), ediff);
446                (rhs_sign * rhs.significand.clone(), rhs.exponent)
447            } else if self.is_limited() && rdigits >= self.precision {
448                // if the rhs already exceeds the desired precision, just align lhs
449                let (lhs_signif, r) = split_digits::<B>(lhs.significand, ediff);
450                low = (r, ediff);
451                match rhs_sign {
452                    Positive => (lhs_signif + &rhs.significand, rhs.exponent),
453                    Negative => (lhs_signif - &rhs.significand, rhs.exponent),
454                }
455            } else if self.is_limited() && ediff + rdigits > self.precision {
456                // if the shifted rhs exceeds the desired precision, align lhs and rhs to precision
457                let lshift = self.precision - rdigits;
458                let rshift = ediff - lshift;
459                let (lhs_signif, r) = split_digits::<B>(lhs.significand, rshift);
460                let rhs_signif = shl_digits::<B>(&rhs.significand, lshift);
461
462                low = (r, rshift);
463                (rhs_sign * rhs_signif + lhs_signif, rhs.exponent - lshift as isize)
464            } else {
465                // otherwise directly shift rhs to required position
466                let rhs_signif = shl_digits::<B>(&rhs.significand, ediff);
467                low = (IBig::ZERO, 0);
468                (rhs_sign * rhs_signif + lhs.significand, lhs.exponent)
469            };
470
471        self.repr_round_sum(significand, exponent, low, is_sub)
472    }
473
474    /// Add two floating point numbers under this context.
475    ///
476    /// # Examples
477    ///
478    /// ```
479    /// # use core::str::FromStr;
480    /// # use dashu_base::ParseError;
481    /// # use dashu_float::DBig;
482    /// use dashu_base::Approximation::*;
483    /// use dashu_float::{Context, round::{mode::HalfAway, Rounding::*}};
484    ///
485    /// let context = Context::<HalfAway>::new(2);
486    /// let a = DBig::from_str("1.234")?;
487    /// let b = DBig::from_str("6.789")?;
488    /// assert_eq!(context.add(&a.repr(), &b.repr()), Inexact(DBig::from_str("8.0")?, NoOp));
489    /// # Ok::<(), ParseError>(())
490    /// ```
491    pub fn add<const B: Word>(&self, lhs: &Repr<B>, rhs: &Repr<B>) -> Rounded<FBig<R, B>> {
492        assert_finite_operands(lhs, rhs);
493
494        let sum = if lhs.is_zero() {
495            self.repr_round_ref(rhs)
496        } else if rhs.is_zero() {
497            self.repr_round_ref(lhs)
498        } else {
499            match lhs.exponent.cmp(&rhs.exponent) {
500                Ordering::Equal => {
501                    self.repr_round(Repr::new(&lhs.significand + &rhs.significand, lhs.exponent))
502                }
503                Ordering::Greater => self.repr_add_large_small(lhs.clone(), rhs, Positive),
504                Ordering::Less => self.repr_add_small_large(lhs.clone(), rhs, Positive),
505            }
506        };
507        sum.map(|v| FBig::new(v, *self))
508    }
509
510    /// Subtract two floating point numbers under this context.
511    ///
512    /// # Examples
513    ///
514    /// ```
515    /// # use core::str::FromStr;
516    /// # use dashu_base::ParseError;
517    /// # use dashu_float::DBig;
518    /// use dashu_base::Approximation::*;
519    /// use dashu_float::{Context, round::{mode::HalfAway, Rounding::*}};
520    ///
521    /// let context = Context::<HalfAway>::new(2);
522    /// let a = DBig::from_str("1.234")?;
523    /// let b = DBig::from_str("6.789")?;
524    /// assert_eq!(
525    ///     context.sub(&a.repr(), &b.repr()),
526    ///     Inexact(DBig::from_str("-5.6")?, SubOne)
527    /// );
528    /// # Ok::<(), ParseError>(())
529    /// ```
530    pub fn sub<const B: Word>(&self, lhs: &Repr<B>, rhs: &Repr<B>) -> Rounded<FBig<R, B>> {
531        assert_finite_operands(lhs, rhs);
532
533        let sum = if lhs.is_zero() {
534            // Round `-rhs` directly rather than negating *after* rounding. For the asymmetric
535            // modes (Up = toward +∞, Down = toward −∞), `round(-x) != -round(x)`: rounding
536            // `rhs` toward +∞ then negating rounds in the wrong direction, so `0 - rhs`
537            // would land one ULP off (e.g. truncated instead of rounded away from the result).
538            self.repr_round_ref(&Repr::new(-&rhs.significand, rhs.exponent))
539        } else if rhs.is_zero() {
540            self.repr_round_ref(lhs)
541        } else {
542            match lhs.exponent.cmp(&rhs.exponent) {
543                Ordering::Equal => {
544                    self.repr_round(Repr::new(&lhs.significand - &rhs.significand, lhs.exponent))
545                }
546                Ordering::Greater => self.repr_add_large_small(lhs.clone(), rhs, Negative),
547                Ordering::Less => self.repr_add_small_large(lhs.clone(), rhs, Negative),
548            }
549        };
550        sum.map(|v| FBig::new(v, *self))
551    }
552}
553
554#[cfg(test)]
555mod tests {
556    use super::*;
557    use crate::round::mode::{HalfAway, HalfEven};
558
559    // Build a normalized Repr from a small integer significand and an exponent.
560    fn r<const B: Word>(sig: i128, exp: isize) -> Repr<B> {
561        Repr::new(IBig::from(sig), exp)
562    }
563
564    // Severe cancellation must not collapse a genuinely small difference to 0.
565    // Pristine returned `0` for the first two rows: the trimmed alignment path
566    // keeps only a bounded low tail and its single re-expand can't recover the
567    // lost leading digits.
568    #[test]
569    fn sub_severe_cancellation_decimal() {
570        let ctx = Context::<HalfAway>::new(3);
571        // 1.00 - 0.99999999 = 1e-8 (exactly representable at precision 3)
572        assert_eq!(
573            ctx.sub(&r::<10>(100, -2), &r::<10>(99999999, -8))
574                .value()
575                .repr(),
576            &r::<10>(1, -8)
577        );
578        // 1.00 - 0.99950001 = 4.9999e-4, rounds to 5.00e-4 (HalfAway)
579        assert_eq!(
580            ctx.sub(&r::<10>(100, -2), &r::<10>(99950001, -8))
581                .value()
582                .repr(),
583            &r::<10>(500, -6)
584        );
585    }
586
587    #[test]
588    fn sub_severe_cancellation_binary() {
589        let ctx = Context::<HalfEven>::new(10);
590        // 2^20 - (2^20 - 1) = 1, with the operands 20 exponent positions apart
591        assert_eq!(
592            ctx.sub(&r::<2>(1, 20), &r::<2>((1i128 << 20) - 1, 0))
593                .value()
594                .repr(),
595            &r::<2>(1, 0)
596        );
597        // same magnitude gap but the smaller-exponent operand is on the left
598        assert_eq!(
599            ctx.sub(&r::<2>((1i128 << 20) - 1, 0), &r::<2>(1, 20))
600                .value()
601                .repr(),
602            &r::<2>(-1, 0)
603        );
604        // 2^30 - (2^30 - 1) = 1
605        assert_eq!(
606            ctx.sub(&r::<2>(1, 30), &r::<2>((1i128 << 30) - 1, 0))
607                .value()
608                .repr(),
609            &r::<2>(1, 0)
610        );
611    }
612
613    // Effective subtraction reached through `Context::add` (opposite signs) must
614    // be fixed as well.
615    #[test]
616    fn add_effective_severe_cancellation() {
617        let ctx = Context::<HalfEven>::new(10);
618        // 2^20 + (-(2^20 - 1)) = 1
619        assert_eq!(
620            ctx.add(&r::<2>(1, 20), &r::<2>(-((1i128 << 20) - 1), 0))
621                .value()
622                .repr(),
623            &r::<2>(1, 0)
624        );
625    }
626
627    // The public operator path (`a - b`) routes through the same kernel.
628    #[test]
629    fn sub_operator_severe_cancellation() {
630        let a = FBig::<HalfEven, 2>::from_parts(IBig::from(1), 20);
631        let b = FBig::<HalfEven, 2>::from_parts(IBig::from((1i128 << 20) - 1), 0);
632        assert_eq!((a - b).repr(), &r::<2>(1, 0));
633    }
634
635    // Mild subtractions — the smaller operand stays below the larger's precision
636    // window — must keep their existing behavior and not be diverted to the
637    // full-width path.
638    #[test]
639    fn sub_mild_unchanged() {
640        let ctx = Context::<HalfAway>::new(3);
641        // 101 - 0.2 = 100.8, kept as 1008 * 10^-1 (one guard digit, as before)
642        assert_eq!(ctx.sub(&r::<10>(101, 0), &r::<10>(2, -1)).value().repr(), &r::<10>(1008, -1));
643    }
644
645    // Regression for the branch-1 signum-proxy bug (SUM-BUG.md §2c): when the larger
646    // operand has fewer digits than the precision and a negligible operand is added,
647    // the sticky proxy must be positioned at the operand's *real* magnitude (`ediff`),
648    // not at `precision - ldigits`. The old positioning let the re-expand drag the
649    // sticky up to the LSB, where for base 2 + round-to-nearest it equals exactly half
650    // and injected a spurious ULP: `1 + 2^-100` at precision 10 gave `513*2^-9` (=
651    // 1.00195…) instead of `1`.
652    #[test]
653    fn add_negligible_short_operand_no_spurious_ulp() {
654        // base 2 + HalfAway: the exact tie case
655        let ctx = Context::<HalfAway>::new(10);
656        assert_eq!(ctx.add(&r::<2>(1, 0), &r::<2>(1, -100)).value().repr(), &r::<2>(1, 0));
657        assert_eq!(ctx.sub(&r::<2>(1, 0), &r::<2>(1, -100)).value().repr(), &r::<2>(1, 0));
658        // larger short operand (digits < precision), negligible addend
659        let ctx = Context::<HalfAway>::new(50);
660        assert_eq!(
661            ctx.add(&r::<2>(0x12345, 0), &r::<2>(1, -200))
662                .value()
663                .repr(),
664            &r::<2>(0x12345, 0)
665        );
666        // base 10 was never affected (1 < ½·10), but check it stays correct
667        let ctx = Context::<HalfAway>::new(10);
668        assert_eq!(ctx.add(&r::<10>(1, 0), &r::<10>(1, -100)).value().repr(), &r::<10>(1, 0));
669    }
670}