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 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 return Rounded::Exact(Repr::new(significand, exponent));
230 }
231
232 let rnd_precision = self.precision + is_sub as usize;
234
235 let digits = digit_len::<B>(&significand);
237 match digits.cmp(&rnd_precision) {
238 Ordering::Equal => {}
239 Ordering::Greater => {
240 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 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 if low.0.is_zero() {
281 Rounded::Exact(Repr::new(significand, exponent))
282 } else {
283 let adjust = R::round_fract::<B>(&significand, low.0, low.1);
286 Rounded::Inexact(Repr::new(significand + adjust, exponent), adjust)
287 }
288 }
289
290 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 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(); let low: (IBig, usize); let (significand, exponent) =
310 if self.is_limited() && is_sub && rdigits_est + self.precision >= ldigits + ediff {
311 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 low = (rhs_sign * rhs.significand.signum(), ediff);
347 (lhs.significand, lhs.exponent)
348 } else if self.is_limited() && ldigits >= self.precision {
349 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 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 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 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 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 let low: (IBig, usize);
427 let (significand, exponent) =
428 if self.is_limited() && is_sub && ldigits_est + self.precision >= rdigits + ediff {
429 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 low = (lhs.significand.signum(), ediff);
446 (rhs_sign * rhs.significand.clone(), rhs.exponent)
447 } else if self.is_limited() && rdigits >= self.precision {
448 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 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 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 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 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 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 fn r<const B: Word>(sig: i128, exp: isize) -> Repr<B> {
561 Repr::new(IBig::from(sig), exp)
562 }
563
564 #[test]
569 fn sub_severe_cancellation_decimal() {
570 let ctx = Context::<HalfAway>::new(3);
571 assert_eq!(
573 ctx.sub(&r::<10>(100, -2), &r::<10>(99999999, -8))
574 .value()
575 .repr(),
576 &r::<10>(1, -8)
577 );
578 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 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 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 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 #[test]
616 fn add_effective_severe_cancellation() {
617 let ctx = Context::<HalfEven>::new(10);
618 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 #[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 #[test]
639 fn sub_mild_unchanged() {
640 let ctx = Context::<HalfAway>::new(3);
641 assert_eq!(ctx.sub(&r::<10>(101, 0), &r::<10>(2, -1)).value().repr(), &r::<10>(1008, -1));
643 }
644
645 #[test]
653 fn add_negligible_short_operand_no_spurious_ulp() {
654 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 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 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}