dashu_float/
div.rs

1use crate::{
2    error::{assert_finite_operands, assert_limited_precision},
3    fbig::FBig,
4    helper_macros::{self, impl_binop_assign_by_taking},
5    repr::{Context, Repr, Word},
6    round::{Round, Rounded},
7    utils::{digit_len, shl_digits_in_place, split_digits},
8};
9use core::ops::{Div, DivAssign, Rem, RemAssign};
10use dashu_base::{Approximation, DivEuclid, DivRem, DivRemEuclid, Inverse, RemEuclid};
11use dashu_int::{fast_div::ConstDivisor, modular::IntoRing, IBig, UBig};
12
13macro_rules! impl_div_or_rem_for_fbig {
14    (impl $op:ident, $method:ident, $repr_method:ident) => {
15        impl<R: Round, const B: Word> $op<FBig<R, B>> for FBig<R, B> {
16            type Output = FBig<R, B>;
17            fn $method(self, rhs: FBig<R, B>) -> Self::Output {
18                let context = Context::max(self.context, rhs.context);
19                FBig::new(context.$repr_method(self.repr, rhs.repr).value(), context)
20            }
21        }
22
23        impl<'l, R: Round, const B: Word> $op<FBig<R, B>> for &'l FBig<R, B> {
24            type Output = FBig<R, B>;
25            fn $method(self, rhs: FBig<R, B>) -> Self::Output {
26                let context = Context::max(self.context, rhs.context);
27                FBig::new(context.$repr_method(self.repr.clone(), rhs.repr).value(), context)
28            }
29        }
30
31        impl<'r, R: Round, const B: Word> $op<&'r FBig<R, B>> for FBig<R, B> {
32            type Output = FBig<R, B>;
33            fn $method(self, rhs: &FBig<R, B>) -> Self::Output {
34                let context = Context::max(self.context, rhs.context);
35                FBig::new(context.$repr_method(self.repr, rhs.repr.clone()).value(), context)
36            }
37        }
38
39        impl<'l, 'r, R: Round, const B: Word> $op<&'r FBig<R, B>> for &'l FBig<R, B> {
40            type Output = FBig<R, B>;
41            fn $method(self, rhs: &FBig<R, B>) -> Self::Output {
42                let context = Context::max(self.context, rhs.context);
43                FBig::new(
44                    context
45                        .$repr_method(self.repr.clone(), rhs.repr.clone())
46                        .value(),
47                    context,
48                )
49            }
50        }
51    };
52}
53impl_div_or_rem_for_fbig!(impl Div, div, repr_div);
54impl_div_or_rem_for_fbig!(impl Rem, rem, repr_rem);
55impl_binop_assign_by_taking!(impl DivAssign<Self>, div_assign, div);
56impl_binop_assign_by_taking!(impl RemAssign<Self>, rem_assign, rem);
57
58impl<R: Round, const B: Word> DivEuclid<FBig<R, B>> for FBig<R, B> {
59    type Output = IBig;
60    #[inline]
61    fn div_euclid(self, rhs: FBig<R, B>) -> Self::Output {
62        let (num, den) = align_as_int(self, rhs);
63        num.div_euclid(den)
64    }
65}
66
67impl<'l, R: Round, const B: Word> DivEuclid<FBig<R, B>> for &'l FBig<R, B> {
68    type Output = IBig;
69    #[inline]
70    fn div_euclid(self, rhs: FBig<R, B>) -> Self::Output {
71        self.clone().div_euclid(rhs)
72    }
73}
74
75impl<'r, R: Round, const B: Word> DivEuclid<&'r FBig<R, B>> for FBig<R, B> {
76    type Output = IBig;
77    #[inline]
78    fn div_euclid(self, rhs: &FBig<R, B>) -> Self::Output {
79        self.div_euclid(rhs.clone())
80    }
81}
82
83impl<'l, 'r, R: Round, const B: Word> DivEuclid<&'r FBig<R, B>> for &'l FBig<R, B> {
84    type Output = IBig;
85    #[inline]
86    fn div_euclid(self, rhs: &FBig<R, B>) -> Self::Output {
87        self.clone().div_euclid(rhs.clone())
88    }
89}
90
91impl<R: Round, const B: Word> RemEuclid<FBig<R, B>> for FBig<R, B> {
92    type Output = FBig<R, B>;
93    #[inline]
94    fn rem_euclid(self, rhs: FBig<R, B>) -> Self::Output {
95        let r_exponent = self.repr.exponent.min(rhs.repr.exponent);
96        let context = Context::max(self.context, rhs.context);
97
98        let (num, den) = align_as_int(self, rhs);
99        let r = num.rem_euclid(den);
100        let mut r = context.convert_int(r.into()).value();
101        if !r.repr.significand.is_zero() {
102            r.repr.exponent += r_exponent;
103        }
104        r
105    }
106}
107
108impl<'l, R: Round, const B: Word> RemEuclid<FBig<R, B>> for &'l FBig<R, B> {
109    type Output = FBig<R, B>;
110    #[inline]
111    fn rem_euclid(self, rhs: FBig<R, B>) -> Self::Output {
112        self.clone().rem_euclid(rhs)
113    }
114}
115
116impl<'r, R: Round, const B: Word> RemEuclid<&'r FBig<R, B>> for FBig<R, B> {
117    type Output = FBig<R, B>;
118    #[inline]
119    fn rem_euclid(self, rhs: &FBig<R, B>) -> Self::Output {
120        self.rem_euclid(rhs.clone())
121    }
122}
123
124impl<'l, 'r, R: Round, const B: Word> RemEuclid<&'r FBig<R, B>> for &'l FBig<R, B> {
125    type Output = FBig<R, B>;
126    #[inline]
127    fn rem_euclid(self, rhs: &FBig<R, B>) -> Self::Output {
128        self.clone().rem_euclid(rhs.clone())
129    }
130}
131
132impl<R: Round, const B: Word> DivRemEuclid<FBig<R, B>> for FBig<R, B> {
133    type OutputDiv = IBig;
134    type OutputRem = FBig<R, B>;
135    #[inline]
136    fn div_rem_euclid(self, rhs: FBig<R, B>) -> (IBig, FBig<R, B>) {
137        let r_exponent = self.repr.exponent.min(rhs.repr.exponent);
138        let context = Context::max(self.context, rhs.context);
139
140        let (num, den) = align_as_int(self, rhs);
141        let (q, r) = num.div_rem_euclid(den);
142        let mut r = context.convert_int(r.into()).value();
143        if !r.repr.significand.is_zero() {
144            r.repr.exponent += r_exponent;
145        }
146        (q, r)
147    }
148}
149
150impl<'l, R: Round, const B: Word> DivRemEuclid<FBig<R, B>> for &'l FBig<R, B> {
151    type OutputDiv = IBig;
152    type OutputRem = FBig<R, B>;
153    #[inline]
154    fn div_rem_euclid(self, rhs: FBig<R, B>) -> (IBig, FBig<R, B>) {
155        self.clone().div_rem_euclid(rhs)
156    }
157}
158
159impl<'r, R: Round, const B: Word> DivRemEuclid<&'r FBig<R, B>> for FBig<R, B> {
160    type OutputDiv = IBig;
161    type OutputRem = FBig<R, B>;
162    #[inline]
163    fn div_rem_euclid(self, rhs: &FBig<R, B>) -> (IBig, FBig<R, B>) {
164        self.div_rem_euclid(rhs.clone())
165    }
166}
167
168impl<'l, 'r, R: Round, const B: Word> DivRemEuclid<&'r FBig<R, B>> for &'l FBig<R, B> {
169    type OutputDiv = IBig;
170    type OutputRem = FBig<R, B>;
171    #[inline]
172    fn div_rem_euclid(self, rhs: &FBig<R, B>) -> (IBig, FBig<R, B>) {
173        self.clone().div_rem_euclid(rhs.clone())
174    }
175}
176
177macro_rules! impl_div_primitive_with_fbig {
178    ($($t:ty)*) => {$(
179        helper_macros::impl_binop_with_primitive!(impl Div<$t>, div);
180        helper_macros::impl_binop_assign_with_primitive!(impl DivAssign<$t>, div_assign);
181    )*};
182}
183impl_div_primitive_with_fbig!(u8 u16 u32 u64 u128 usize UBig i8 i16 i32 i64 i128 isize IBig);
184// TODO: we should specialize FBig / UBig or FBig / IBig for better efficiency
185
186impl<R: Round, const B: Word> Inverse for FBig<R, B> {
187    type Output = FBig<R, B>;
188
189    #[inline]
190    fn inv(self) -> Self::Output {
191        self.context.inv(&self.repr).value()
192    }
193}
194
195impl<R: Round, const B: Word> Inverse for &FBig<R, B> {
196    type Output = FBig<R, B>;
197
198    #[inline]
199    fn inv(self) -> Self::Output {
200        self.context.inv(&self.repr).value()
201    }
202}
203
204// Align two float by exponent such that they are both turned into integers
205fn align_as_int<R: Round, const B: Word>(lhs: FBig<R, B>, rhs: FBig<R, B>) -> (IBig, IBig) {
206    let ediff = lhs.repr.exponent - rhs.repr.exponent;
207    let (mut num, mut den) = (lhs.repr.significand, rhs.repr.significand);
208    if ediff >= 0 {
209        shl_digits_in_place::<B>(&mut num, ediff as _);
210    } else {
211        shl_digits_in_place::<B>(&mut den, (-ediff) as _);
212    }
213    (num, den)
214}
215
216impl<R: Round> Context<R> {
217    pub(crate) fn repr_div<const B: Word>(&self, lhs: Repr<B>, rhs: Repr<B>) -> Rounded<Repr<B>> {
218        assert_finite_operands(&lhs, &rhs);
219        assert_limited_precision(self.precision);
220
221        // this method don't deal with the case where lhs significand is too large
222        debug_assert!(lhs.digits() <= self.precision + rhs.digits());
223
224        let (mut q, mut r) = lhs.significand.div_rem(&rhs.significand);
225        let mut e = lhs.exponent - rhs.exponent;
226        if r.is_zero() {
227            return Approximation::Exact(Repr::new(q, e));
228        }
229
230        let ddigits = digit_len::<B>(&rhs.significand);
231        if q.is_zero() {
232            // lhs.significand < rhs.significand
233            let rdigits = digit_len::<B>(&r); // rdigits <= ddigits
234            let shift = ddigits + self.precision - rdigits;
235            shl_digits_in_place::<B>(&mut r, shift);
236            e -= shift as isize;
237            let (q0, r0) = r.div_rem(&rhs.significand);
238            q = q0;
239            r = r0;
240        } else {
241            let ndigits = digit_len::<B>(&q) + ddigits;
242            if ndigits < ddigits + self.precision {
243                // TODO: here the operations can be optimized: 1. prevent double power, 2. q += q0 can be |= if B is power of 2
244                let shift = ddigits + self.precision - ndigits;
245                shl_digits_in_place::<B>(&mut q, shift);
246                shl_digits_in_place::<B>(&mut r, shift);
247                e -= shift as isize;
248
249                let (q0, r0) = r.div_rem(&rhs.significand);
250                q += q0;
251                r = r0;
252            }
253        }
254
255        if r.is_zero() {
256            Approximation::Exact(Repr::new(q, e))
257        } else {
258            let adjust = R::round_ratio(&q, r, &rhs.significand);
259            Approximation::Inexact(Repr::new(q + adjust, e), adjust)
260        }
261    }
262
263    pub(crate) fn repr_rem<const B: Word>(&self, lhs: Repr<B>, rhs: Repr<B>) -> Rounded<Repr<B>> {
264        assert_finite_operands(&lhs, &rhs);
265
266        let (lhs_sign, lhs_signif) = lhs.significand.into_parts();
267        let (_, rhs_signif) = rhs.significand.into_parts();
268
269        use core::cmp::Ordering;
270        let significand = match lhs.exponent.cmp(&rhs.exponent) {
271            Ordering::Equal => {
272                let r1 = lhs_signif % &rhs_signif;
273                let r2 = rhs_signif - &r1;
274                if r1 < r2 {
275                    IBig::from_parts(lhs_sign, r1)
276                } else {
277                    IBig::from_parts(-lhs_sign, r2)
278                }
279            }
280            Ordering::Greater => {
281                // if the least significant digit of lhs is higher than rhs, then we can
282                // align lhs to rhs and do simple modulo operations
283                let modulo = ConstDivisor::new(rhs_signif);
284                let shift = (lhs.exponent - rhs.exponent) as usize;
285                let scaling = if B == 2 {
286                    (UBig::ONE << shift).into_ring(&modulo)
287                } else {
288                    UBig::from_word(B).into_ring(&modulo).pow(&shift.into())
289                };
290                let r = lhs_signif.into_ring(&modulo) * scaling;
291                let r1 = r.residue();
292                let r2 = (-r).residue();
293                if r1 < r2 {
294                    IBig::from_parts(lhs_sign, r1)
295                } else {
296                    IBig::from_parts(-lhs_sign, r2)
297                }
298            }
299            Ordering::Less => {
300                // otherwise we have to split lhs into two parts
301                let shift = (rhs.exponent - lhs.exponent) as usize;
302                let (hi, lo) = split_digits::<B>(lhs_signif.into(), shift);
303
304                let mut r1 = hi % &rhs_signif;
305                let mut r2 = rhs_signif - &r1;
306
307                shl_digits_in_place::<B>(&mut r1, shift);
308                r1 += &lo;
309
310                shl_digits_in_place::<B>(&mut r2, shift);
311                r2 -= lo;
312
313                if r1 < r2 {
314                    lhs_sign * r1
315                } else {
316                    (-lhs_sign) * r2
317                }
318            }
319        };
320
321        let exponent = lhs.exponent.min(rhs.exponent);
322        if significand.is_zero() {
323            Approximation::Exact(Repr::zero())
324        } else {
325            self.repr_round(Repr::new(significand, exponent))
326        }
327    }
328
329    /// Divide two floating point numbers under this context.
330    ///
331    /// # Examples
332    ///
333    /// ```
334    /// # use core::str::FromStr;
335    /// # use dashu_base::ParseError;
336    /// # use dashu_float::DBig;
337    /// use dashu_base::Approximation::*;
338    /// use dashu_float::{Context, round::{mode::HalfAway, Rounding::*}};
339    ///
340    /// let context = Context::<HalfAway>::new(2);
341    /// let a = DBig::from_str("-1.234")?;
342    /// let b = DBig::from_str("6.789")?;
343    /// assert_eq!(context.div(&a.repr(), &b.repr()), Inexact(DBig::from_str("-0.18")?, NoOp));
344    /// # Ok::<(), ParseError>(())
345    /// ```
346    ///
347    /// # Euclidean Division
348    ///
349    /// To do euclidean division on the float numbers (get an integer quotient and remainder, equivalent to C99's
350    /// `fmod` and `remquo`), please use the methods provided by traits [DivEuclid], [RemEuclid] and [DivRemEuclid].
351    ///
352    pub fn div<const B: Word>(&self, lhs: &Repr<B>, rhs: &Repr<B>) -> Rounded<FBig<R, B>> {
353        assert_finite_operands(lhs, rhs);
354
355        let lhs_repr = if !lhs.is_zero() && lhs.digits_ub() > rhs.digits_lb() + self.precision {
356            // shrink lhs if it's larger than necessary
357            Self::new(rhs.digits() + self.precision)
358                .repr_round_ref(lhs)
359                .value()
360        } else {
361            lhs.clone()
362        };
363        self.repr_div(lhs_repr, rhs.clone())
364            .map(|v| FBig::new(v, *self))
365    }
366
367    /// Calculate the remainder of `⌈lhs / rhs⌋`.
368    ///
369    /// The remainder is calculated as `r = lhs - ⌈lhs / rhs⌋ * rhs`, the division rounds to the nearest and ties to away.
370    /// So if `n = (lhs / rhs).round()`, then `lhs == n * rhs + r` (given enough precision).
371    ///
372    /// # Examples
373    ///
374    /// ```
375    /// # use core::str::FromStr;
376    /// # use dashu_base::ParseError;
377    /// # use dashu_float::DBig;
378    /// use dashu_base::Approximation::*;
379    /// use dashu_float::{Context, round::{mode::HalfAway, Rounding::*}};
380    ///
381    /// let context = Context::<HalfAway>::new(3);
382    /// let a = DBig::from_str("6.789")?;
383    /// let b = DBig::from_str("-1.234")?;
384    /// assert_eq!(context.rem(&a.repr(), &b.repr()), Exact(DBig::from_str("-0.615")?));
385    /// # Ok::<(), ParseError>(())
386    /// ```
387    pub fn rem<const B: Word>(&self, lhs: &Repr<B>, rhs: &Repr<B>) -> Rounded<FBig<R, B>> {
388        assert_finite_operands(lhs, rhs);
389        self.repr_rem(lhs.clone(), rhs.clone())
390            .map(|v| FBig::new(v, *self))
391    }
392
393    /// Compute the multiplicative inverse of an `FBig`
394    ///
395    /// ```
396    /// # use core::str::FromStr;
397    /// # use dashu_base::ParseError;
398    /// # use dashu_float::DBig;
399    /// use dashu_base::Approximation::*;
400    /// use dashu_float::{Context, round::{mode::HalfAway, Rounding::*}};
401    ///
402    /// let context = Context::<HalfAway>::new(2);
403    /// let a = DBig::from_str("-1.234")?;
404    /// assert_eq!(context.inv(&a.repr()), Inexact(DBig::from_str("-0.81")?, NoOp));
405    /// # Ok::<(), ParseError>(())
406    /// ```
407    #[inline]
408    pub fn inv<const B: Word>(&self, f: &Repr<B>) -> Rounded<FBig<R, B>> {
409        self.repr_div(Repr::one(), f.clone())
410            .map(|v| FBig::new(v, *self))
411    }
412}
dashu_float/div.rs

dashu_float/
div.rs
