1use super::addition::__add2;
2use super::{cmp_slice, BigUint};
3
4use crate::big_digit::{self, BigDigit, DoubleBigDigit};
5use crate::UsizePromotion;
6
7use core::cmp::Ordering::{Equal, Greater, Less};
8use core::mem;
9use core::ops::{Div, DivAssign, Rem, RemAssign};
10use num_integer::Integer;
11use num_traits::{CheckedDiv, CheckedEuclid, Euclid, One, ToPrimitive, Zero};
12
13pub(super) const FAST_DIV_WIDE: bool = cfg!(any(target_arch = "x86", target_arch = "x86_64"));
14
15/// Divide a two digit numerator by a one digit divisor, returns quotient and remainder:
16///
17/// Note: the caller must ensure that both the quotient and remainder will fit into a single digit.
18/// This is _not_ true for an arbitrary numerator/denominator.
19///
20/// (This function also matches what the x86 divide instruction does).
21#[cfg(any(miri, not(any(target_arch = "x86", target_arch = "x86_64"))))]
22#[inline]
23fn div_wide(hi: BigDigit, lo: BigDigit, divisor: BigDigit) -> (BigDigit, BigDigit) {
24 debug_assert!(hi < divisor);
25
26 let lhs = big_digit::to_doublebigdigit(hi, lo);
27 let rhs = DoubleBigDigit::from(divisor);
28 ((lhs / rhs) as BigDigit, (lhs % rhs) as BigDigit)
29}
30
31/// x86 and x86_64 can use a real `div` instruction.
32#[cfg(all(not(miri), any(target_arch = "x86", target_arch = "x86_64")))]
33#[inline]
34fn div_wide(hi: BigDigit, lo: BigDigit, divisor: BigDigit) -> (BigDigit, BigDigit) {
35 // This debug assertion covers the potential #DE for divisor==0 or a quotient too large for one
36 // register, otherwise in release mode it will become a target-specific fault like SIGFPE.
37 // This should never occur with the inputs from our few `div_wide` callers.
38 debug_assert!(hi < divisor);
39
40 // SAFETY: The `div` instruction only affects registers, reading the explicit operand as the
41 // divisor, and implicitly reading RDX:RAX or EDX:EAX as the dividend. The result is implicitly
42 // written back to RAX or EAX for the quotient and RDX or EDX for the remainder. No memory is
43 // used, and flags are not preserved.
44 unsafe {
45 let (div, rem);
46
47 cfg_digit!(
48 macro_rules! div {
49 () => {
50 "div {0:e}"
51 };
52 }
53 macro_rules! div {
54 () => {
55 "div {0:r}"
56 };
57 }
58 );
59
60 core::arch::asm!(
61 div!(),
62 in(reg) divisor,
63 inout("dx") hi => rem,
64 inout("ax") lo => div,
65 options(pure, nomem, nostack),
66 );
67
68 (div, rem)
69 }
70}
71
72/// For small divisors, we can divide without promoting to `DoubleBigDigit` by
73/// using half-size pieces of digit, like long-division.
74#[inline]
75fn div_half(rem: BigDigit, digit: BigDigit, divisor: BigDigit) -> (BigDigit, BigDigit) {
76 use crate::big_digit::{HALF, HALF_BITS};
77
78 debug_assert!(rem < divisor && divisor <= HALF);
79 let (hi: u64, rem: u64) = ((rem << HALF_BITS) | (digit >> HALF_BITS)).div_rem(&divisor);
80 let (lo: u64, rem: u64) = ((rem << HALF_BITS) | (digit & HALF)).div_rem(&divisor);
81 ((hi << HALF_BITS) | lo, rem)
82}
83
84#[inline]
85pub(super) fn div_rem_digit(mut a: BigUint, b: BigDigit) -> (BigUint, BigDigit) {
86 if b == 0 {
87 panic!("attempt to divide by zero")
88 }
89
90 let mut rem: u64 = 0;
91
92 if !FAST_DIV_WIDE && b <= big_digit::HALF {
93 for d: &mut u64 in a.data.iter_mut().rev() {
94 let (q: u64, r: u64) = div_half(rem, *d, divisor:b);
95 *d = q;
96 rem = r;
97 }
98 } else {
99 for d: &mut u64 in a.data.iter_mut().rev() {
100 let (q: u64, r: u64) = div_wide(hi:rem, *d, divisor:b);
101 *d = q;
102 rem = r;
103 }
104 }
105
106 (a.normalized(), rem)
107}
108
109#[inline]
110fn rem_digit(a: &BigUint, b: BigDigit) -> BigDigit {
111 if b == 0 {
112 panic!("attempt to divide by zero")
113 }
114
115 let mut rem: u64 = 0;
116
117 if !FAST_DIV_WIDE && b <= big_digit::HALF {
118 for &digit: u64 in a.data.iter().rev() {
119 let (_, r: u64) = div_half(rem, digit, divisor:b);
120 rem = r;
121 }
122 } else {
123 for &digit: u64 in a.data.iter().rev() {
124 let (_, r: u64) = div_wide(hi:rem, lo:digit, divisor:b);
125 rem = r;
126 }
127 }
128
129 rem
130}
131
132/// Subtract a multiple.
133/// a -= b * c
134/// Returns a borrow (if a < b then borrow > 0).
135fn sub_mul_digit_same_len(a: &mut [BigDigit], b: &[BigDigit], c: BigDigit) -> BigDigit {
136 debug_assert!(a.len() == b.len());
137
138 // carry is between -big_digit::MAX and 0, so to avoid overflow we store
139 // offset_carry = carry + big_digit::MAX
140 let mut offset_carry = big_digit::MAX;
141
142 for (x, y) in a.iter_mut().zip(b) {
143 // We want to calculate sum = x - y * c + carry.
144 // sum >= -(big_digit::MAX * big_digit::MAX) - big_digit::MAX
145 // sum <= big_digit::MAX
146 // Offsetting sum by (big_digit::MAX << big_digit::BITS) puts it in DoubleBigDigit range.
147 let offset_sum = big_digit::to_doublebigdigit(big_digit::MAX, *x)
148 - big_digit::MAX as DoubleBigDigit
149 + offset_carry as DoubleBigDigit
150 - *y as DoubleBigDigit * c as DoubleBigDigit;
151
152 let (new_offset_carry, new_x) = big_digit::from_doublebigdigit(offset_sum);
153 offset_carry = new_offset_carry;
154 *x = new_x;
155 }
156
157 // Return the borrow.
158 big_digit::MAX - offset_carry
159}
160
161fn div_rem(mut u: BigUint, mut d: BigUint) -> (BigUint, BigUint) {
162 if d.is_zero() {
163 panic!("attempt to divide by zero")
164 }
165 if u.is_zero() {
166 return (BigUint::ZERO, BigUint::ZERO);
167 }
168
169 if d.data.len() == 1 {
170 if d.data == [1] {
171 return (u, BigUint::ZERO);
172 }
173 let (div, rem) = div_rem_digit(u, d.data[0]);
174 // reuse d
175 d.data.clear();
176 d += rem;
177 return (div, d);
178 }
179
180 // Required or the q_len calculation below can underflow:
181 match u.cmp(&d) {
182 Less => return (BigUint::ZERO, u),
183 Equal => {
184 u.set_one();
185 return (u, BigUint::ZERO);
186 }
187 Greater => {} // Do nothing
188 }
189
190 // This algorithm is from Knuth, TAOCP vol 2 section 4.3, algorithm D:
191 //
192 // First, normalize the arguments so the highest bit in the highest digit of the divisor is
193 // set: the main loop uses the highest digit of the divisor for generating guesses, so we
194 // want it to be the largest number we can efficiently divide by.
195 //
196 let shift = d.data.last().unwrap().leading_zeros() as usize;
197
198 if shift == 0 {
199 // no need to clone d
200 div_rem_core(u, &d.data)
201 } else {
202 let (q, r) = div_rem_core(u << shift, &(d << shift).data);
203 // renormalize the remainder
204 (q, r >> shift)
205 }
206}
207
208pub(super) fn div_rem_ref(u: &BigUint, d: &BigUint) -> (BigUint, BigUint) {
209 if d.is_zero() {
210 panic!("attempt to divide by zero")
211 }
212 if u.is_zero() {
213 return (BigUint::ZERO, BigUint::ZERO);
214 }
215
216 if d.data.len() == 1 {
217 if d.data == [1] {
218 return (u.clone(), BigUint::ZERO);
219 }
220
221 let (div, rem) = div_rem_digit(u.clone(), d.data[0]);
222 return (div, rem.into());
223 }
224
225 // Required or the q_len calculation below can underflow:
226 match u.cmp(d) {
227 Less => return (BigUint::ZERO, u.clone()),
228 Equal => return (One::one(), BigUint::ZERO),
229 Greater => {} // Do nothing
230 }
231
232 // This algorithm is from Knuth, TAOCP vol 2 section 4.3, algorithm D:
233 //
234 // First, normalize the arguments so the highest bit in the highest digit of the divisor is
235 // set: the main loop uses the highest digit of the divisor for generating guesses, so we
236 // want it to be the largest number we can efficiently divide by.
237 //
238 let shift = d.data.last().unwrap().leading_zeros() as usize;
239
240 if shift == 0 {
241 // no need to clone d
242 div_rem_core(u.clone(), &d.data)
243 } else {
244 let (q, r) = div_rem_core(u << shift, &(d << shift).data);
245 // renormalize the remainder
246 (q, r >> shift)
247 }
248}
249
250/// An implementation of the base division algorithm.
251/// Knuth, TAOCP vol 2 section 4.3.1, algorithm D, with an improvement from exercises 19-21.
252fn div_rem_core(mut a: BigUint, b: &[BigDigit]) -> (BigUint, BigUint) {
253 debug_assert!(a.data.len() >= b.len() && b.len() > 1);
254 debug_assert!(b.last().unwrap().leading_zeros() == 0);
255
256 // The algorithm works by incrementally calculating "guesses", q0, for the next digit of the
257 // quotient. Once we have any number q0 such that (q0 << j) * b <= a, we can set
258 //
259 // q += q0 << j
260 // a -= (q0 << j) * b
261 //
262 // and then iterate until a < b. Then, (q, a) will be our desired quotient and remainder.
263 //
264 // q0, our guess, is calculated by dividing the last three digits of a by the last two digits of
265 // b - this will give us a guess that is close to the actual quotient, but is possibly greater.
266 // It can only be greater by 1 and only in rare cases, with probability at most
267 // 2^-(big_digit::BITS-1) for random a, see TAOCP 4.3.1 exercise 21.
268 //
269 // If the quotient turns out to be too large, we adjust it by 1:
270 // q -= 1 << j
271 // a += b << j
272
273 // a0 stores an additional extra most significant digit of the dividend, not stored in a.
274 let mut a0 = 0;
275
276 // [b1, b0] are the two most significant digits of the divisor. They never change.
277 let b0 = b[b.len() - 1];
278 let b1 = b[b.len() - 2];
279
280 let q_len = a.data.len() - b.len() + 1;
281 let mut q = BigUint {
282 data: vec![0; q_len],
283 };
284
285 for j in (0..q_len).rev() {
286 debug_assert!(a.data.len() == b.len() + j);
287
288 let a1 = *a.data.last().unwrap();
289 let a2 = a.data[a.data.len() - 2];
290
291 // The first q0 estimate is [a1,a0] / b0. It will never be too small, it may be too large
292 // by at most 2.
293 let (mut q0, mut r) = if a0 < b0 {
294 let (q0, r) = div_wide(a0, a1, b0);
295 (q0, r as DoubleBigDigit)
296 } else {
297 debug_assert!(a0 == b0);
298 // Avoid overflowing q0, we know the quotient fits in BigDigit.
299 // [a1,a0] = b0 * (1<<BITS - 1) + (a0 + a1)
300 (big_digit::MAX, a0 as DoubleBigDigit + a1 as DoubleBigDigit)
301 };
302
303 // r = [a1,a0] - q0 * b0
304 //
305 // Now we want to compute a more precise estimate [a2,a1,a0] / [b1,b0] which can only be
306 // less or equal to the current q0.
307 //
308 // q0 is too large if:
309 // [a2,a1,a0] < q0 * [b1,b0]
310 // (r << BITS) + a2 < q0 * b1
311 while r <= big_digit::MAX as DoubleBigDigit
312 && big_digit::to_doublebigdigit(r as BigDigit, a2)
313 < q0 as DoubleBigDigit * b1 as DoubleBigDigit
314 {
315 q0 -= 1;
316 r += b0 as DoubleBigDigit;
317 }
318
319 // q0 is now either the correct quotient digit, or in rare cases 1 too large.
320 // Subtract (q0 << j) from a. This may overflow, in which case we will have to correct.
321
322 let mut borrow = sub_mul_digit_same_len(&mut a.data[j..], b, q0);
323 if borrow > a0 {
324 // q0 is too large. We need to add back one multiple of b.
325 q0 -= 1;
326 borrow -= __add2(&mut a.data[j..], b);
327 }
328 // The top digit of a, stored in a0, has now been zeroed.
329 debug_assert!(borrow == a0);
330
331 q.data[j] = q0;
332
333 // Pop off the next top digit of a.
334 a0 = a.data.pop().unwrap();
335 }
336
337 a.data.push(a0);
338 a.normalize();
339
340 debug_assert_eq!(cmp_slice(&a.data, b), Less);
341
342 (q.normalized(), a)
343}
344
345forward_val_ref_binop!(impl Div for BigUint, div);
346forward_ref_val_binop!(impl Div for BigUint, div);
347forward_val_assign!(impl DivAssign for BigUint, div_assign);
348
349impl Div<BigUint> for BigUint {
350 type Output = BigUint;
351
352 #[inline]
353 fn div(self, other: BigUint) -> BigUint {
354 let (q: BigUint, _) = div_rem(self, d:other);
355 q
356 }
357}
358
359impl Div<&BigUint> for &BigUint {
360 type Output = BigUint;
361
362 #[inline]
363 fn div(self, other: &BigUint) -> BigUint {
364 let (q: BigUint, _) = self.div_rem(other);
365 q
366 }
367}
368impl DivAssign<&BigUint> for BigUint {
369 #[inline]
370 fn div_assign(&mut self, other: &BigUint) {
371 *self = &*self / other;
372 }
373}
374
375promote_unsigned_scalars!(impl Div for BigUint, div);
376promote_unsigned_scalars_assign!(impl DivAssign for BigUint, div_assign);
377forward_all_scalar_binop_to_val_val!(impl Div<u32> for BigUint, div);
378forward_all_scalar_binop_to_val_val!(impl Div<u64> for BigUint, div);
379forward_all_scalar_binop_to_val_val!(impl Div<u128> for BigUint, div);
380
381impl Div<u32> for BigUint {
382 type Output = BigUint;
383
384 #[inline]
385 fn div(self, other: u32) -> BigUint {
386 let (q: BigUint, _) = div_rem_digit(self, b:other as BigDigit);
387 q
388 }
389}
390impl DivAssign<u32> for BigUint {
391 #[inline]
392 fn div_assign(&mut self, other: u32) {
393 *self = &*self / other;
394 }
395}
396
397impl Div<BigUint> for u32 {
398 type Output = BigUint;
399
400 #[inline]
401 fn div(self, other: BigUint) -> BigUint {
402 match other.data.len() {
403 0 => panic!("attempt to divide by zero"),
404 1 => From::from(self as BigDigit / other.data[0]),
405 _ => BigUint::ZERO,
406 }
407 }
408}
409
410impl Div<u64> for BigUint {
411 type Output = BigUint;
412
413 #[inline]
414 fn div(self, other: u64) -> BigUint {
415 let (q: BigUint, _) = div_rem(self, d:From::from(other));
416 q
417 }
418}
419impl DivAssign<u64> for BigUint {
420 #[inline]
421 fn div_assign(&mut self, other: u64) {
422 // a vec of size 0 does not allocate, so this is fairly cheap
423 let temp: BigUint = mem::replace(self, Self::ZERO);
424 *self = temp / other;
425 }
426}
427
428impl Div<BigUint> for u64 {
429 type Output = BigUint;
430
431 cfg_digit!(
432 #[inline]
433 fn div(self, other: BigUint) -> BigUint {
434 match other.data.len() {
435 0 => panic!("attempt to divide by zero"),
436 1 => From::from(self / u64::from(other.data[0])),
437 2 => From::from(self / big_digit::to_doublebigdigit(other.data[1], other.data[0])),
438 _ => BigUint::ZERO,
439 }
440 }
441
442 #[inline]
443 fn div(self, other: BigUint) -> BigUint {
444 match other.data.len() {
445 0 => panic!("attempt to divide by zero"),
446 1 => From::from(self / other.data[0]),
447 _ => BigUint::ZERO,
448 }
449 }
450 );
451}
452
453impl Div<u128> for BigUint {
454 type Output = BigUint;
455
456 #[inline]
457 fn div(self, other: u128) -> BigUint {
458 let (q: BigUint, _) = div_rem(self, d:From::from(other));
459 q
460 }
461}
462
463impl DivAssign<u128> for BigUint {
464 #[inline]
465 fn div_assign(&mut self, other: u128) {
466 *self = &*self / other;
467 }
468}
469
470impl Div<BigUint> for u128 {
471 type Output = BigUint;
472
473 cfg_digit!(
474 #[inline]
475 fn div(self, other: BigUint) -> BigUint {
476 use super::u32_to_u128;
477 match other.data.len() {
478 0 => panic!("attempt to divide by zero"),
479 1 => From::from(self / u128::from(other.data[0])),
480 2 => From::from(
481 self / u128::from(big_digit::to_doublebigdigit(other.data[1], other.data[0])),
482 ),
483 3 => From::from(self / u32_to_u128(0, other.data[2], other.data[1], other.data[0])),
484 4 => From::from(
485 self / u32_to_u128(other.data[3], other.data[2], other.data[1], other.data[0]),
486 ),
487 _ => BigUint::ZERO,
488 }
489 }
490
491 #[inline]
492 fn div(self, other: BigUint) -> BigUint {
493 match other.data.len() {
494 0 => panic!("attempt to divide by zero"),
495 1 => From::from(self / other.data[0] as u128),
496 2 => From::from(self / big_digit::to_doublebigdigit(other.data[1], other.data[0])),
497 _ => BigUint::ZERO,
498 }
499 }
500 );
501}
502
503forward_val_ref_binop!(impl Rem for BigUint, rem);
504forward_ref_val_binop!(impl Rem for BigUint, rem);
505forward_val_assign!(impl RemAssign for BigUint, rem_assign);
506
507impl Rem<BigUint> for BigUint {
508 type Output = BigUint;
509
510 #[inline]
511 fn rem(self, other: BigUint) -> BigUint {
512 if let Some(other: u32) = other.to_u32() {
513 &self % other
514 } else {
515 let (_, r: BigUint) = div_rem(self, d:other);
516 r
517 }
518 }
519}
520
521impl Rem<&BigUint> for &BigUint {
522 type Output = BigUint;
523
524 #[inline]
525 fn rem(self, other: &BigUint) -> BigUint {
526 if let Some(other: u32) = other.to_u32() {
527 self % other
528 } else {
529 let (_, r: BigUint) = self.div_rem(other);
530 r
531 }
532 }
533}
534impl RemAssign<&BigUint> for BigUint {
535 #[inline]
536 fn rem_assign(&mut self, other: &BigUint) {
537 *self = &*self % other;
538 }
539}
540
541promote_unsigned_scalars!(impl Rem for BigUint, rem);
542promote_unsigned_scalars_assign!(impl RemAssign for BigUint, rem_assign);
543forward_all_scalar_binop_to_ref_val!(impl Rem<u32> for BigUint, rem);
544forward_all_scalar_binop_to_val_val!(impl Rem<u64> for BigUint, rem);
545forward_all_scalar_binop_to_val_val!(impl Rem<u128> for BigUint, rem);
546
547impl Rem<u32> for &BigUint {
548 type Output = BigUint;
549
550 #[inline]
551 fn rem(self, other: u32) -> BigUint {
552 rem_digit(self, b:other as BigDigit).into()
553 }
554}
555impl RemAssign<u32> for BigUint {
556 #[inline]
557 fn rem_assign(&mut self, other: u32) {
558 *self = &*self % other;
559 }
560}
561
562impl Rem<&BigUint> for u32 {
563 type Output = BigUint;
564
565 #[inline]
566 fn rem(mut self, other: &BigUint) -> BigUint {
567 self %= other;
568 From::from(self)
569 }
570}
571
572macro_rules! impl_rem_assign_scalar {
573 ($scalar:ty, $to_scalar:ident) => {
574 forward_val_assign_scalar!(impl RemAssign for BigUint, $scalar, rem_assign);
575 impl RemAssign<&BigUint> for $scalar {
576 #[inline]
577 fn rem_assign(&mut self, other: &BigUint) {
578 *self = match other.$to_scalar() {
579 None => *self,
580 Some(0) => panic!("attempt to divide by zero"),
581 Some(v) => *self % v
582 };
583 }
584 }
585 }
586}
587
588// we can scalar %= BigUint for any scalar, including signed types
589impl_rem_assign_scalar!(u128, to_u128);
590impl_rem_assign_scalar!(usize, to_usize);
591impl_rem_assign_scalar!(u64, to_u64);
592impl_rem_assign_scalar!(u32, to_u32);
593impl_rem_assign_scalar!(u16, to_u16);
594impl_rem_assign_scalar!(u8, to_u8);
595impl_rem_assign_scalar!(i128, to_i128);
596impl_rem_assign_scalar!(isize, to_isize);
597impl_rem_assign_scalar!(i64, to_i64);
598impl_rem_assign_scalar!(i32, to_i32);
599impl_rem_assign_scalar!(i16, to_i16);
600impl_rem_assign_scalar!(i8, to_i8);
601
602impl Rem<u64> for BigUint {
603 type Output = BigUint;
604
605 #[inline]
606 fn rem(self, other: u64) -> BigUint {
607 let (_, r: BigUint) = div_rem(self, d:From::from(other));
608 r
609 }
610}
611impl RemAssign<u64> for BigUint {
612 #[inline]
613 fn rem_assign(&mut self, other: u64) {
614 *self = &*self % other;
615 }
616}
617
618impl Rem<BigUint> for u64 {
619 type Output = BigUint;
620
621 #[inline]
622 fn rem(mut self, other: BigUint) -> BigUint {
623 self %= other;
624 From::from(self)
625 }
626}
627
628impl Rem<u128> for BigUint {
629 type Output = BigUint;
630
631 #[inline]
632 fn rem(self, other: u128) -> BigUint {
633 let (_, r: BigUint) = div_rem(self, d:From::from(other));
634 r
635 }
636}
637
638impl RemAssign<u128> for BigUint {
639 #[inline]
640 fn rem_assign(&mut self, other: u128) {
641 *self = &*self % other;
642 }
643}
644
645impl Rem<BigUint> for u128 {
646 type Output = BigUint;
647
648 #[inline]
649 fn rem(mut self, other: BigUint) -> BigUint {
650 self %= other;
651 From::from(self)
652 }
653}
654
655impl CheckedDiv for BigUint {
656 #[inline]
657 fn checked_div(&self, v: &BigUint) -> Option<BigUint> {
658 if v.is_zero() {
659 return None;
660 }
661 Some(self.div(v))
662 }
663}
664
665impl CheckedEuclid for BigUint {
666 #[inline]
667 fn checked_div_euclid(&self, v: &BigUint) -> Option<BigUint> {
668 if v.is_zero() {
669 return None;
670 }
671 Some(self.div_euclid(v))
672 }
673
674 #[inline]
675 fn checked_rem_euclid(&self, v: &BigUint) -> Option<BigUint> {
676 if v.is_zero() {
677 return None;
678 }
679 Some(self.rem_euclid(v))
680 }
681
682 fn checked_div_rem_euclid(&self, v: &Self) -> Option<(Self, Self)> {
683 Some(self.div_rem_euclid(v))
684 }
685}
686
687impl Euclid for BigUint {
688 #[inline]
689 fn div_euclid(&self, v: &BigUint) -> BigUint {
690 // trivially same as regular division
691 self / v
692 }
693
694 #[inline]
695 fn rem_euclid(&self, v: &BigUint) -> BigUint {
696 // trivially same as regular remainder
697 self % v
698 }
699
700 fn div_rem_euclid(&self, v: &Self) -> (Self, Self) {
701 // trivially same as regular division and remainder
702 self.div_rem(v)
703 }
704}
705