1 //! Building-blocks for arbitrary-precision math.
2 //!
3 //! These algorithms assume little-endian order for the large integer
4 //! buffers, so for a `vec![0, 1, 2, 3]`, `3` is the most significant limb,
5 //! and `0` is the least significant limb.
6
7 // We have a lot of higher-order, efficient math routines that may
8 // come in handy later. These aren't trivial to implement, so keep them.
9 #![allow(dead_code)]
10
11 use crate::util::*;
12
13 // ALIASES
14 // -------
15
16 // Type for a single limb of the big integer.
17 //
18 // A limb is analogous to a digit in base10, except, it stores 32-bit
19 // or 64-bit numbers instead.
20 //
21 // This should be all-known 64-bit platforms supported by Rust.
22 // https://forge.rust-lang.org/platform-support.html
23 //
24 // Platforms where native 128-bit multiplication is explicitly supported:
25 // - x86_64 (Supported via `MUL`).
26 // - mips64 (Supported via `DMULTU`, which `HI` and `LO` can be read-from).
27 //
28 // Platforms where native 64-bit multiplication is supported and
29 // you can extract hi-lo for 64-bit multiplications.
30 // aarch64 (Requires `UMULH` and `MUL` to capture high and low bits).
31 // powerpc64 (Requires `MULHDU` and `MULLD` to capture high and low bits).
32 //
33 // Platforms where native 128-bit multiplication is not supported,
34 // requiring software emulation.
35 // sparc64 (`UMUL` only supported double-word arguments).
36 cfg_if! {
37 if #[cfg(limb_width_64)] {
38 pub type Limb = u64;
39 type Wide = u128;
40 type SignedWide = i128;
41 } else {
42 pub type Limb = u32;
43 type Wide = u64;
44 type SignedWide = i64;
45 }} // cfg_if
46
47 perftools_inline!{
48 /// Cast to limb type.
49 pub(super) fn as_limb<T: Integer>(t: T)
50 -> Limb
51 {
52 as_cast(t)
53 }}
54
55 perftools_inline!{
56 /// Cast to wide type.
57 fn as_wide<T: Integer>(t: T)
58 -> Wide
59 {
60 as_cast(t)
61 }}
62
63 perftools_inline!{
64 /// Cast tosigned wide type.
65 fn as_signed_wide<T: Integer>(t: T)
66 -> SignedWide
67 {
68 as_cast(t)
69 }}
70
71 // SPLIT
72 // -----
73
74 perftools_inline!{
75 /// Split u16 into limbs, in little-endian order.
76 fn split_u16(x: u16) -> [Limb; 1] {
77 [as_limb(x)]
78 }}
79
80 perftools_inline!{
81 /// Split u32 into limbs, in little-endian order.
82 fn split_u32(x: u32) -> [Limb; 1] {
83 [as_limb(x)]
84 }}
85
86 perftools_inline!{
87 /// Split u64 into limbs, in little-endian order.
88 #[cfg(limb_width_32)]
89 fn split_u64(x: u64) -> [Limb; 2] {
90 [as_limb(x), as_limb(x >> 32)]
91 }}
92
93 perftools_inline!{
94 /// Split u64 into limbs, in little-endian order.
95 #[cfg(limb_width_64)]
96 fn split_u64(x: u64) -> [Limb; 1] {
97 [as_limb(x)]
98 }}
99
100 perftools_inline!{
101 /// Split u128 into limbs, in little-endian order.
102 #[cfg(limb_width_32)]
103 fn split_u128(x: u128) -> [Limb; 4] {
104 [as_limb(x), as_limb(x >> 32), as_limb(x >> 64), as_limb(x >> 96)]
105 }}
106
107 perftools_inline!{
108 /// Split u128 into limbs, in little-endian order.
109 #[cfg(limb_width_64)]
110 fn split_u128(x: u128) -> [Limb; 2] {
111 [as_limb(x), as_limb(x >> 64)]
112 }}
113
114 // HI BITS
115 // -------
116
117 // NONZERO
118
119 perftools_inline!{
120 /// Check if any of the remaining bits are non-zero.
121 pub fn nonzero<T: Integer>(x: &[T], rindex: usize) -> bool {
122 let len = x.len();
123 let slc = &x[..len-rindex];
124 slc.iter().rev().any(|&x| x != T::ZERO)
125 }}
126
127 // HI16
128
129 perftools_inline!{
130 /// Shift 16-bit integer to high 16-bits.
131 fn u16_to_hi16_1(r0: u16) -> (u16, bool) {
132 debug_assert!(r0 != 0);
133 let ls = r0.leading_zeros();
134 (r0 << ls, false)
135 }}
136
137 perftools_inline!{
138 /// Shift 2 16-bit integers to high 16-bits.
139 fn u16_to_hi16_2(r0: u16, r1: u16) -> (u16, bool) {
140 debug_assert!(r0 != 0);
141 let ls = r0.leading_zeros();
142 let rs = 16 - ls;
143 let v = match ls {
144 0 => r0,
145 _ => (r0 << ls) | (r1 >> rs),
146 };
147 let n = r1 << ls != 0;
148 (v, n)
149 }}
150
151 perftools_inline!{
152 /// Shift 32-bit integer to high 16-bits.
153 fn u32_to_hi16_1(r0: u32) -> (u16, bool) {
154 let r0 = u32_to_hi32_1(r0).0;
155 ((r0 >> 16).as_u16(), r0.as_u16() != 0)
156 }}
157
158 perftools_inline!{
159 /// Shift 2 32-bit integers to high 16-bits.
160 fn u32_to_hi16_2(r0: u32, r1: u32) -> (u16, bool) {
161 let (r0, n) = u32_to_hi32_2(r0, r1);
162 ((r0 >> 16).as_u16(), n || r0.as_u16() != 0)
163 }}
164
165 perftools_inline!{
166 /// Shift 64-bit integer to high 16-bits.
167 fn u64_to_hi16_1(r0: u64) -> (u16, bool) {
168 let r0 = u64_to_hi64_1(r0).0;
169 ((r0 >> 48).as_u16(), r0.as_u16() != 0)
170 }}
171
172 perftools_inline!{
173 /// Shift 2 64-bit integers to high 16-bits.
174 fn u64_to_hi16_2(r0: u64, r1: u64) -> (u16, bool) {
175 let (r0, n) = u64_to_hi64_2(r0, r1);
176 ((r0 >> 48).as_u16(), n || r0.as_u16() != 0)
177 }}
178
179 /// Trait to export the high 16-bits from a little-endian slice.
180 trait Hi16<T>: SliceLike<T> {
181 /// Get the hi16 bits from a 1-limb slice.
hi16_1(&self) -> (u16, bool)182 fn hi16_1(&self) -> (u16, bool);
183
184 /// Get the hi16 bits from a 2-limb slice.
hi16_2(&self) -> (u16, bool)185 fn hi16_2(&self) -> (u16, bool);
186
187 perftools_inline!{
188 /// High-level exporter to extract the high 16 bits from a little-endian slice.
189 fn hi16(&self) -> (u16, bool) {
190 match self.len() {
191 0 => (0, false),
192 1 => self.hi16_1(),
193 _ => self.hi16_2(),
194 }
195 }}
196 }
197
198 impl Hi16<u16> for [u16] {
199 perftools_inline!{
200 fn hi16_1(&self) -> (u16, bool) {
201 debug_assert!(self.len() == 1);
202 let rview = self.rview();
203 let r0 = rview[0];
204 u16_to_hi16_1(r0)
205 }}
206
207 perftools_inline!{
208 fn hi16_2(&self) -> (u16, bool) {
209 debug_assert!(self.len() == 2);
210 let rview = self.rview();
211 let r0 = rview[0];
212 let r1 = rview[1];
213 let (v, n) = u16_to_hi16_2(r0, r1);
214 (v, n || nonzero(self, 2))
215 }}
216 }
217
218 impl Hi16<u32> for [u32] {
219 perftools_inline!{
220 fn hi16_1(&self) -> (u16, bool) {
221 debug_assert!(self.len() == 1);
222 let rview = self.rview();
223 let r0 = rview[0];
224 u32_to_hi16_1(r0)
225 }}
226
227 perftools_inline!{
228 fn hi16_2(&self) -> (u16, bool) {
229 debug_assert!(self.len() == 2);
230 let rview = self.rview();
231 let r0 = rview[0];
232 let r1 = rview[1];
233 let (v, n) = u32_to_hi16_2(r0, r1);
234 (v, n || nonzero(self, 2))
235 }}
236 }
237
238 impl Hi16<u64> for [u64] {
239 perftools_inline!{
240 fn hi16_1(&self) -> (u16, bool) {
241 debug_assert!(self.len() == 1);
242 let rview = self.rview();
243 let r0 = rview[0];
244 u64_to_hi16_1(r0)
245 }}
246
247 perftools_inline!{
248 fn hi16_2(&self) -> (u16, bool) {
249 debug_assert!(self.len() == 2);
250 let rview = self.rview();
251 let r0 = rview[0];
252 let r1 = rview[1];
253 let (v, n) = u64_to_hi16_2(r0, r1);
254 (v, n || nonzero(self, 2))
255 }}
256 }
257
258 // HI32
259
260 perftools_inline!{
261 /// Shift 32-bit integer to high 32-bits.
262 fn u32_to_hi32_1(r0: u32) -> (u32, bool) {
263 debug_assert!(r0 != 0);
264 let ls = r0.leading_zeros();
265 (r0 << ls, false)
266 }}
267
268 perftools_inline!{
269 /// Shift 2 32-bit integers to high 32-bits.
270 fn u32_to_hi32_2(r0: u32, r1: u32) -> (u32, bool) {
271 debug_assert!(r0 != 0);
272 let ls = r0.leading_zeros();
273 let rs = 32 - ls;
274 let v = match ls {
275 0 => r0,
276 _ => (r0 << ls) | (r1 >> rs),
277 };
278 let n = r1 << ls != 0;
279 (v, n)
280 }}
281
282 perftools_inline!{
283 /// Shift 64-bit integer to high 32-bits.
284 fn u64_to_hi32_1(r0: u64) -> (u32, bool) {
285 let r0 = u64_to_hi64_1(r0).0;
286 ((r0 >> 32).as_u32(), r0.as_u32() != 0)
287 }}
288
289 perftools_inline!{
290 /// Shift 2 64-bit integers to high 32-bits.
291 fn u64_to_hi32_2(r0: u64, r1: u64) -> (u32, bool) {
292 let (r0, n) = u64_to_hi64_2(r0, r1);
293 ((r0 >> 32).as_u32(), n || r0.as_u32() != 0)
294 }}
295
296 /// Trait to export the high 32-bits from a little-endian slice.
297 trait Hi32<T>: SliceLike<T> {
298 /// Get the hi32 bits from a 1-limb slice.
hi32_1(&self) -> (u32, bool)299 fn hi32_1(&self) -> (u32, bool);
300
301 /// Get the hi32 bits from a 2-limb slice.
hi32_2(&self) -> (u32, bool)302 fn hi32_2(&self) -> (u32, bool);
303
304 /// Get the hi32 bits from a 3-limb slice.
hi32_3(&self) -> (u32, bool)305 fn hi32_3(&self) -> (u32, bool);
306
307 perftools_inline!{
308 /// High-level exporter to extract the high 32 bits from a little-endian slice.
309 fn hi32(&self) -> (u32, bool) {
310 match self.len() {
311 0 => (0, false),
312 1 => self.hi32_1(),
313 2 => self.hi32_2(),
314 _ => self.hi32_3(),
315 }
316 }}
317 }
318
319 impl Hi32<u16> for [u16] {
320 perftools_inline!{
321 fn hi32_1(&self) -> (u32, bool) {
322 debug_assert!(self.len() == 1);
323 let rview = self.rview();
324 u32_to_hi32_1(rview[0].as_u32())
325 }}
326
327 perftools_inline!{
328 fn hi32_2(&self) -> (u32, bool) {
329 debug_assert!(self.len() == 2);
330 let rview = self.rview();
331 let r0 = rview[0].as_u32() << 16;
332 let r1 = rview[1].as_u32();
333 u32_to_hi32_1(r0 | r1)
334 }}
335
336 perftools_inline!{
337 fn hi32_3(&self) -> (u32, bool) {
338 debug_assert!(self.len() >= 3);
339 let rview = self.rview();
340 let r0 = rview[0].as_u32();
341 let r1 = rview[1].as_u32() << 16;
342 let r2 = rview[2].as_u32();
343 let (v, n) = u32_to_hi32_2( r0, r1 | r2);
344 (v, n || nonzero(self, 3))
345 }}
346 }
347
348 impl Hi32<u32> for [u32] {
349 perftools_inline!{
350 fn hi32_1(&self) -> (u32, bool) {
351 debug_assert!(self.len() == 1);
352 let rview = self.rview();
353 let r0 = rview[0];
354 u32_to_hi32_1(r0)
355 }}
356
357 perftools_inline!{
358 fn hi32_2(&self) -> (u32, bool) {
359 debug_assert!(self.len() >= 2);
360 let rview = self.rview();
361 let r0 = rview[0];
362 let r1 = rview[1];
363 let (v, n) = u32_to_hi32_2(r0, r1);
364 (v, n || nonzero(self, 2))
365 }}
366
367 perftools_inline!{
368 fn hi32_3(&self) -> (u32, bool) {
369 self.hi32_2()
370 }}
371 }
372
373 impl Hi32<u64> for [u64] {
374 perftools_inline!{
375 fn hi32_1(&self) -> (u32, bool) {
376 debug_assert!(self.len() == 1);
377 let rview = self.rview();
378 let r0 = rview[0];
379 u64_to_hi32_1(r0)
380 }}
381
382 perftools_inline!{
383 fn hi32_2(&self) -> (u32, bool) {
384 debug_assert!(self.len() >= 2);
385 let rview = self.rview();
386 let r0 = rview[0];
387 let r1 = rview[1];
388 let (v, n) = u64_to_hi32_2(r0, r1);
389 (v, n || nonzero(self, 2))
390 }}
391
392 perftools_inline!{
393 fn hi32_3(&self) -> (u32, bool) {
394 self.hi32_2()
395 }}
396 }
397
398 // HI64
399
400 perftools_inline!{
401 /// Shift 64-bit integer to high 64-bits.
402 fn u64_to_hi64_1(r0: u64) -> (u64, bool) {
403 debug_assert!(r0 != 0);
404 let ls = r0.leading_zeros();
405 (r0 << ls, false)
406 }}
407
408 perftools_inline!{
409 /// Shift 2 64-bit integers to high 64-bits.
410 fn u64_to_hi64_2(r0: u64, r1: u64) -> (u64, bool) {
411 debug_assert!(r0 != 0);
412 let ls = r0.leading_zeros();
413 let rs = 64 - ls;
414 let v = match ls {
415 0 => r0,
416 _ => (r0 << ls) | (r1 >> rs),
417 };
418 let n = r1 << ls != 0;
419 (v, n)
420 }}
421
422 /// Trait to export the high 64-bits from a little-endian slice.
423 trait Hi64<T>: SliceLike<T> {
424 /// Get the hi64 bits from a 1-limb slice.
hi64_1(&self) -> (u64, bool)425 fn hi64_1(&self) -> (u64, bool);
426
427 /// Get the hi64 bits from a 2-limb slice.
hi64_2(&self) -> (u64, bool)428 fn hi64_2(&self) -> (u64, bool);
429
430 /// Get the hi64 bits from a 3-limb slice.
hi64_3(&self) -> (u64, bool)431 fn hi64_3(&self) -> (u64, bool);
432
433 /// Get the hi64 bits from a 4-limb slice.
hi64_4(&self) -> (u64, bool)434 fn hi64_4(&self) -> (u64, bool);
435
436 /// Get the hi64 bits from a 5-limb slice.
hi64_5(&self) -> (u64, bool)437 fn hi64_5(&self) -> (u64, bool);
438
439 perftools_inline!{
440 /// High-level exporter to extract the high 64 bits from a little-endian slice.
441 fn hi64(&self) -> (u64, bool) {
442 match self.len() {
443 0 => (0, false),
444 1 => self.hi64_1(),
445 2 => self.hi64_2(),
446 3 => self.hi64_3(),
447 4 => self.hi64_4(),
448 _ => self.hi64_5(),
449 }
450 }}
451 }
452
453 impl Hi64<u16> for [u16] {
454 perftools_inline!{
455 fn hi64_1(&self) -> (u64, bool) {
456 debug_assert!(self.len() == 1);
457 let rview = self.rview();
458 let r0 = rview[0].as_u64();
459 u64_to_hi64_1(r0)
460 }}
461
462 perftools_inline!{
463 fn hi64_2(&self) -> (u64, bool) {
464 debug_assert!(self.len() == 2);
465 let rview = self.rview();
466 let r0 = rview[0].as_u64() << 16;
467 let r1 = rview[1].as_u64();
468 u64_to_hi64_1(r0 | r1)
469 }}
470
471 perftools_inline!{
472 fn hi64_3(&self) -> (u64, bool) {
473 debug_assert!(self.len() == 3);
474 let rview = self.rview();
475 let r0 = rview[0].as_u64() << 32;
476 let r1 = rview[1].as_u64() << 16;
477 let r2 = rview[2].as_u64();
478 u64_to_hi64_1(r0 | r1 | r2)
479 }}
480
481 perftools_inline!{
482 fn hi64_4(&self) -> (u64, bool) {
483 debug_assert!(self.len() == 4);
484 let rview = self.rview();
485 let r0 = rview[0].as_u64() << 48;
486 let r1 = rview[1].as_u64() << 32;
487 let r2 = rview[2].as_u64() << 16;
488 let r3 = rview[3].as_u64();
489 u64_to_hi64_1(r0 | r1 | r2 | r3)
490 }}
491
492 perftools_inline!{
493 fn hi64_5(&self) -> (u64, bool) {
494 debug_assert!(self.len() >= 5);
495 let rview = self.rview();
496 let r0 = rview[0].as_u64();
497 let r1 = rview[1].as_u64() << 48;
498 let r2 = rview[2].as_u64() << 32;
499 let r3 = rview[3].as_u64() << 16;
500 let r4 = rview[4].as_u64();
501 let (v, n) = u64_to_hi64_2(r0, r1 | r2 | r3 | r4);
502 (v, n || nonzero(self, 5))
503 }}
504 }
505
506 impl Hi64<u32> for [u32] {
507 perftools_inline!{
508 fn hi64_1(&self) -> (u64, bool) {
509 debug_assert!(self.len() == 1);
510 let rview = self.rview();
511 let r0 = rview[0].as_u64();
512 u64_to_hi64_1(r0)
513 }}
514
515 perftools_inline!{
516 fn hi64_2(&self) -> (u64, bool) {
517 debug_assert!(self.len() == 2);
518 let rview = self.rview();
519 let r0 = rview[0].as_u64() << 32;
520 let r1 = rview[1].as_u64();
521 u64_to_hi64_1(r0 | r1)
522 }}
523
524 perftools_inline!{
525 fn hi64_3(&self) -> (u64, bool) {
526 debug_assert!(self.len() >= 3);
527 let rview = self.rview();
528 let r0 = rview[0].as_u64();
529 let r1 = rview[1].as_u64() << 32;
530 let r2 = rview[2].as_u64();
531 let (v, n) = u64_to_hi64_2(r0, r1 | r2);
532 (v, n || nonzero(self, 3))
533 }}
534
535 perftools_inline!{
536 fn hi64_4(&self) -> (u64, bool) {
537 self.hi64_3()
538 }}
539
540 perftools_inline!{
541 fn hi64_5(&self) -> (u64, bool) {
542 self.hi64_3()
543 }}
544 }
545
546 impl Hi64<u64> for [u64] {
547 perftools_inline!{
548 fn hi64_1(&self) -> (u64, bool) {
549 debug_assert!(self.len() == 1);
550 let rview = self.rview();
551 let r0 = rview[0];
552 u64_to_hi64_1(r0)
553 }}
554
555 perftools_inline!{
556 fn hi64_2(&self) -> (u64, bool) {
557 debug_assert!(self.len() >= 2);
558 let rview = self.rview();
559 let r0 = rview[0];
560 let r1 = rview[1];
561 let (v, n) = u64_to_hi64_2(r0, r1);
562 (v, n || nonzero(self, 2))
563 }}
564
565 perftools_inline!{
566 fn hi64_3(&self) -> (u64, bool) {
567 self.hi64_2()
568 }}
569
570 perftools_inline!{
571 fn hi64_4(&self) -> (u64, bool) {
572 self.hi64_2()
573 }}
574
575 perftools_inline!{
576 fn hi64_5(&self) -> (u64, bool) {
577 self.hi64_2()
578 }}
579 }
580
581 // HI128
582
583 perftools_inline!{
584 /// Shift 128-bit integer to high 128-bits.
585 fn u128_to_hi128_1(r0: u128) -> (u128, bool) {
586 let ls = r0.leading_zeros();
587 (r0 << ls, false)
588 }}
589
590 perftools_inline!{
591 /// Shift 2 128-bit integers to high 128-bits.
592 fn u128_to_hi128_2(r0: u128, r1: u128) -> (u128, bool) {
593 let ls = r0.leading_zeros();
594 let rs = 128 - ls;
595 let v = (r0 << ls) | (r1 >> rs);
596 let n = r1 << ls != 0;
597 (v, n)
598 }}
599
600 /// Trait to export the high 128-bits from a little-endian slice.
601 trait Hi128<T>: SliceLike<T> {
602 /// Get the hi128 bits from a 1-limb slice.
hi128_1(&self) -> (u128, bool)603 fn hi128_1(&self) -> (u128, bool);
604
605 /// Get the hi128 bits from a 2-limb slice.
hi128_2(&self) -> (u128, bool)606 fn hi128_2(&self) -> (u128, bool);
607
608 /// Get the hi128 bits from a 3-limb slice.
hi128_3(&self) -> (u128, bool)609 fn hi128_3(&self) -> (u128, bool);
610
611 /// Get the hi128 bits from a 4-limb slice.
hi128_4(&self) -> (u128, bool)612 fn hi128_4(&self) -> (u128, bool);
613
614 /// Get the hi128 bits from a 5-limb slice.
hi128_5(&self) -> (u128, bool)615 fn hi128_5(&self) -> (u128, bool);
616
617 /// Get the hi128 bits from a 5-limb slice.
hi128_6(&self) -> (u128, bool)618 fn hi128_6(&self) -> (u128, bool);
619
620 /// Get the hi128 bits from a 5-limb slice.
hi128_7(&self) -> (u128, bool)621 fn hi128_7(&self) -> (u128, bool);
622
623 /// Get the hi128 bits from a 5-limb slice.
hi128_8(&self) -> (u128, bool)624 fn hi128_8(&self) -> (u128, bool);
625
626 /// Get the hi128 bits from a 5-limb slice.
hi128_9(&self) -> (u128, bool)627 fn hi128_9(&self) -> (u128, bool);
628
629 perftools_inline!{
630 /// High-level exporter to extract the high 128 bits from a little-endian slice.
631 fn hi128(&self) -> (u128, bool) {
632 match self.len() {
633 0 => (0, false),
634 1 => self.hi128_1(),
635 2 => self.hi128_2(),
636 3 => self.hi128_3(),
637 4 => self.hi128_4(),
638 6 => self.hi128_6(),
639 7 => self.hi128_7(),
640 8 => self.hi128_8(),
641 _ => self.hi128_9(),
642 }
643 }}
644 }
645
646 impl Hi128<u16> for [u16] {
647 perftools_inline!{
648 fn hi128_1(&self) -> (u128, bool) {
649 debug_assert!(self.len() == 1);
650 let rview = self.rview();
651 let r0 = rview[0].as_u128();
652 u128_to_hi128_1(r0)
653 }}
654
655 perftools_inline!{
656 fn hi128_2(&self) -> (u128, bool) {
657 debug_assert!(self.len() == 2);
658 let rview = self.rview();
659 let r0 = rview[0].as_u128() << 16;
660 let r1 = rview[1].as_u128();
661 u128_to_hi128_1(r0 | r1)
662 }}
663
664 perftools_inline!{
665 fn hi128_3(&self) -> (u128, bool) {
666 debug_assert!(self.len() == 3);
667 let rview = self.rview();
668 let r0 = rview[0].as_u128() << 32;
669 let r1 = rview[1].as_u128() << 16;
670 let r2 = rview[2].as_u128();
671 u128_to_hi128_1(r0 | r1 | r2)
672 }}
673
674 perftools_inline!{
675 fn hi128_4(&self) -> (u128, bool) {
676 debug_assert!(self.len() == 4);
677 let rview = self.rview();
678 let r0 = rview[0].as_u128() << 48;
679 let r1 = rview[1].as_u128() << 32;
680 let r2 = rview[2].as_u128() << 16;
681 let r3 = rview[3].as_u128();
682 u128_to_hi128_1(r0 | r1 | r2 | r3)
683 }}
684
685 perftools_inline!{
686 fn hi128_5(&self) -> (u128, bool) {
687 debug_assert!(self.len() == 5);
688 let rview = self.rview();
689 let r0 = rview[0].as_u128() << 64;
690 let r1 = rview[1].as_u128() << 48;
691 let r2 = rview[2].as_u128() << 32;
692 let r3 = rview[3].as_u128() << 16;
693 let r4 = rview[4].as_u128();
694 u128_to_hi128_1(r0 | r1 | r2 | r3 | r4)
695 }}
696
697 perftools_inline!{
698 fn hi128_6(&self) -> (u128, bool) {
699 debug_assert!(self.len() == 6);
700 let rview = self.rview();
701 let r0 = rview[0].as_u128() << 80;
702 let r1 = rview[1].as_u128() << 64;
703 let r2 = rview[2].as_u128() << 48;
704 let r3 = rview[3].as_u128() << 32;
705 let r4 = rview[4].as_u128() << 16;
706 let r5 = rview[5].as_u128();
707 u128_to_hi128_1(r0 | r1 | r2 | r3 | r4 | r5)
708 }}
709
710 perftools_inline!{
711 fn hi128_7(&self) -> (u128, bool) {
712 debug_assert!(self.len() == 7);
713 let rview = self.rview();
714 let r0 = rview[0].as_u128() << 96;
715 let r1 = rview[1].as_u128() << 80;
716 let r2 = rview[2].as_u128() << 64;
717 let r3 = rview[3].as_u128() << 48;
718 let r4 = rview[4].as_u128() << 32;
719 let r5 = rview[5].as_u128() << 16;
720 let r6 = rview[6].as_u128();
721 u128_to_hi128_1(r0 | r1 | r2 | r3 | r4 | r5 | r6)
722 }}
723
724 perftools_inline!{
725 fn hi128_8(&self) -> (u128, bool) {
726 debug_assert!(self.len() == 8);
727 let rview = self.rview();
728 let r0 = rview[0].as_u128() << 112;
729 let r1 = rview[1].as_u128() << 96;
730 let r2 = rview[2].as_u128() << 80;
731 let r3 = rview[3].as_u128() << 64;
732 let r4 = rview[4].as_u128() << 48;
733 let r5 = rview[5].as_u128() << 32;
734 let r6 = rview[6].as_u128() << 16;
735 let r7 = rview[7].as_u128();
736 u128_to_hi128_1(r0 | r1 | r2 | r3 | r4 | r5 | r6 | r7)
737 }}
738
739 perftools_inline!{
740 fn hi128_9(&self) -> (u128, bool) {
741 debug_assert!(self.len() >= 9);
742 let rview = self.rview();
743 let r0 = rview[0].as_u128();
744 let r1 = rview[1].as_u128() << 112;
745 let r2 = rview[2].as_u128() << 96;
746 let r3 = rview[3].as_u128() << 80;
747 let r4 = rview[4].as_u128() << 64;
748 let r5 = rview[5].as_u128() << 48;
749 let r6 = rview[6].as_u128() << 32;
750 let r7 = rview[7].as_u128() << 16;
751 let r8 = rview[8].as_u128();
752 let (v, n) = u128_to_hi128_2(r0, r1 | r2 | r3 | r4 | r5 | r6 | r7 | r8);
753 (v, n || nonzero(self, 9))
754 }}
755 }
756
757 impl Hi128<u32> for [u32] {
758 perftools_inline!{
759 fn hi128_1(&self) -> (u128, bool) {
760 debug_assert!(self.len() == 1);
761 let rview = self.rview();
762 let r0 = rview[0].as_u128();
763 u128_to_hi128_1(r0)
764 }}
765
766 perftools_inline!{
767 fn hi128_2(&self) -> (u128, bool) {
768 debug_assert!(self.len() == 2);
769 let rview = self.rview();
770 let r0 = rview[0].as_u128() << 32;
771 let r1 = rview[1].as_u128();
772 u128_to_hi128_1(r0 | r1)
773 }}
774
775 perftools_inline!{
776 fn hi128_3(&self) -> (u128, bool) {
777 debug_assert!(self.len() == 3);
778 let rview = self.rview();
779 let r0 = rview[0].as_u128() << 64;
780 let r1 = rview[1].as_u128() << 32;
781 let r2 = rview[2].as_u128();
782 u128_to_hi128_1(r0 | r1 | r2)
783 }}
784
785 perftools_inline!{
786 fn hi128_4(&self) -> (u128, bool) {
787 debug_assert!(self.len() == 4);
788 let rview = self.rview();
789 let r0 = rview[0].as_u128() << 96;
790 let r1 = rview[1].as_u128() << 64;
791 let r2 = rview[2].as_u128() << 32;
792 let r3 = rview[3].as_u128();
793 u128_to_hi128_1(r0 | r1 | r2 | r3)
794 }}
795
796 perftools_inline!{
797 fn hi128_5(&self) -> (u128, bool) {
798 debug_assert!(self.len() >= 5);
799 let rview = self.rview();
800 let r0 = rview[0].as_u128();
801 let r1 = rview[1].as_u128() << 96;
802 let r2 = rview[2].as_u128() << 64;
803 let r3 = rview[3].as_u128() << 32;
804 let r4 = rview[4].as_u128();
805 let (v, n) = u128_to_hi128_2(r0, r1 | r2 | r3 | r4);
806 (v, n || nonzero(self, 5))
807 }}
808
809 perftools_inline!{
810 fn hi128_6(&self) -> (u128, bool) {
811 self.hi128_5()
812 }}
813
814 perftools_inline!{
815 fn hi128_7(&self) -> (u128, bool) {
816 self.hi128_5()
817 }}
818
819 perftools_inline!{
820 fn hi128_8(&self) -> (u128, bool) {
821 self.hi128_5()
822 }}
823
824 perftools_inline!{
825 fn hi128_9(&self) -> (u128, bool) {
826 self.hi128_5()
827 }}
828 }
829
830 impl Hi128<u64> for [u64] {
831 perftools_inline!{
832 fn hi128_1(&self) -> (u128, bool) {
833 debug_assert!(self.len() == 1);
834 let rview = self.rview();
835 let r0 = rview[0].as_u128();
836 u128_to_hi128_1(r0)
837 }}
838
839 perftools_inline!{
840 fn hi128_2(&self) -> (u128, bool) {
841 debug_assert!(self.len() == 2);
842 let rview = self.rview();
843 let r0 = rview[0].as_u128() << 64;
844 let r1 = rview[1].as_u128();
845 u128_to_hi128_1(r0 | r1)
846 }}
847
848 perftools_inline!{
849 fn hi128_3(&self) -> (u128, bool) {
850 debug_assert!(self.len() >= 3);
851 let rview = self.rview();
852 let r0 = rview[0].as_u128();
853 let r1 = rview[1].as_u128() << 64;
854 let r2 = rview[2].as_u128();
855 let (v, n) = u128_to_hi128_2(r0, r1 | r2);
856 (v, n || nonzero(self, 3))
857 }}
858
859 perftools_inline!{
860 fn hi128_4(&self) -> (u128, bool) {
861 self.hi128_3()
862 }}
863
864 perftools_inline!{
865 fn hi128_5(&self) -> (u128, bool) {
866 self.hi128_3()
867 }}
868
869 perftools_inline!{
870 fn hi128_6(&self) -> (u128, bool) {
871 self.hi128_3()
872 }}
873
874 perftools_inline!{
875 fn hi128_7(&self) -> (u128, bool) {
876 self.hi128_3()
877 }}
878
879 perftools_inline!{
880 fn hi128_8(&self) -> (u128, bool) {
881 self.hi128_3()
882 }}
883
884 perftools_inline!{
885 fn hi128_9(&self) -> (u128, bool) {
886 self.hi128_3()
887 }}
888 }
889
890 // SCALAR
891 // ------
892
893 // Scalar-to-scalar operations, for building-blocks for arbitrary-precision
894 // operations.
895
896 pub(in crate::atof::algorithm) mod scalar {
897
898 use super::*;
899
900 // ADDITION
901
902 perftools_inline!{
903 /// Add two small integers and return the resulting value and if overflow happens.
904 pub fn add(x: Limb, y: Limb)
905 -> (Limb, bool)
906 {
907 x.overflowing_add(y)
908 }}
909
910 perftools_inline!{
911 /// AddAssign two small integers and return if overflow happens.
912 pub fn iadd(x: &mut Limb, y: Limb)
913 -> bool
914 {
915 let t = add(*x, y);
916 *x = t.0;
917 t.1
918 }}
919
920 // SUBTRACTION
921
922 perftools_inline!{
923 /// Subtract two small integers and return the resulting value and if overflow happens.
924 pub fn sub(x: Limb, y: Limb)
925 -> (Limb, bool)
926 {
927 x.overflowing_sub(y)
928 }}
929
930 perftools_inline!{
931 /// SubAssign two small integers and return if overflow happens.
932 pub fn isub(x: &mut Limb, y: Limb)
933 -> bool
934 {
935 let t = sub(*x, y);
936 *x = t.0;
937 t.1
938 }}
939
940 // MULTIPLICATION
941
942 perftools_inline!{
943 /// Multiply two small integers (with carry) (and return the overflow contribution).
944 ///
945 /// Returns the (low, high) components.
946 pub fn mul(x: Limb, y: Limb, carry: Limb)
947 -> (Limb, Limb)
948 {
949 // Cannot overflow, as long as wide is 2x as wide. This is because
950 // the following is always true:
951 // `Wide::max_value() - (Narrow::max_value() * Narrow::max_value()) >= Narrow::max_value()`
952 let z: Wide = as_wide(x) * as_wide(y) + as_wide(carry);
953 (as_limb(z), as_limb(z >> <Limb as Integer>::BITS))
954 }}
955
956 perftools_inline!{
957 /// Multiply two small integers (with carry) (and return if overflow happens).
958 pub fn imul(x: &mut Limb, y: Limb, carry: Limb)
959 -> Limb
960 {
961 let t = mul(*x, y, carry);
962 *x = t.0;
963 t.1
964 }}
965
966 // DIVISION
967
968 perftools_inline!{
969 /// Divide two small integers (with remainder) (and return the remainder contribution).
970 ///
971 /// Returns the (value, remainder) components.
972 pub fn div(x: Limb, y: Limb, rem: Limb)
973 -> (Limb, Limb)
974 {
975 // Cannot overflow, as long as wide is 2x as wide.
976 let x = as_wide(x) | (as_wide(rem) << <Limb as Integer>::BITS);
977 let y = as_wide(y);
978 (as_limb(x / y), as_limb(x % y))
979 }}
980
981 perftools_inline!{
982 /// DivAssign two small integers and return the remainder.
983 pub fn idiv(x: &mut Limb, y: Limb, rem: Limb)
984 -> Limb
985 {
986 let t = div(*x, y, rem);
987 *x = t.0;
988 t.1
989 }}
990
991 } // scalar
992
993 // SMALL
994 // -----
995
996 // Large-to-small operations, to modify a big integer from a native scalar.
997
998 pub(in crate::atof::algorithm) mod small {
999
1000 use crate::lib::iter;
1001 use super::*;
1002 use super::super::small_powers::*;
1003 use super::super::large_powers::*;
1004
1005 // PROPERTIES
1006
1007 perftools_inline!{
1008 /// Get the number of leading zero values in the storage.
1009 /// Assumes the value is normalized.
1010 pub fn leading_zero_limbs(_: &[Limb]) -> usize {
1011 0
1012 }}
1013
1014 perftools_inline!{
1015 /// Get the number of trailing zero values in the storage.
1016 /// Assumes the value is normalized.
1017 pub fn trailing_zero_limbs(x: &[Limb]) -> usize {
1018 let mut iter = x.iter().enumerate();
1019 let opt = iter.find(|&tup| !tup.1.is_zero());
1020 let value = opt
1021 .map(|t| t.0)
1022 .unwrap_or(x.len());
1023
1024 value
1025 }}
1026
1027 perftools_inline!{
1028 /// Get number of leading zero bits in the storage.
1029 pub fn leading_zeros(x: &[Limb]) -> usize {
1030 if x.is_empty() {
1031 0
1032 } else {
1033 x.rindex(0).leading_zeros().as_usize()
1034 }
1035 }}
1036
1037 perftools_inline!{
1038 /// Get number of trailing zero bits in the storage.
1039 /// Assumes the value is normalized.
1040 pub fn trailing_zeros(x: &[Limb]) -> usize {
1041 // Get the index of the last non-zero value
1042 let index = trailing_zero_limbs(x);
1043 let mut count = index.saturating_mul(<Limb as Integer>::BITS);
1044 if let Some(value) = x.get(index) {
1045 count = count.saturating_add(value.trailing_zeros().as_usize());
1046 }
1047 count
1048 }}
1049
1050 // BIT LENGTH
1051
1052 perftools_inline!{
1053 /// Calculate the bit-length of the big-integer.
1054 pub fn bit_length(x: &[Limb]) -> usize {
1055 // Avoid overflowing, calculate via total number of bits
1056 // minus leading zero bits.
1057 let nlz = leading_zeros(x);
1058 <Limb as Integer>::BITS.checked_mul(x.len())
1059 .map(|v| v - nlz)
1060 .unwrap_or(usize::max_value())
1061 }}
1062
1063 // BIT LENGTH
1064
1065 perftools_inline!{
1066 /// Calculate the limb-length of the big-integer.
1067 pub fn limb_length(x: &[Limb]) -> usize {
1068 x.len()
1069 }}
1070
1071 // SHR
1072
1073 perftools_inline!{
1074 /// Shift-right bits inside a buffer and returns the truncated bits.
1075 ///
1076 /// Returns the truncated bits.
1077 ///
1078 /// Assumes `n < <Limb as Integer>::BITS`, IE, internally shifting bits.
1079 pub fn ishr_bits<T>(x: &mut T, n: usize)
1080 -> Limb
1081 where T: CloneableVecLike<Limb>
1082 {
1083 // Need to shift by the number of `bits % <Limb as Integer>::BITS`.
1084 let bits = <Limb as Integer>::BITS;
1085 debug_assert!(n < bits && n != 0);
1086
1087 // Internally, for each item, we shift left by n, and add the previous
1088 // right shifted limb-bits.
1089 // For example, we transform (for u8) shifted right 2, to:
1090 // b10100100 b01000010
1091 // b101001 b00010000
1092 let lshift = bits - n;
1093 let rshift = n;
1094 let mut prev: Limb = 0;
1095 for xi in x.iter_mut().rev() {
1096 let tmp = *xi;
1097 *xi >>= rshift;
1098 *xi |= prev << lshift;
1099 prev = tmp;
1100 }
1101
1102 prev & lower_n_mask(as_limb(rshift))
1103 }}
1104
1105 perftools_inline!{
1106 /// Shift-right `n` limbs inside a buffer and returns if all the truncated limbs are zero.
1107 ///
1108 /// Assumes `n` is not 0.
1109 pub fn ishr_limbs<T>(x: &mut T, n: usize)
1110 -> bool
1111 where T: CloneableVecLike<Limb>
1112 {
1113 debug_assert!(n != 0);
1114
1115 if n >= x.len() {
1116 x.clear();
1117 false
1118 } else {
1119 let is_zero = (&x[..n]).iter().all(|v| v.is_zero());
1120 x.remove_many(0..n);
1121 is_zero
1122 }
1123 }}
1124
1125 perftools_inline!{
1126 /// Shift-left buffer by n bits and return if we should round-up.
1127 pub fn ishr<T>(x: &mut T, n: usize)
1128 -> bool
1129 where T: CloneableVecLike<Limb>
1130 {
1131 let bits = <Limb as Integer>::BITS;
1132 // Need to pad with zeros for the number of `bits / <Limb as Integer>::BITS`,
1133 // and shift-left with carry for `bits % <Limb as Integer>::BITS`.
1134 let rem = n % bits;
1135 let div = n / bits;
1136 let is_zero = match div.is_zero() {
1137 true => true,
1138 false => ishr_limbs(x, div),
1139 };
1140 let truncated = match rem.is_zero() {
1141 true => 0,
1142 false => ishr_bits(x, rem),
1143 };
1144
1145 // Calculate if we need to roundup.
1146 let roundup = {
1147 let halfway = lower_n_halfway(as_limb(rem));
1148 if truncated > halfway {
1149 // Above halfway
1150 true
1151 } else if truncated == halfway {
1152 // Exactly halfway, if !is_zero, we have a tie-breaker,
1153 // otherwise, we follow round-to-nearest, tie-even rules.
1154 // Cannot be empty, since truncated is non-zero.
1155 !is_zero || x[0].is_odd()
1156 } else {
1157 // Below halfway
1158 false
1159 }
1160 };
1161
1162 // Normalize the data
1163 normalize(x);
1164
1165 roundup
1166 }}
1167
1168 perftools_inline!{
1169 /// Shift-left buffer by n bits.
1170 pub fn shr<T>(x: &[Limb], n: usize)
1171 -> (T, bool)
1172 where T: CloneableVecLike<Limb>
1173 {
1174 let mut z = T::default();
1175 z.extend_from_slice(x);
1176 let roundup = ishr(&mut z, n);
1177 (z, roundup)
1178 }}
1179
1180 // SHL
1181
1182 perftools_inline!{
1183 /// Shift-left bits inside a buffer.
1184 ///
1185 /// Assumes `n < <Limb as Integer>::BITS`, IE, internally shifting bits.
1186 pub fn ishl_bits<T>(x: &mut T, n: usize)
1187 where T: CloneableVecLike<Limb>
1188 {
1189 // Need to shift by the number of `bits % <Limb as Integer>::BITS)`.
1190 let bits = <Limb as Integer>::BITS;
1191 debug_assert!(n < bits);
1192 if n.is_zero() {
1193 return;
1194 }
1195
1196 // Internally, for each item, we shift left by n, and add the previous
1197 // right shifted limb-bits.
1198 // For example, we transform (for u8) shifted left 2, to:
1199 // b10100100 b01000010
1200 // b10 b10010001 b00001000
1201 let rshift = bits - n;
1202 let lshift = n;
1203 let mut prev: Limb = 0;
1204 for xi in x.iter_mut() {
1205 let tmp = *xi;
1206 *xi <<= lshift;
1207 *xi |= prev >> rshift;
1208 prev = tmp;
1209 }
1210
1211 // Always push the carry, even if it creates a non-normal result.
1212 let carry = prev >> rshift;
1213 if carry != 0 {
1214 x.push(carry);
1215 }
1216 }}
1217
1218 perftools_inline!{
1219 /// Shift-left bits inside a buffer.
1220 ///
1221 /// Assumes `n < <Limb as Integer>::BITS`, IE, internally shifting bits.
1222 pub fn shl_bits<T>(x: &[Limb], n: usize)
1223 -> T
1224 where T: CloneableVecLike<Limb>
1225 {
1226 let mut z = T::default();
1227 z.extend_from_slice(x);
1228 ishl_bits(&mut z, n);
1229 z
1230 }}
1231
1232 perftools_inline!{
1233 /// Shift-left `n` digits inside a buffer.
1234 ///
1235 /// Assumes `n` is not 0.
1236 pub fn ishl_limbs<T>(x: &mut T, n: usize)
1237 where T: CloneableVecLike<Limb>
1238 {
1239 debug_assert!(n != 0);
1240 if !x.is_empty() {
1241 x.insert_many(0, iter::repeat(0).take(n));
1242 }
1243 }}
1244
1245 perftools_inline!{
1246 /// Shift-left buffer by n bits.
1247 pub fn ishl<T>(x: &mut T, n: usize)
1248 where T: CloneableVecLike<Limb>
1249 {
1250 let bits = <Limb as Integer>::BITS;
1251 // Need to pad with zeros for the number of `bits / <Limb as Integer>::BITS`,
1252 // and shift-left with carry for `bits % <Limb as Integer>::BITS`.
1253 let rem = n % bits;
1254 let div = n / bits;
1255 ishl_bits(x, rem);
1256 if !div.is_zero() {
1257 ishl_limbs(x, div);
1258 }
1259 }}
1260
1261 perftools_inline!{
1262 /// Shift-left buffer by n bits.
1263 pub fn shl<T>(x: &[Limb], n: usize)
1264 -> T
1265 where T: CloneableVecLike<Limb>
1266 {
1267 let mut z = T::default();
1268 z.extend_from_slice(x);
1269 ishl(&mut z, n);
1270 z
1271 }}
1272
1273 // NORMALIZE
1274
1275 perftools_inline!{
1276 /// Normalize the container by popping any leading zeros.
1277 pub fn normalize<T>(x: &mut T)
1278 where T: CloneableVecLike<Limb>
1279 {
1280 // Remove leading zero if we cause underflow. Since we're dividing
1281 // by a small power, we have at max 1 int removed.
1282 while !x.is_empty() && x.rindex(0).is_zero() {
1283 x.pop();
1284 }
1285 }}
1286
1287 // ADDITION
1288
1289 perftools_inline!{
1290 /// Implied AddAssign implementation for adding a small integer to bigint.
1291 ///
1292 /// Allows us to choose a start-index in x to store, to allow incrementing
1293 /// from a non-zero start.
1294 pub fn iadd_impl<T>(x: &mut T, y: Limb, xstart: usize)
1295 where T: CloneableVecLike<Limb>
1296 {
1297 if x.len() <= xstart {
1298 x.push(y);
1299 } else {
1300 // Initial add
1301 let mut carry = scalar::iadd(&mut x[xstart], y);
1302
1303 // Increment until overflow stops occurring.
1304 let mut size = xstart + 1;
1305 while carry && size < x.len() {
1306 carry = scalar::iadd(&mut x[size], 1);
1307 size += 1;
1308 }
1309
1310 // If we overflowed the buffer entirely, need to add 1 to the end
1311 // of the buffer.
1312 if carry {
1313 x.push(1);
1314 }
1315 }
1316 }}
1317
1318 perftools_inline!{
1319 /// AddAssign small integer to bigint.
1320 pub fn iadd<T>(x: &mut T, y: Limb)
1321 where T: CloneableVecLike<Limb>
1322 {
1323 iadd_impl(x, y, 0);
1324 }}
1325
1326 perftools_inline!{
1327 /// Add small integer to bigint.
1328 pub fn add<T>(x: &[Limb], y: Limb)
1329 -> T
1330 where T: CloneableVecLike<Limb>
1331 {
1332 let mut z = T::default();
1333 z.extend_from_slice(x);
1334 iadd(&mut z, y);
1335 z
1336 }}
1337
1338 // SUBTRACTION
1339
1340 perftools_inline!{
1341 /// SubAssign small integer to bigint.
1342 /// Does not do overflowing subtraction.
1343 pub fn isub_impl<T>(x: &mut T, y: Limb, xstart: usize)
1344 where T: CloneableVecLike<Limb>
1345 {
1346 debug_assert!(x.len() > xstart && (x[xstart] >= y || x.len() > xstart+1));
1347
1348 // Initial subtraction
1349 let mut carry = scalar::isub(&mut x[xstart], y);
1350
1351 // Increment until overflow stops occurring.
1352 let mut size = xstart + 1;
1353 while carry && size < x.len() {
1354 carry = scalar::isub(&mut x[size], 1);
1355 size += 1;
1356 }
1357 normalize(x);
1358 }}
1359
1360 perftools_inline!{
1361 /// SubAssign small integer to bigint.
1362 /// Does not do overflowing subtraction.
1363 pub fn isub<T>(x: &mut T, y: Limb)
1364 where T: CloneableVecLike<Limb>
1365 {
1366 isub_impl(x, y, 0);
1367 }}
1368
1369 perftools_inline!{
1370 /// Sub small integer to bigint.
1371 pub fn sub<T>(x: &[Limb], y: Limb)
1372 -> T
1373 where T: CloneableVecLike<Limb>
1374 {
1375 let mut z = T::default();
1376 z.extend_from_slice(x);
1377 isub(&mut z, y);
1378 z
1379 }}
1380
1381 // MULTIPLICATION
1382
1383 perftools_inline!{
1384 /// MulAssign small integer to bigint.
1385 pub fn imul<T>(x: &mut T, y: Limb)
1386 where T: CloneableVecLike<Limb>
1387 {
1388 // Multiply iteratively over all elements, adding the carry each time.
1389 let mut carry: Limb = 0;
1390 for xi in x.iter_mut() {
1391 carry = scalar::imul(xi, y, carry);
1392 }
1393
1394 // Overflow of value, add to end.
1395 if carry != 0 {
1396 x.push(carry);
1397 }
1398 }}
1399
1400 perftools_inline!{
1401 /// Mul small integer to bigint.
1402 pub fn mul<T>(x: &[Limb], y: Limb)
1403 -> T
1404 where T: CloneableVecLike<Limb>
1405 {
1406 let mut z = T::default();
1407 z.extend_from_slice(x);
1408 imul(&mut z, y);
1409 z
1410 }}
1411
1412 /// MulAssign by a power.
1413 ///
1414 /// Theoretically...
1415 ///
1416 /// Use an exponentiation by squaring method, since it reduces the time
1417 /// complexity of the multiplication to ~`O(log(n))` for the squaring,
1418 /// and `O(n*m)` for the result. Since `m` is typically a lower-order
1419 /// factor, this significantly reduces the number of multiplications
1420 /// we need to do. Iteratively multiplying by small powers follows
1421 /// the nth triangular number series, which scales as `O(p^2)`, but
1422 /// where `p` is `n+m`. In short, it scales very poorly.
1423 ///
1424 /// Practically....
1425 ///
1426 /// Exponentiation by Squaring:
1427 /// running 2 tests
1428 /// test bigcomp_f32_lexical ... bench: 1,018 ns/iter (+/- 78)
1429 /// test bigcomp_f64_lexical ... bench: 3,639 ns/iter (+/- 1,007)
1430 ///
1431 /// Exponentiation by Iterative Small Powers:
1432 /// running 2 tests
1433 /// test bigcomp_f32_lexical ... bench: 518 ns/iter (+/- 31)
1434 /// test bigcomp_f64_lexical ... bench: 583 ns/iter (+/- 47)
1435 ///
1436 /// Exponentiation by Iterative Large Powers (of 2):
1437 /// running 2 tests
1438 /// test bigcomp_f32_lexical ... bench: 671 ns/iter (+/- 31)
1439 /// test bigcomp_f64_lexical ... bench: 1,394 ns/iter (+/- 47)
1440 ///
1441 /// Even using worst-case scenarios, exponentiation by squaring is
1442 /// significantly slower for our workloads. Just multiply by small powers,
1443 /// in simple cases, and use precalculated large powers in other cases.
imul_power<T>(x: &mut T, radix: u32, n: u32) where T: CloneableVecLike<Limb>1444 pub fn imul_power<T>(x: &mut T, radix: u32, n: u32)
1445 where T: CloneableVecLike<Limb>
1446 {
1447 use super::large::KARATSUBA_CUTOFF;
1448
1449 let small_powers = get_small_powers(radix);
1450 let large_powers = get_large_powers(radix);
1451
1452 if n == 0 {
1453 // No exponent, just return.
1454 // The 0-index of the large powers is `2^0`, which is 1, so we want
1455 // to make sure we don't take that path with a literal 0.
1456 return;
1457 }
1458
1459 // We want to use the asymptotically faster algorithm if we're going
1460 // to be using Karabatsu multiplication sometime during the result,
1461 // otherwise, just use exponentiation by squaring.
1462 let bit_length = 32 - n.leading_zeros().as_usize();
1463 debug_assert!(bit_length != 0 && bit_length <= large_powers.len());
1464 if x.len() + large_powers[bit_length-1].len() < 2*KARATSUBA_CUTOFF {
1465 // We can use iterative small powers to make this faster for the
1466 // easy cases.
1467
1468 // Multiply by the largest small power until n < step.
1469 let step = small_powers.len() - 1;
1470 let power = small_powers[step];
1471 let mut n = n.as_usize();
1472 while n >= step {
1473 imul(x, power);
1474 n -= step;
1475 }
1476
1477 // Multiply by the remainder.
1478 imul(x, small_powers[n]);
1479 } else {
1480 // In theory, this code should be asymptotically a lot faster,
1481 // in practice, our small::imul seems to be the limiting step,
1482 // and large imul is slow as well.
1483
1484 // Multiply by higher order powers.
1485 let mut idx: usize = 0;
1486 let mut bit: usize = 1;
1487 let mut n = n.as_usize();
1488 while n != 0 {
1489 if n & bit != 0 {
1490 debug_assert!(idx < large_powers.len());
1491 large::imul(x, large_powers[idx]);
1492 n ^= bit;
1493 }
1494 idx += 1;
1495 bit <<= 1;
1496 }
1497 }
1498 }
1499
1500 perftools_inline!{
1501 /// Mul by a power.
1502 pub fn mul_power<T>(x: &[Limb], radix: u32, n: u32)
1503 -> T
1504 where T: CloneableVecLike<Limb>
1505 {
1506 let mut z = T::default();
1507 z.extend_from_slice(x);
1508 imul_power(&mut z, radix, n);
1509 z
1510 }}
1511
1512 // DIVISION
1513
1514 perftools_inline!{
1515 /// DivAssign small integer to bigint and get the remainder.
1516 pub fn idiv<T>(x: &mut T, y: Limb)
1517 -> Limb
1518 where T: CloneableVecLike<Limb>
1519 {
1520 // Divide iteratively over all elements, adding the carry each time.
1521 let mut rem: Limb = 0;
1522 for xi in x.iter_mut().rev() {
1523 rem = scalar::idiv(xi, y, rem);
1524 }
1525 normalize(x);
1526
1527 rem
1528 }}
1529
1530 perftools_inline!{
1531 /// Div small integer to bigint and get the remainder.
1532 pub fn div<T>(x: &[Limb], y: Limb)
1533 -> (T, Limb)
1534 where T: CloneableVecLike<Limb>
1535 {
1536 let mut z = T::default();
1537 z.extend_from_slice(x);
1538 let rem = idiv(&mut z, y);
1539 (z, rem)
1540 }}
1541
1542 // POWER
1543
1544 perftools_inline!{
1545 /// Calculate x^n, using exponentiation by squaring.
1546 ///
1547 /// This algorithm is slow, using `mul_power` should generally be preferred,
1548 /// as although it's not asymptotically faster, it precalculates a lot
1549 /// of results.
1550 pub fn ipow<T>(x: &mut T, mut n: Limb)
1551 where T: CloneableVecLike<Limb>
1552 {
1553 // Store `x` as 1, and switch `base` to `x`.
1554 let mut base = T::default();
1555 base.push(1);
1556 mem::swap(x, &mut base);
1557
1558 // Do main algorithm.
1559 loop {
1560 if n.is_odd() {
1561 large::imul(x, &base);
1562 }
1563 n /= 2;
1564
1565 // We need to break as a post-condition, since the real work
1566 // is in the `imul` and `mul` algorithms.
1567 if n.is_zero() {
1568 break;
1569 } else {
1570 base = large::mul(&base, &base);
1571 }
1572 }
1573 }}
1574
1575 perftools_inline!{
1576 /// Calculate x^n, using exponentiation by squaring.
1577 pub fn pow<T>(x: &[Limb], n: Limb)
1578 -> T
1579 where T: CloneableVecLike<Limb>
1580 {
1581 let mut z = T::default();
1582 z.extend_from_slice(x);
1583 ipow(&mut z, n);
1584 z
1585 }}
1586
1587 } // small
1588
1589 // LARGE
1590 // -----
1591
1592 // Large-to-large operations, to modify a big integer from a native scalar.
1593
1594 pub(in crate::atof::algorithm) mod large {
1595
1596 use crate::lib::cmp;
1597 use super::*;
1598
1599 // RELATIVE OPERATORS
1600
1601 perftools_inline!{
1602 /// Compare `x` to `y`, in little-endian order.
1603 pub fn compare(x: &[Limb], y: &[Limb]) -> cmp::Ordering {
1604 if x.len() > y.len() {
1605 return cmp::Ordering::Greater;
1606 } else if x.len() < y.len() {
1607 return cmp::Ordering::Less;
1608 } else {
1609 let iter = x.iter().rev().zip(y.iter().rev());
1610 for (&xi, &yi) in iter {
1611 if xi > yi {
1612 return cmp::Ordering::Greater;
1613 } else if xi < yi {
1614 return cmp::Ordering::Less;
1615 }
1616 }
1617 // Equal case.
1618 return cmp::Ordering::Equal;
1619 }
1620 }}
1621
1622 perftools_inline!{
1623 /// Check if x is greater than y.
1624 pub fn greater(x: &[Limb], y: &[Limb]) -> bool {
1625 compare(x, y) == cmp::Ordering::Greater
1626 }}
1627
1628 perftools_inline!{
1629 /// Check if x is greater than or equal to y.
1630 pub fn greater_equal(x: &[Limb], y: &[Limb]) -> bool {
1631 !less(x, y)
1632 }}
1633
1634 perftools_inline!{
1635 /// Check if x is less than y.
1636 pub fn less(x: &[Limb], y: &[Limb]) -> bool {
1637 compare(x, y) == cmp::Ordering::Less
1638 }}
1639
1640 perftools_inline!{
1641 /// Check if x is less than or equal to y.
1642 pub fn less_equal(x: &[Limb], y: &[Limb]) -> bool {
1643 !greater(x, y)
1644 }}
1645
1646 perftools_inline!{
1647 /// Check if x is equal to y.
1648 /// Slightly optimized for equality comparisons, since it reduces the number
1649 /// of comparisons relative to `compare`.
1650 pub fn equal(x: &[Limb], y: &[Limb]) -> bool {
1651 let mut iter = x.iter().rev().zip(y.iter().rev());
1652 x.len() == y.len() && iter.all(|(&xi, &yi)| xi == yi)
1653 }}
1654
1655 /// ADDITION
1656
1657 /// Implied AddAssign implementation for bigints.
1658 ///
1659 /// Allows us to choose a start-index in x to store, so we can avoid
1660 /// padding the buffer with zeros when not needed, optimized for vectors.
iadd_impl<T>(x: &mut T, y: &[Limb], xstart: usize) where T: CloneableVecLike<Limb>1661 pub fn iadd_impl<T>(x: &mut T, y: &[Limb], xstart: usize)
1662 where T: CloneableVecLike<Limb>
1663 {
1664 // The effective x buffer is from `xstart..x.len()`, so we need to treat
1665 // that as the current range. If the effective y buffer is longer, need
1666 // to resize to that, + the start index.
1667 if y.len() > x.len() - xstart {
1668 x.resize(y.len() + xstart, 0);
1669 }
1670
1671 // Iteratively add elements from y to x.
1672 let mut carry = false;
1673 for (xi, yi) in (&mut x[xstart..]).iter_mut().zip(y.iter()) {
1674 // Only one op of the two can overflow, since we added at max
1675 // Limb::max_value() + Limb::max_value(). Add the previous carry,
1676 // and store the current carry for the next.
1677 let mut tmp = scalar::iadd(xi, *yi);
1678 if carry {
1679 tmp |= scalar::iadd(xi, 1);
1680 }
1681 carry = tmp;
1682 }
1683
1684 // Overflow from the previous bit.
1685 if carry {
1686 small::iadd_impl(x, 1, y.len()+xstart);
1687 }
1688 }
1689
1690 perftools_inline!{
1691 /// AddAssign bigint to bigint.
1692 pub fn iadd<T>(x: &mut T, y: &[Limb])
1693 where T: CloneableVecLike<Limb>
1694 {
1695 iadd_impl(x, y, 0)
1696 }}
1697
1698 perftools_inline!{
1699 /// Add bigint to bigint.
1700 pub fn add<T>(x: &[Limb], y: &[Limb])
1701 -> T
1702 where T: CloneableVecLike<Limb>
1703 {
1704 let mut z = T::default();
1705 z.extend_from_slice(x);
1706 iadd(&mut z, y);
1707 z
1708 }}
1709
1710 // SUBTRACTION
1711
1712 /// SubAssign bigint to bigint.
isub<T>(x: &mut T, y: &[Limb]) where T: CloneableVecLike<Limb>1713 pub fn isub<T>(x: &mut T, y: &[Limb])
1714 where T: CloneableVecLike<Limb>
1715 {
1716 // Basic underflow checks.
1717 debug_assert!(greater_equal(x, y));
1718
1719 // Iteratively add elements from y to x.
1720 let mut carry = false;
1721 for (xi, yi) in x.iter_mut().zip(y.iter()) {
1722 // Only one op of the two can overflow, since we added at max
1723 // Limb::max_value() + Limb::max_value(). Add the previous carry,
1724 // and store the current carry for the next.
1725 let mut tmp = scalar::isub(xi, *yi);
1726 if carry {
1727 tmp |= scalar::isub(xi, 1);
1728 }
1729 carry = tmp;
1730 }
1731
1732 if carry {
1733 small::isub_impl(x, 1, y.len());
1734 } else {
1735 small::normalize(x);
1736 }
1737 }
1738
1739 perftools_inline!{
1740 /// Sub bigint to bigint.
1741 pub fn sub<T>(x: &[Limb], y: &[Limb])
1742 -> T
1743 where T: CloneableVecLike<Limb>
1744 {
1745 let mut z = T::default();
1746 z.extend_from_slice(x);
1747 isub(&mut z, y);
1748 z
1749 }}
1750
1751 // MULTIPLICATIION
1752
1753 /// Number of digits to bottom-out to asymptotically slow algorithms.
1754 ///
1755 /// Karatsuba tends to out-perform long-multiplication at ~320-640 bits,
1756 /// so we go halfway, while Newton division tends to out-perform
1757 /// Algorithm D at ~1024 bits. We can toggle this for optimal performance.
1758 pub const KARATSUBA_CUTOFF: usize = 32;
1759
1760 /// Grade-school multiplication algorithm.
1761 ///
1762 /// Slow, naive algorithm, using limb-bit bases and just shifting left for
1763 /// each iteration. This could be optimized with numerous other algorithms,
1764 /// but it's extremely simple, and works in O(n*m) time, which is fine
1765 /// by me. Each iteration, of which there are `m` iterations, requires
1766 /// `n` multiplications, and `n` additions, or grade-school multiplication.
long_mul<T>(x: &[Limb], y: &[Limb]) -> T where T: CloneableVecLike<Limb>1767 fn long_mul<T>(x: &[Limb], y: &[Limb])
1768 -> T
1769 where T: CloneableVecLike<Limb>
1770 {
1771 // Using the immutable value, multiply by all the scalars in y, using
1772 // the algorithm defined above. Use a single buffer to avoid
1773 // frequent reallocations. Handle the first case to avoid a redundant
1774 // addition, since we know y.len() >= 1.
1775 let mut z: T = small::mul(x, y[0]);
1776 z.resize(x.len() + y.len(), 0);
1777
1778 // Handle the iterative cases.
1779 for (i, &yi) in y[1..].iter().enumerate() {
1780 let zi: T = small::mul(x, yi);
1781 iadd_impl(&mut z, &zi, i+1);
1782 }
1783
1784 small::normalize(&mut z);
1785
1786 z
1787 }
1788
1789 perftools_inline!{
1790 /// Split two buffers into halfway, into (lo, hi).
1791 pub fn karatsuba_split<'a>(z: &'a [Limb], m: usize)
1792 -> (&'a [Limb], &'a [Limb])
1793 {
1794 (&z[..m], &z[m..])
1795 }}
1796
1797 /// Karatsuba multiplication algorithm with roughly equal input sizes.
1798 ///
1799 /// Assumes `y.len() >= x.len()`.
karatsuba_mul<T>(x: &[Limb], y: &[Limb]) -> T where T: CloneableVecLike<Limb>1800 fn karatsuba_mul<T>(x: &[Limb], y: &[Limb])
1801 -> T
1802 where T: CloneableVecLike<Limb>
1803 {
1804 if y.len() <= KARATSUBA_CUTOFF {
1805 // Bottom-out to long division for small cases.
1806 long_mul(x, y)
1807 } else if x.len() < y.len() / 2 {
1808 karatsuba_uneven_mul(x, y)
1809 } else {
1810 // Do our 3 multiplications.
1811 let m = y.len() / 2;
1812 let (xl, xh) = karatsuba_split(x, m);
1813 let (yl, yh) = karatsuba_split(y, m);
1814 let sumx: T = add(xl, xh);
1815 let sumy: T = add(yl, yh);
1816 let z0: T = karatsuba_mul(xl, yl);
1817 let mut z1: T = karatsuba_mul(&sumx, &sumy);
1818 let z2: T = karatsuba_mul(xh, yh);
1819 // Properly scale z1, which is `z1 - z2 - zo`.
1820 isub(&mut z1, &z2);
1821 isub(&mut z1, &z0);
1822
1823 // Create our result, which is equal to, in little-endian order:
1824 // [z0, z1 - z2 - z0, z2]
1825 // z1 must be shifted m digits (2^(32m)) over.
1826 // z2 must be shifted 2*m digits (2^(64m)) over.
1827 let mut result = T::default();
1828 let len = z0.len().max(m + z1.len()).max(2*m + z2.len());
1829 result.reserve_exact(len);
1830 result.extend_from_slice(&z0);
1831 iadd_impl(&mut result, &z1, m);
1832 iadd_impl(&mut result, &z2, 2*m);
1833
1834 result
1835 }
1836 }
1837
1838 /// Karatsuba multiplication algorithm where y is substantially larger than x.
1839 ///
1840 /// Assumes `y.len() >= x.len()`.
karatsuba_uneven_mul<T>(x: &[Limb], mut y: &[Limb]) -> T where T: CloneableVecLike<Limb>1841 fn karatsuba_uneven_mul<T>(x: &[Limb], mut y: &[Limb])
1842 -> T
1843 where T: CloneableVecLike<Limb>
1844 {
1845 let mut result = T::default();
1846 result.resize(x.len() + y.len(), 0);
1847
1848 // This effectively is like grade-school multiplication between
1849 // two numbers, except we're using splits on `y`, and the intermediate
1850 // step is a Karatsuba multiplication.
1851 let mut start = 0;
1852 while y.len() != 0 {
1853 let m = x.len().min(y.len());
1854 let (yl, yh) = karatsuba_split(y, m);
1855 let prod: T = karatsuba_mul(x, yl);
1856 iadd_impl(&mut result, &prod, start);
1857 y = yh;
1858 start += m;
1859 }
1860 small::normalize(&mut result);
1861
1862 result
1863 }
1864
1865 perftools_inline!{
1866 /// Forwarder to the proper Karatsuba algorithm.
1867 fn karatsuba_mul_fwd<T>(x: &[Limb], y: &[Limb])
1868 -> T
1869 where T: CloneableVecLike<Limb>
1870 {
1871 if x.len() < y.len() {
1872 karatsuba_mul(x, y)
1873 } else {
1874 karatsuba_mul(y, x)
1875 }
1876 }}
1877
1878 perftools_inline!{
1879 /// MulAssign bigint to bigint.
1880 pub fn imul<T>(x: &mut T, y: &[Limb])
1881 where T: CloneableVecLike<Limb>
1882 {
1883 if y.len() == 1 {
1884 small::imul(x, y[0]);
1885 } else {
1886 // We're not really in a condition where using Karatsuba
1887 // multiplication makes sense, so we're just going to use long
1888 // division. ~20% speedup compared to:
1889 // *x = karatsuba_mul_fwd(x, y);
1890 *x = karatsuba_mul_fwd(x, y);
1891 }
1892 }}
1893
1894 perftools_inline!{
1895 /// Mul bigint to bigint.
1896 pub fn mul<T>(x: &[Limb], y: &[Limb])
1897 -> T
1898 where T: CloneableVecLike<Limb>
1899 {
1900 let mut z = T::default();
1901 z.extend_from_slice(x);
1902 imul(&mut z, y);
1903 z
1904 }}
1905
1906 // DIVISION
1907
1908 /// Constants for algorithm D.
1909 const ALGORITHM_D_B: Wide = 1 << <Limb as Integer>::BITS;
1910 const ALGORITHM_D_M: Wide = ALGORITHM_D_B - 1;
1911
1912 /// Calculate qhat (an estimate for the quotient).
1913 ///
1914 /// This is step D3 in Algorithm D in "The Art of Computer Programming".
1915 /// Assumes `x.len() > y.len()` and `y.len() >= 2`.
1916 ///
1917 /// * `j` - Current index on the iteration of the loop.
calculate_qhat(x: &[Limb], y: &[Limb], j: usize) -> Wide1918 fn calculate_qhat(x: &[Limb], y: &[Limb], j: usize)
1919 -> Wide
1920 {
1921 let n = y.len();
1922
1923 // Estimate qhat of q[j]
1924 // Original Code:
1925 // qhat = (x[j+n]*B + x[j+n-1])/y[n-1];
1926 // rhat = (x[j+n]*B + x[j+n-1]) - qhat*y[n-1];
1927 let x_jn = as_wide(x[j+n]);
1928 let x_jn1 = as_wide(x[j+n-1]);
1929 let num = (x_jn << <Limb as Integer>::BITS) + x_jn1;
1930 let den = as_wide(y[n-1]);
1931 let mut qhat = num / den;
1932 let mut rhat = num - qhat * den;
1933
1934 // Scale qhat and rhat
1935 // Original Code:
1936 // again:
1937 // if (qhat >= B || qhat*y[n-2] > B*rhat + x[j+n-2])
1938 // { qhat = qhat - 1;
1939 // rhat = rhat + y[n-1];
1940 // if (rhat < B) goto again;
1941 // }
1942 let x_jn2 = as_wide(x[j+n-2]);
1943 let y_n2 = as_wide(y[n-2]);
1944 let y_n1 = as_wide(y[n-1]);
1945 // This only happens when the leading bit of qhat is set.
1946 while qhat >= ALGORITHM_D_B || qhat * y_n2 > (rhat << <Limb as Integer>::BITS) + x_jn2 {
1947 qhat -= 1;
1948 rhat += y_n1;
1949 if rhat >= ALGORITHM_D_B {
1950 break;
1951 }
1952 }
1953
1954 qhat
1955 }
1956
1957 /// Multiply and subtract.
1958 ///
1959 /// This is step D4 in Algorithm D in "The Art of Computer Programming",
1960 /// and returns the remainder.
multiply_and_subtract<T>(x: &mut T, y: &T, qhat: Wide, j: usize) -> SignedWide where T: CloneableVecLike<Limb>1961 fn multiply_and_subtract<T>(x: &mut T, y: &T, qhat: Wide, j: usize)
1962 -> SignedWide
1963 where T: CloneableVecLike<Limb>
1964 {
1965 let n = y.len();
1966
1967 // Multiply and subtract
1968 // Original Code:
1969 // k = 0;
1970 // for (i = 0; i < n; i++) {
1971 // p = qhat*y[i];
1972 // t = x[i+j] - k - (p & 0xFFFFFFFFLL);
1973 // x[i+j] = t;
1974 // k = (p >> 32) - (t >> 32);
1975 // }
1976 // t = x[j+n] - k;
1977 // x[j+n] = t;
1978 let mut k: SignedWide = 0;
1979 let mut t: SignedWide;
1980 for i in 0..n {
1981 let x_ij = as_signed_wide(x[i+j]);
1982 let y_i = as_wide(y[i]);
1983 let p = qhat * y_i;
1984 t = x_ij.wrapping_sub(k).wrapping_sub(as_signed_wide(p & ALGORITHM_D_M));
1985 x[i+j] = as_limb(t);
1986 k = as_signed_wide(p >> <Limb as Integer>::BITS) - (t >> <Limb as Integer>::BITS);
1987 }
1988 t = as_signed_wide(x[j+n]) - k;
1989 x[j+n] = as_limb(t);
1990
1991 t
1992 }
1993
1994 perftools_inline!{
1995 /// Calculate the quotient from the estimate and the test.
1996 ///
1997 /// This is a mix of step D5 and D6 in Algorithm D, so the algorithm
1998 /// may work for single passes, without a quotient buffer.
1999 fn test_quotient(qhat: Wide, t: SignedWide)
2000 -> Wide
2001 {
2002 if t < 0 {
2003 qhat - 1
2004 } else {
2005 qhat
2006 }
2007 }}
2008
2009 /// Add back.
2010 ///
2011 /// This is step D6 in Algorithm D in "The Art of Computer Programming",
2012 /// and adds back the remainder on the very unlikely scenario we overestimated
2013 /// the quotient by 1. Subtract 1 from the quotient, and add back the
2014 /// remainder.
2015 ///
2016 /// This step should be specifically debugged, due to its low likelihood,
2017 /// since the probability is ~2/b, where b in this case is 2^32 or 2^64.
add_back<T>(x: &mut T, y: &T, mut t: SignedWide, j: usize) where T: CloneableVecLike<Limb>2018 fn add_back<T>(x: &mut T, y: &T, mut t: SignedWide, j: usize)
2019 where T: CloneableVecLike<Limb>
2020 {
2021 let n = y.len();
2022
2023 // Store quotient digits
2024 // If we subtracted too much, add back.
2025 // Original Code:
2026 // q[j] = qhat; // Store quotient digit.
2027 // if (t < 0) { // If we subtracted too
2028 // q[j] = q[j] - 1; // much, add back.
2029 // k = 0;
2030 // for (i = 0; i < n; i++) {
2031 // t = (unsigned long long)x[i+j] + y[i] + k;
2032 // x[i+j] = t;
2033 // k = t >> 32;
2034 // }
2035 // x[j+n] = x[j+n] + k;
2036 // }
2037 if t < 0 {
2038 let mut k: SignedWide = 0;
2039 for i in 0..n {
2040 t = as_signed_wide(as_wide(x[i+j]) + as_wide(y[i])) + k;
2041 x[i+j] = as_limb(t);
2042 k = t >> <Limb as Integer>::BITS;
2043 }
2044 let x_jn = as_signed_wide(x[j+n]) + k;
2045 x[j+n] = as_limb(x_jn);
2046 }
2047 }
2048
2049 /// Calculate the remainder from the quotient.
2050 ///
2051 /// This is step D8 in Algorithm D in "The Art of Computer Programming",
2052 /// and "unnormalizes" to calculate the remainder from the quotient.
calculate_remainder<T>(x: &[Limb], y: &[Limb], s: usize) -> T where T: CloneableVecLike<Limb>2053 fn calculate_remainder<T>(x: &[Limb], y: &[Limb], s: usize)
2054 -> T
2055 where T: CloneableVecLike<Limb>
2056 {
2057 // Calculate the remainder.
2058 // Original Code:
2059 // for (i = 0; i < n-1; i++)
2060 // r[i] = (x[i] >> s) | ((unsigned long long)x[i+1] << (32-s));
2061 // r[n-1] = x[n-1] >> s;
2062 let n = y.len();
2063 let mut r = T::default();
2064 r.reserve_exact(n);
2065 let rs = <Limb as Integer>::BITS - s;
2066 for i in 0..n-1 {
2067 let xi = as_wide(x[i]) >> s;
2068 let xi1 = as_wide(x[i+1]) << rs;
2069 let ri = xi | xi1;
2070 r.push(as_limb(ri));
2071 }
2072 let x_n1 = x[n-1] >> s;
2073 r.push(as_limb(x_n1));
2074
2075 r
2076 }
2077
2078 /// Implementation of Knuth's Algorithm D, and return the quotient and remainder.
2079 ///
2080 /// `x` is the dividend, and `y` is the divisor.
2081 /// Assumes `x.len() > y.len()` and `y.len() >= 2`.
2082 ///
2083 /// Based off the Hacker's Delight implementation of Knuth's Algorithm D
2084 /// in "The Art of Computer Programming".
2085 /// http://www.hackersdelight.org/hdcodetxt/divmnu64.c.txt
2086 ///
2087 /// All Hacker's Delight code is public domain, so this routine shall
2088 /// also be placed in the public domain. See:
2089 /// https://www.hackersdelight.org/permissions.htm
algorithm_d_div<T>(x: &[Limb], y: &[Limb]) -> (T, T) where T: CloneableVecLike<Limb>2090 fn algorithm_d_div<T>(x: &[Limb], y: &[Limb])
2091 -> (T, T)
2092 where T: CloneableVecLike<Limb>
2093 {
2094 // Normalize the divisor so the leading-bit is set to 1.
2095 // x is the dividend, y is the divisor.
2096 // Need a leading zero on the numerator.
2097 let s = y.rindex(0).leading_zeros().as_usize();
2098 let m = x.len();
2099 let n = y.len();
2100 let mut xn: T = small::shl_bits(x, s);
2101 let yn: T = small::shl_bits(y, s);
2102 xn.push(0);
2103
2104 // Store certain variables for the algorithm.
2105 let mut q = T::default();
2106 q.resize(m-n+1, 0);
2107 for j in (0..m-n+1).rev() {
2108 // Estimate the quotient
2109 let mut qhat = calculate_qhat(&xn, &yn, j);
2110 if qhat != 0 {
2111 let t = multiply_and_subtract(&mut xn, &yn, qhat, j);
2112 qhat = test_quotient(qhat, t);
2113 add_back(&mut xn, &yn, t, j);
2114 }
2115 q[j] = as_limb(qhat);
2116 }
2117 let mut r = calculate_remainder(&xn, &yn, s);
2118
2119 // Normalize our results
2120 small::normalize(&mut q);
2121 small::normalize(&mut r);
2122
2123 (q, r)
2124 }
2125
2126 perftools_inline!{
2127 /// DivAssign bigint to bigint.
2128 pub fn idiv<T>(x: &mut T, y: &[Limb])
2129 -> T
2130 where T: CloneableVecLike<Limb>
2131 {
2132 debug_assert!(y.len() != 0);
2133
2134 if x.len() < y.len() {
2135 // Can optimize easily, since the quotient is 0,
2136 // and the remainder is x. Put before `y.len() == 1`, since
2137 // it optimizes when `x.len() == 0` nicely.
2138 let mut r = T::default();
2139 mem::swap(x, &mut r);
2140 r
2141 } else if y.len() == 1 {
2142 // Can optimize for division by a small value.
2143 let mut r = T::default();
2144 r.push(small::idiv(x, y[0]));
2145 r
2146 } else {
2147 let (q, r) = algorithm_d_div(x, y);
2148 *x = q;
2149 r
2150 }
2151 }}
2152
2153 perftools_inline!{
2154 /// Div bigint to bigint.
2155 pub fn div<T>(x: &[Limb], y: &[Limb])
2156 -> (T, T)
2157 where T: CloneableVecLike<Limb>
2158 {
2159 let mut z = T::default();
2160 z.extend_from_slice(x);
2161 let rem = idiv(&mut z, y);
2162 (z, rem)
2163 }}
2164
2165 /// Emit a single digit for the quotient and store the remainder in-place.
2166 ///
2167 /// An extremely efficient division algorithm for small quotients, requiring
2168 /// you to know the full range of the quotient prior to use. For example,
2169 /// with a quotient that can range from [0, 10), you must have 4 leading
2170 /// zeros in the divisor, so we can use a single-limb division to get
2171 /// an accurate estimate of the quotient. Since we always underestimate
2172 /// the quotient, we can add 1 and then emit the digit.
2173 ///
2174 /// Requires a non-normalized denominator, with at least [1-6] leading
2175 /// zeros, depending on the base (for example, 1 for base2, 6 for base36).
2176 ///
2177 /// Adapted from David M. Gay's dtoa, and therefore under an MIT license:
2178 /// www.netlib.org/fp/dtoa.c
quorem<T>(x: &mut T, y: &T) -> Limb where T: CloneableVecLike<Limb>2179 pub fn quorem<T>(x: &mut T, y: &T)
2180 -> Limb
2181 where T: CloneableVecLike<Limb>
2182 {
2183 debug_assert!(y.len() > 0);
2184 let mask = as_wide(Limb::max_value());
2185
2186 // Numerator is smaller the denominator, quotient always 0.
2187 let m = x.len();
2188 let n = y.len();
2189 if m < n {
2190 return 0;
2191 }
2192
2193 // Calculate our initial estimate for q
2194 let mut q = x[m-1] / (y[n-1] + 1);
2195
2196 // Need to calculate the remainder if we don't have a 0 quotient.
2197 if q != 0 {
2198 let mut borrow: Wide = 0;
2199 let mut carry: Wide = 0;
2200 for j in 0..m {
2201 let p = as_wide(y[j]) * as_wide(q) + carry;
2202 carry = p >> <Limb as Integer>::BITS;
2203 let t = as_wide(x[j]).wrapping_sub(p & mask).wrapping_sub(borrow);
2204 borrow = (t >> <Limb as Integer>::BITS) & 1;
2205 x[j] = as_limb(t);
2206 }
2207 small::normalize(x);
2208 }
2209
2210 // Check if we under-estimated x.
2211 if greater_equal(x, y) {
2212 q += 1;
2213 let mut borrow: Wide = 0;
2214 let mut carry: Wide = 0;
2215 for j in 0..m {
2216 let p = as_wide(y[j]) + carry;
2217 carry = p >> <Limb as Integer>::BITS;
2218 let t = as_wide(x[j]).wrapping_sub(p & mask).wrapping_sub(borrow);
2219 borrow = (t >> <Limb as Integer>::BITS) & 1;
2220 x[j] = as_limb(t);
2221 }
2222 small::normalize(x);
2223 }
2224
2225 q
2226 }
2227
2228 } // large
2229
2230 use crate::lib::cmp;
2231 use super::small_powers::*;
2232 use super::large_powers::*;
2233
2234 /// Generate the imul_pown wrappers.
2235 macro_rules! imul_power {
2236 ($name:ident, $base:expr) => (
2237 perftools_inline!{
2238 /// Multiply by a power of $base.
2239 fn $name(&mut self, n: u32) {
2240 self.imul_power_impl($base, n)
2241 }}
2242 );
2243 }
2244
2245 // TRAITS
2246 // ------
2247
2248 /// Traits for shared operations for big integers.
2249 ///
2250 /// None of these are implemented using normal traits, since these
2251 /// are very expensive operations, and we want to deliberately
2252 /// and explicitly use these functions.
2253 pub(in crate::atof::algorithm) trait SharedOps: Clone + Sized + Default {
2254 /// Underlying storage type for a SmallOps.
2255 type StorageType: CloneableVecLike<Limb>;
2256
2257 // DATA
2258
2259 /// Get access to the underlying data
data<'a>(&'a self) -> &'a Self::StorageType2260 fn data<'a>(&'a self) -> &'a Self::StorageType;
2261
2262 /// Get access to the underlying data
data_mut<'a>(&'a mut self) -> &'a mut Self::StorageType2263 fn data_mut<'a>(&'a mut self) -> &'a mut Self::StorageType;
2264
2265 // ZERO
2266
2267 perftools_inline!{
2268 /// Check if the value is a normalized 0.
2269 fn is_zero(&self) -> bool {
2270 self.limb_length() == 0
2271 }}
2272
2273 // RELATIVE OPERATIONS
2274
2275 perftools_inline!{
2276 /// Compare self to y.
2277 fn compare(&self, y: &Self) -> cmp::Ordering {
2278 large::compare(self.data(), y.data())
2279 }}
2280
2281 perftools_inline!{
2282 /// Check if self is greater than y.
2283 fn greater(&self, y: &Self) -> bool {
2284 large::greater(self.data(), y.data())
2285 }}
2286
2287 perftools_inline!{
2288 /// Check if self is greater than or equal to y.
2289 fn greater_equal(&self, y: &Self) -> bool {
2290 large::greater_equal(self.data(), y.data())
2291 }}
2292
2293 perftools_inline!{
2294 /// Check if self is less than y.
2295 fn less(&self, y: &Self) -> bool {
2296 large::less(self.data(), y.data())
2297 }}
2298
2299 perftools_inline!{
2300 /// Check if self is less than or equal to y.
2301 fn less_equal(&self, y: &Self) -> bool {
2302 large::less_equal(self.data(), y.data())
2303 }}
2304
2305 perftools_inline!{
2306 /// Check if self is equal to y.
2307 fn equal(&self, y: &Self) -> bool {
2308 large::equal(self.data(), y.data())
2309 }}
2310
2311 // PROPERTIES
2312
2313 perftools_inline!{
2314 /// Get the number of leading zero digits in the storage.
2315 /// Assumes the value is normalized.
2316 fn leading_zero_limbs(&self) -> usize {
2317 small::leading_zero_limbs(self.data())
2318 }}
2319
2320 perftools_inline!{
2321 /// Get the number of trailing zero digits in the storage.
2322 /// Assumes the value is normalized.
2323 fn trailing_zero_limbs(&self) -> usize {
2324 small::trailing_zero_limbs(self.data())
2325 }}
2326
2327 perftools_inline!{
2328 /// Get number of leading zero bits in the storage.
2329 /// Assumes the value is normalized.
2330 fn leading_zeros(&self) -> usize {
2331 small::leading_zeros(self.data())
2332 }}
2333
2334 perftools_inline!{
2335 /// Get number of trailing zero bits in the storage.
2336 /// Assumes the value is normalized.
2337 fn trailing_zeros(&self) -> usize {
2338 small::trailing_zeros(self.data())
2339 }}
2340
2341 perftools_inline!{
2342 /// Calculate the bit-length of the big-integer.
2343 /// Returns usize::max_value() if the value overflows,
2344 /// IE, if `self.data().len() > usize::max_value() / 8`.
2345 fn bit_length(&self) -> usize {
2346 small::bit_length(self.data())
2347 }}
2348
2349 perftools_inline!{
2350 /// Calculate the digit-length of the big-integer.
2351 fn limb_length(&self) -> usize {
2352 small::limb_length(self.data())
2353 }}
2354
2355 perftools_inline!{
2356 /// Get the high 16-bits from the bigint and if there are remaining bits.
2357 fn hi16(&self) -> (u16, bool) {
2358 self.data().as_slice().hi16()
2359 }}
2360
2361 perftools_inline!{
2362 /// Get the high 32-bits from the bigint and if there are remaining bits.
2363 fn hi32(&self) -> (u32, bool) {
2364 self.data().as_slice().hi32()
2365 }}
2366
2367 perftools_inline!{
2368 /// Get the high 64-bits from the bigint and if there are remaining bits.
2369 fn hi64(&self) -> (u64, bool) {
2370 self.data().as_slice().hi64()
2371 }}
2372
2373 perftools_inline!{
2374 /// Get the high 128-bits from the bigint and if there are remaining bits.
2375 fn hi128(&self) -> (u128, bool) {
2376 self.data().as_slice().hi128()
2377 }}
2378
2379 perftools_inline!{
2380 /// Pad the buffer with zeros to the least-significant bits.
2381 fn pad_zero_digits(&mut self, n: usize) -> usize {
2382 small::ishl_limbs(self.data_mut(), n);
2383 n
2384 }}
2385
2386 // INTEGER CONVERSIONS
2387
2388 // CREATION
2389
2390 perftools_inline!{
2391 /// Create new big integer from u16.
2392 fn from_u16(x: u16) -> Self {
2393 let mut v = Self::default();
2394 let slc = split_u16(x);
2395 v.data_mut().extend_from_slice(&slc);
2396 v.normalize();
2397 v
2398 }}
2399
2400 perftools_inline!{
2401 /// Create new big integer from u32.
2402 fn from_u32(x: u32) -> Self {
2403 let mut v = Self::default();
2404 let slc = split_u32(x);
2405 v.data_mut().extend_from_slice(&slc);
2406 v.normalize();
2407 v
2408 }}
2409
2410 perftools_inline!{
2411 /// Create new big integer from u64.
2412 fn from_u64(x: u64) -> Self {
2413 let mut v = Self::default();
2414 let slc = split_u64(x);
2415 v.data_mut().extend_from_slice(&slc);
2416 v.normalize();
2417 v
2418 }}
2419
2420 perftools_inline!{
2421 /// Create new big integer from u128.
2422 fn from_u128(x: u128) -> Self {
2423 let mut v = Self::default();
2424 let slc = split_u128(x);
2425 v.data_mut().extend_from_slice(&slc);
2426 v.normalize();
2427 v
2428 }}
2429
2430 // NORMALIZE
2431
2432 perftools_inline!{
2433 /// Normalize the integer, so any leading zero values are removed.
2434 fn normalize(&mut self) {
2435 small::normalize(self.data_mut());
2436 }}
2437
2438 perftools_inline!{
2439 /// Get if the big integer is normalized.
2440 fn is_normalized(&self) -> bool {
2441 self.data().is_empty() || !self.data().rindex(0).is_zero()
2442 }}
2443
2444 // SHIFTS
2445
2446 perftools_inline!{
2447 /// Shift-left the entire buffer n bits.
2448 fn ishl(&mut self, n: usize) {
2449 small::ishl(self.data_mut(), n);
2450 }}
2451
2452 perftools_inline!{
2453 /// Shift-left the entire buffer n bits.
2454 fn shl(&self, n: usize) -> Self {
2455 let mut x = self.clone();
2456 x.ishl(n);
2457 x
2458 }}
2459
2460 perftools_inline!{
2461 /// Shift-right the entire buffer n bits.
2462 fn ishr(&mut self, n: usize, mut roundup: bool) {
2463 roundup &= small::ishr(self.data_mut(), n);
2464
2465 // Round-up the least significant bit.
2466 if roundup {
2467 if self.data().is_empty() {
2468 self.data_mut().push(1);
2469 } else {
2470 self.data_mut()[0] += 1;
2471 }
2472 }
2473 }}
2474
2475 perftools_inline!{
2476 /// Shift-right the entire buffer n bits.
2477 fn shr(&self, n: usize, roundup: bool) -> Self {
2478 let mut x = self.clone();
2479 x.ishr(n, roundup);
2480 x
2481 }}
2482 }
2483
2484 /// Trait for small operations for arbitrary-precision numbers.
2485 pub(in crate::atof::algorithm) trait SmallOps: SharedOps {
2486 // SMALL POWERS
2487
2488 perftools_inline!{
2489 /// Get the small powers from the radix.
2490 fn small_powers(radix: u32) -> &'static [Limb] {
2491 get_small_powers(radix)
2492 }}
2493
2494 perftools_inline!{
2495 /// Get the large powers from the radix.
2496 fn large_powers(radix: u32) -> &'static [&'static [Limb]] {
2497 get_large_powers(radix)
2498 }}
2499
2500 // ADDITION
2501
2502 perftools_inline!{
2503 /// AddAssign small integer.
2504 fn iadd_small(&mut self, y: Limb) {
2505 small::iadd(self.data_mut(), y);
2506 }}
2507
2508 perftools_inline!{
2509 /// Add small integer to a copy of self.
2510 fn add_small(&self, y: Limb) -> Self {
2511 let mut x = self.clone();
2512 x.iadd_small(y);
2513 x
2514 }}
2515
2516 // SUBTRACTION
2517
2518 perftools_inline!{
2519 /// SubAssign small integer.
2520 /// Warning: Does no overflow checking, x must be >= y.
2521 fn isub_small(&mut self, y: Limb) {
2522 small::isub(self.data_mut(), y);
2523 }}
2524
2525 perftools_inline!{
2526 /// Sub small integer to a copy of self.
2527 /// Warning: Does no overflow checking, x must be >= y.
2528 fn sub_small(&mut self, y: Limb) -> Self {
2529 let mut x = self.clone();
2530 x.isub_small(y);
2531 x
2532 }}
2533
2534 // MULTIPLICATION
2535
2536 perftools_inline!{
2537 /// MulAssign small integer.
2538 fn imul_small(&mut self, y: Limb) {
2539 small::imul(self.data_mut(), y);
2540 }}
2541
2542 perftools_inline!{
2543 /// Mul small integer to a copy of self.
2544 fn mul_small(&self, y: Limb) -> Self {
2545 let mut x = self.clone();
2546 x.imul_small(y);
2547 x
2548 }}
2549
2550 perftools_inline!{
2551 /// MulAssign by a power.
2552 fn imul_power_impl(&mut self, radix: u32, n: u32) {
2553 small::imul_power(self.data_mut(), radix, n);
2554 }}
2555
2556 perftools_inline!{
2557 fn imul_power(&mut self, radix: u32, n: u32) {
2558 match radix {
2559 2 => self.imul_pow2(n),
2560 3 => self.imul_pow3(n),
2561 4 => self.imul_pow4(n),
2562 5 => self.imul_pow5(n),
2563 6 => self.imul_pow6(n),
2564 7 => self.imul_pow7(n),
2565 8 => self.imul_pow8(n),
2566 9 => self.imul_pow9(n),
2567 10 => self.imul_pow10(n),
2568 11 => self.imul_pow11(n),
2569 12 => self.imul_pow12(n),
2570 13 => self.imul_pow13(n),
2571 14 => self.imul_pow14(n),
2572 15 => self.imul_pow15(n),
2573 16 => self.imul_pow16(n),
2574 17 => self.imul_pow17(n),
2575 18 => self.imul_pow18(n),
2576 19 => self.imul_pow19(n),
2577 20 => self.imul_pow20(n),
2578 21 => self.imul_pow21(n),
2579 22 => self.imul_pow22(n),
2580 23 => self.imul_pow23(n),
2581 24 => self.imul_pow24(n),
2582 25 => self.imul_pow25(n),
2583 26 => self.imul_pow26(n),
2584 27 => self.imul_pow27(n),
2585 28 => self.imul_pow28(n),
2586 29 => self.imul_pow29(n),
2587 30 => self.imul_pow30(n),
2588 31 => self.imul_pow31(n),
2589 32 => self.imul_pow32(n),
2590 33 => self.imul_pow33(n),
2591 34 => self.imul_pow34(n),
2592 35 => self.imul_pow35(n),
2593 36 => self.imul_pow36(n),
2594 _ => unreachable!()
2595 }
2596 }}
2597
2598 perftools_inline!{
2599 /// Multiply by a power of 2.
2600 fn imul_pow2(&mut self, n: u32) {
2601 self.ishl(n.as_usize())
2602 }}
2603
2604 imul_power!(imul_pow3, 3);
2605
2606 perftools_inline!{
2607 /// Multiply by a power of 4.
2608 fn imul_pow4(&mut self, n: u32) {
2609 self.imul_pow2(2*n);
2610 }}
2611
2612 imul_power!(imul_pow5, 5);
2613
2614 perftools_inline!{
2615 /// Multiply by a power of 6.
2616 fn imul_pow6(&mut self, n: u32) {
2617 self.imul_pow3(n);
2618 self.imul_pow2(n);
2619 }}
2620
2621 imul_power!(imul_pow7, 7);
2622
2623 perftools_inline!{
2624 /// Multiply by a power of 8.
2625 fn imul_pow8(&mut self, n: u32) {
2626 self.imul_pow2(3*n);
2627 }}
2628
2629 perftools_inline!{
2630 /// Multiply by a power of 9.
2631 fn imul_pow9(&mut self, n: u32) {
2632 self.imul_pow3(n);
2633 self.imul_pow3(n);
2634 }}
2635
2636 perftools_inline!{
2637 /// Multiply by a power of 10.
2638 fn imul_pow10(&mut self, n: u32) {
2639 self.imul_pow5(n);
2640 self.imul_pow2(n);
2641 }}
2642
2643 imul_power!(imul_pow11, 11);
2644
2645 perftools_inline!{
2646 /// Multiply by a power of 12.
2647 fn imul_pow12(&mut self, n: u32) {
2648 self.imul_pow3(n);
2649 self.imul_pow4(n);
2650 }}
2651
2652 imul_power!(imul_pow13, 13);
2653
2654 perftools_inline!{
2655 /// Multiply by a power of 14.
2656 fn imul_pow14(&mut self, n: u32) {
2657 self.imul_pow7(n);
2658 self.imul_pow2(n);
2659 }}
2660
2661 perftools_inline!{
2662 /// Multiply by a power of 15.
2663 fn imul_pow15(&mut self, n: u32) {
2664 self.imul_pow3(n);
2665 self.imul_pow5(n);
2666 }}
2667
2668 perftools_inline!{
2669 /// Multiply by a power of 16.
2670 fn imul_pow16(&mut self, n: u32) {
2671 self.imul_pow2(4*n);
2672 }}
2673
2674 imul_power!(imul_pow17, 17);
2675
2676 perftools_inline!{
2677 /// Multiply by a power of 18.
2678 fn imul_pow18(&mut self, n: u32) {
2679 self.imul_pow9(n);
2680 self.imul_pow2(n);
2681 }}
2682
2683 imul_power!(imul_pow19, 19);
2684
2685 perftools_inline!{
2686 /// Multiply by a power of 20.
2687 fn imul_pow20(&mut self, n: u32) {
2688 self.imul_pow5(n);
2689 self.imul_pow4(n);
2690 }}
2691
2692 perftools_inline!{
2693 /// Multiply by a power of 21.
2694 fn imul_pow21(&mut self, n: u32) {
2695 self.imul_pow3(n);
2696 self.imul_pow7(n);
2697 }}
2698
2699 perftools_inline!{
2700 /// Multiply by a power of 22.
2701 fn imul_pow22(&mut self, n: u32) {
2702 self.imul_pow11(n);
2703 self.imul_pow2(n);
2704 }}
2705
2706 imul_power!(imul_pow23, 23);
2707
2708 perftools_inline!{
2709 /// Multiply by a power of 24.
2710 fn imul_pow24(&mut self, n: u32) {
2711 self.imul_pow3(n);
2712 self.imul_pow8(n);
2713 }}
2714
2715 perftools_inline!{
2716 /// Multiply by a power of 25.
2717 fn imul_pow25(&mut self, n: u32) {
2718 self.imul_pow5(n);
2719 self.imul_pow5(n);
2720 }}
2721
2722 perftools_inline!{
2723 /// Multiply by a power of 26.
2724 fn imul_pow26(&mut self, n: u32) {
2725 self.imul_pow13(n);
2726 self.imul_pow2(n);
2727 }}
2728
2729 perftools_inline!{
2730 /// Multiply by a power of 27.
2731 fn imul_pow27(&mut self, n: u32) {
2732 self.imul_pow9(n);
2733 self.imul_pow3(n);
2734 }}
2735
2736 perftools_inline!{
2737 /// Multiply by a power of 28.
2738 fn imul_pow28(&mut self, n: u32) {
2739 self.imul_pow7(n);
2740 self.imul_pow4(n);
2741 }}
2742
2743 imul_power!(imul_pow29, 29);
2744
2745 perftools_inline!{
2746 /// Multiply by a power of 30.
2747 fn imul_pow30(&mut self, n: u32) {
2748 self.imul_pow15(n);
2749 self.imul_pow2(n);
2750 }}
2751
2752 imul_power!(imul_pow31, 31);
2753
2754 perftools_inline!{
2755 /// Multiply by a power of 32.
2756 fn imul_pow32(&mut self, n: u32) {
2757 self.imul_pow2(5*n);
2758 }}
2759
2760 perftools_inline!{
2761 /// Multiply by a power of 33.
2762 fn imul_pow33(&mut self, n: u32) {
2763 self.imul_pow3(n);
2764 self.imul_pow11(n);
2765 }}
2766
2767 perftools_inline!{
2768 /// Multiply by a power of 34.
2769 fn imul_pow34(&mut self, n: u32) {
2770 self.imul_pow17(n);
2771 self.imul_pow2(n);
2772 }}
2773
2774 perftools_inline!{
2775 /// Multiply by a power of 35.
2776 fn imul_pow35(&mut self, n: u32) {
2777 self.imul_pow5(n);
2778 self.imul_pow7(n);
2779 }}
2780
2781 perftools_inline!{
2782 /// Multiply by a power of 36.
2783 fn imul_pow36(&mut self, n: u32) {
2784 self.imul_pow9(n);
2785 self.imul_pow4(n);
2786 }}
2787
2788 // DIVISION
2789
2790 perftools_inline!{
2791 /// DivAssign small integer, and return the remainder.
2792 fn idiv_small(&mut self, y: Limb) -> Limb {
2793 small::idiv(self.data_mut(), y)
2794 }}
2795
2796 perftools_inline!{
2797 /// Div small integer to a copy of self, and return the remainder.
2798 fn div_small(&self, y: Limb) -> (Self, Limb) {
2799 let mut x = self.clone();
2800 let rem = x.idiv_small(y);
2801 (x, rem)
2802 }}
2803
2804 // POWER
2805
2806 perftools_inline!{
2807 /// Calculate self^n
2808 fn ipow(&mut self, n: Limb) {
2809 small::ipow(self.data_mut(), n);
2810 }}
2811
2812 perftools_inline!{
2813 /// Calculate self^n
2814 fn pow(&self, n: Limb) -> Self {
2815 let mut x = self.clone();
2816 x.ipow(n);
2817 x
2818 }}
2819 }
2820
2821 /// Trait for large operations for arbitrary-precision numbers.
2822 pub(in crate::atof::algorithm) trait LargeOps: SmallOps {
2823 // ADDITION
2824
2825 perftools_inline!{
2826 /// AddAssign large integer.
2827 fn iadd_large(&mut self, y: &Self) {
2828 large::iadd(self.data_mut(), y.data());
2829 }}
2830
2831 perftools_inline!{
2832 /// Add large integer to a copy of self.
2833 fn add_large(&mut self, y: &Self) -> Self {
2834 let mut x = self.clone();
2835 x.iadd_large(y);
2836 x
2837 }}
2838
2839 // SUBTRACTION
2840
2841 perftools_inline!{
2842 /// SubAssign large integer.
2843 /// Warning: Does no overflow checking, x must be >= y.
2844 fn isub_large(&mut self, y: &Self) {
2845 large::isub(self.data_mut(), y.data());
2846 }}
2847
2848 perftools_inline!{
2849 /// Sub large integer to a copy of self.
2850 /// Warning: Does no overflow checking, x must be >= y.
2851 fn sub_large(&mut self, y: &Self) -> Self {
2852 let mut x = self.clone();
2853 x.isub_large(y);
2854 x
2855 }}
2856
2857 // MULTIPLICATION
2858
2859 perftools_inline!{
2860 /// MulAssign large integer.
2861 fn imul_large(&mut self, y: &Self) {
2862 large::imul(self.data_mut(), y.data());
2863 }}
2864
2865 perftools_inline!{
2866 /// Mul large integer to a copy of self.
2867 fn mul_large(&mut self, y: &Self) -> Self {
2868 let mut x = self.clone();
2869 x.imul_large(y);
2870 x
2871 }}
2872
2873 // DIVISION
2874
2875 perftools_inline!{
2876 /// DivAssign large integer and get remainder.
2877 fn idiv_large(&mut self, y: &Self) -> Self {
2878 let mut rem = Self::default();
2879 *rem.data_mut() = large::idiv(self.data_mut(), y.data());
2880 rem
2881 }}
2882
2883 perftools_inline!{
2884 /// Div large integer to a copy of self and get quotient and remainder.
2885 fn div_large(&mut self, y: &Self) -> (Self, Self) {
2886 let mut x = self.clone();
2887 let rem = x.idiv_large(y);
2888 (x, rem)
2889 }}
2890
2891 perftools_inline!{
2892 /// Calculate the fast quotient for a single limb-bit quotient.
2893 ///
2894 /// This requires a non-normalized divisor, where there at least
2895 /// `integral_binary_factor` 0 bits set, to ensure at maximum a single
2896 /// digit will be produced for a single base.
2897 ///
2898 /// Warning: This is not a general-purpose division algorithm,
2899 /// it is highly specialized for peeling off singular digits.
2900 fn quorem(&mut self, y: &Self) -> Limb {
2901 large::quorem(self.data_mut(), y.data())
2902 }}
2903 }
2904
2905 #[cfg(test)]
2906 mod tests {
2907 use crate::util::test::*;
2908 use super::*;
2909
2910 #[derive(Clone, Default)]
2911 struct Bigint {
2912 data: DataType,
2913 }
2914
2915 impl Bigint {
2916 #[inline]
new() -> Bigint2917 pub fn new() -> Bigint {
2918 Bigint { data: arrvec![] }
2919 }
2920 }
2921
2922 impl SharedOps for Bigint {
2923 type StorageType = DataType;
2924
2925 #[inline]
data<'a>(&'a self) -> &'a Self::StorageType2926 fn data<'a>(&'a self) -> &'a Self::StorageType {
2927 &self.data
2928 }
2929
2930 #[inline]
data_mut<'a>(&'a mut self) -> &'a mut Self::StorageType2931 fn data_mut<'a>(&'a mut self) -> &'a mut Self::StorageType {
2932 &mut self.data
2933 }
2934 }
2935
2936 impl SmallOps for Bigint {
2937 }
2938
2939 impl LargeOps for Bigint {
2940 }
2941
2942 // SHARED OPS
2943
2944 #[test]
greater_test()2945 fn greater_test() {
2946 // Simple
2947 let x = Bigint { data: from_u32(&[1]) };
2948 let y = Bigint { data: from_u32(&[2]) };
2949 assert!(!x.greater(&y));
2950 assert!(!x.greater(&x));
2951 assert!(y.greater(&x));
2952
2953 // Check asymmetric
2954 let x = Bigint { data: from_u32(&[5, 1]) };
2955 let y = Bigint { data: from_u32(&[2]) };
2956 assert!(x.greater(&y));
2957 assert!(!x.greater(&x));
2958 assert!(!y.greater(&x));
2959
2960 // Check when we use reverse ordering properly.
2961 let x = Bigint { data: from_u32(&[5, 1, 9]) };
2962 let y = Bigint { data: from_u32(&[6, 2, 8]) };
2963 assert!(x.greater(&y));
2964 assert!(!x.greater(&x));
2965 assert!(!y.greater(&x));
2966
2967 // Complex scenario, check it properly uses reverse ordering.
2968 let x = Bigint { data: from_u32(&[0, 1, 9]) };
2969 let y = Bigint { data: from_u32(&[4294967295, 0, 9]) };
2970 assert!(x.greater(&y));
2971 assert!(!x.greater(&x));
2972 assert!(!y.greater(&x));
2973 }
2974
2975 #[test]
greater_equal_test()2976 fn greater_equal_test() {
2977 // Simple
2978 let x = Bigint { data: from_u32(&[1]) };
2979 let y = Bigint { data: from_u32(&[2]) };
2980 assert!(!x.greater_equal(&y));
2981 assert!(x.greater_equal(&x));
2982 assert!(y.greater_equal(&x));
2983
2984 // Check asymmetric
2985 let x = Bigint { data: from_u32(&[5, 1]) };
2986 let y = Bigint { data: from_u32(&[2]) };
2987 assert!(x.greater_equal(&y));
2988 assert!(x.greater_equal(&x));
2989 assert!(!y.greater_equal(&x));
2990
2991 // Check when we use reverse ordering properly.
2992 let x = Bigint { data: from_u32(&[5, 1, 9]) };
2993 let y = Bigint { data: from_u32(&[6, 2, 8]) };
2994 assert!(x.greater_equal(&y));
2995 assert!(x.greater_equal(&x));
2996 assert!(!y.greater_equal(&x));
2997
2998 // Complex scenario, check it properly uses reverse ordering.
2999 let x = Bigint { data: from_u32(&[0, 1, 9]) };
3000 let y = Bigint { data: from_u32(&[4294967295, 0, 9]) };
3001 assert!(x.greater_equal(&y));
3002 assert!(x.greater_equal(&x));
3003 assert!(!y.greater_equal(&x));
3004 }
3005
3006 #[test]
equal_test()3007 fn equal_test() {
3008 // Simple
3009 let x = Bigint { data: from_u32(&[1]) };
3010 let y = Bigint { data: from_u32(&[2]) };
3011 assert!(!x.equal(&y));
3012 assert!(x.equal(&x));
3013 assert!(!y.equal(&x));
3014
3015 // Check asymmetric
3016 let x = Bigint { data: from_u32(&[5, 1]) };
3017 let y = Bigint { data: from_u32(&[2]) };
3018 assert!(!x.equal(&y));
3019 assert!(x.equal(&x));
3020 assert!(!y.equal(&x));
3021
3022 // Check when we use reverse ordering properly.
3023 let x = Bigint { data: from_u32(&[5, 1, 9]) };
3024 let y = Bigint { data: from_u32(&[6, 2, 8]) };
3025 assert!(!x.equal(&y));
3026 assert!(x.equal(&x));
3027 assert!(!y.equal(&x));
3028
3029 // Complex scenario, check it properly uses reverse ordering.
3030 let x = Bigint { data: from_u32(&[0, 1, 9]) };
3031 let y = Bigint { data: from_u32(&[4294967295, 0, 9]) };
3032 assert!(!x.equal(&y));
3033 assert!(x.equal(&x));
3034 assert!(!y.equal(&x));
3035 }
3036
3037 #[test]
less_test()3038 fn less_test() {
3039 // Simple
3040 let x = Bigint { data: from_u32(&[1]) };
3041 let y = Bigint { data: from_u32(&[2]) };
3042 assert!(x.less(&y));
3043 assert!(!x.less(&x));
3044 assert!(!y.less(&x));
3045
3046 // Check asymmetric
3047 let x = Bigint { data: from_u32(&[5, 1]) };
3048 let y = Bigint { data: from_u32(&[2]) };
3049 assert!(!x.less(&y));
3050 assert!(!x.less(&x));
3051 assert!(y.less(&x));
3052
3053 // Check when we use reverse ordering properly.
3054 let x = Bigint { data: from_u32(&[5, 1, 9]) };
3055 let y = Bigint { data: from_u32(&[6, 2, 8]) };
3056 assert!(!x.less(&y));
3057 assert!(!x.less(&x));
3058 assert!(y.less(&x));
3059
3060 // Complex scenario, check it properly uses reverse ordering.
3061 let x = Bigint { data: from_u32(&[0, 1, 9]) };
3062 let y = Bigint { data: from_u32(&[4294967295, 0, 9]) };
3063 assert!(!x.less(&y));
3064 assert!(!x.less(&x));
3065 assert!(y.less(&x));
3066 }
3067
3068 #[test]
less_equal_test()3069 fn less_equal_test() {
3070 // Simple
3071 let x = Bigint { data: from_u32(&[1]) };
3072 let y = Bigint { data: from_u32(&[2]) };
3073 assert!(x.less_equal(&y));
3074 assert!(x.less_equal(&x));
3075 assert!(!y.less_equal(&x));
3076
3077 // Check asymmetric
3078 let x = Bigint { data: from_u32(&[5, 1]) };
3079 let y = Bigint { data: from_u32(&[2]) };
3080 assert!(!x.less_equal(&y));
3081 assert!(x.less_equal(&x));
3082 assert!(y.less_equal(&x));
3083
3084 // Check when we use reverse ordering properly.
3085 let x = Bigint { data: from_u32(&[5, 1, 9]) };
3086 let y = Bigint { data: from_u32(&[6, 2, 8]) };
3087 assert!(!x.less_equal(&y));
3088 assert!(x.less_equal(&x));
3089 assert!(y.less_equal(&x));
3090
3091 // Complex scenario, check it properly uses reverse ordering.
3092 let x = Bigint { data: from_u32(&[0, 1, 9]) };
3093 let y = Bigint { data: from_u32(&[4294967295, 0, 9]) };
3094 assert!(!x.less_equal(&y));
3095 assert!(x.less_equal(&x));
3096 assert!(y.less_equal(&x));
3097 }
3098
3099 #[test]
leading_zero_limbs_test()3100 fn leading_zero_limbs_test() {
3101 assert_eq!(Bigint::new().leading_zero_limbs(), 0);
3102
3103 assert_eq!(Bigint::from_u16(0xF).leading_zero_limbs(), 0);
3104 assert_eq!(Bigint::from_u32(0xFF).leading_zero_limbs(), 0);
3105 assert_eq!(Bigint::from_u64(0xFF00000000).leading_zero_limbs(), 0);
3106 assert_eq!(Bigint::from_u128(0xFF000000000000000000000000).leading_zero_limbs(), 0);
3107
3108 assert_eq!(Bigint::from_u16(0xF).leading_zero_limbs(), 0);
3109 assert_eq!(Bigint::from_u32(0xF).leading_zero_limbs(), 0);
3110 assert_eq!(Bigint::from_u64(0xF00000000).leading_zero_limbs(), 0);
3111 assert_eq!(Bigint::from_u128(0xF000000000000000000000000).leading_zero_limbs(), 0);
3112
3113 assert_eq!(Bigint::from_u16(0xF0).leading_zero_limbs(), 0);
3114 assert_eq!(Bigint::from_u32(0xF0).leading_zero_limbs(), 0);
3115 assert_eq!(Bigint::from_u64(0xF000000000).leading_zero_limbs(), 0);
3116 assert_eq!(Bigint::from_u128(0xF0000000000000000000000000).leading_zero_limbs(), 0);
3117 }
3118
3119 #[test]
trailing_zero_limbs_test()3120 fn trailing_zero_limbs_test() {
3121 assert_eq!(Bigint::new().trailing_zero_limbs(), 0);
3122
3123 assert_eq!(Bigint { data: arrvec![0xFF] }.trailing_zero_limbs(), 0);
3124 assert_eq!(Bigint { data: arrvec![0, 0xFF000] }.trailing_zero_limbs(), 1);
3125 assert_eq!(Bigint { data: arrvec![0, 0, 0, 0xFF000] }.trailing_zero_limbs(), 3);
3126 }
3127
3128 #[test]
leading_zeros_test()3129 fn leading_zeros_test() {
3130 assert_eq!(Bigint::new().leading_zeros(), 0);
3131
3132 assert_eq!(Bigint::from_u16(0xFF).leading_zeros(), <Limb as Integer>::BITS-8);
3133 assert_eq!(Bigint::from_u32(0xFF).leading_zeros(), <Limb as Integer>::BITS-8);
3134 assert_eq!(Bigint::from_u64(0xFF00000000).leading_zeros(), 24);
3135 assert_eq!(Bigint::from_u128(0xFF000000000000000000000000).leading_zeros(), 24);
3136
3137 assert_eq!(Bigint::from_u16(0xF).leading_zeros(), <Limb as Integer>::BITS-4);
3138 assert_eq!(Bigint::from_u32(0xF).leading_zeros(), <Limb as Integer>::BITS-4);
3139 assert_eq!(Bigint::from_u64(0xF00000000).leading_zeros(), 28);
3140 assert_eq!(Bigint::from_u128(0xF000000000000000000000000).leading_zeros(), 28);
3141
3142 assert_eq!(Bigint::from_u16(0xF0).leading_zeros(), <Limb as Integer>::BITS-8);
3143 assert_eq!(Bigint::from_u32(0xF0).leading_zeros(), <Limb as Integer>::BITS-8);
3144 assert_eq!(Bigint::from_u64(0xF000000000).leading_zeros(), 24);
3145 assert_eq!(Bigint::from_u128(0xF0000000000000000000000000).leading_zeros(), 24);
3146 }
3147
3148 #[test]
trailing_zeros_test()3149 fn trailing_zeros_test() {
3150 assert_eq!(Bigint::new().trailing_zeros(), 0);
3151
3152 assert_eq!(Bigint::from_u16(0xFF).trailing_zeros(), 0);
3153 assert_eq!(Bigint::from_u32(0xFF).trailing_zeros(), 0);
3154 assert_eq!(Bigint::from_u64(0xFF00000000).trailing_zeros(), 32);
3155 assert_eq!(Bigint::from_u128(0xFF000000000000000000000000).trailing_zeros(), 96);
3156
3157 assert_eq!(Bigint::from_u16(0xF).trailing_zeros(), 0);
3158 assert_eq!(Bigint::from_u32(0xF).trailing_zeros(), 0);
3159 assert_eq!(Bigint::from_u64(0xF00000000).trailing_zeros(), 32);
3160 assert_eq!(Bigint::from_u128(0xF000000000000000000000000).trailing_zeros(), 96);
3161
3162 assert_eq!(Bigint::from_u16(0xF0).trailing_zeros(), 4);
3163 assert_eq!(Bigint::from_u32(0xF0).trailing_zeros(), 4);
3164 assert_eq!(Bigint::from_u64(0xF000000000).trailing_zeros(), 36);
3165 assert_eq!(Bigint::from_u128(0xF0000000000000000000000000).trailing_zeros(), 100);
3166 }
3167
3168 #[test]
hi32_test()3169 fn hi32_test() {
3170 assert_eq!(Bigint::from_u16(0xA).hi32(), (0xA0000000, false));
3171 assert_eq!(Bigint::from_u32(0xAB).hi32(), (0xAB000000, false));
3172 assert_eq!(Bigint::from_u64(0xAB00000000).hi32(), (0xAB000000, false));
3173 assert_eq!(Bigint::from_u64(0xA23456789A).hi32(), (0xA2345678, true));
3174 }
3175
3176 #[test]
hi64_test()3177 fn hi64_test() {
3178 assert_eq!(Bigint::from_u16(0xA).hi64(), (0xA000000000000000, false));
3179 assert_eq!(Bigint::from_u32(0xAB).hi64(), (0xAB00000000000000, false));
3180 assert_eq!(Bigint::from_u64(0xAB00000000).hi64(), (0xAB00000000000000, false));
3181 assert_eq!(Bigint::from_u64(0xA23456789A).hi64(), (0xA23456789A000000, false));
3182 assert_eq!(Bigint::from_u128(0xABCDEF0123456789ABCDEF0123).hi64(), (0xABCDEF0123456789, true));
3183 }
3184
3185 #[test]
hi128_test()3186 fn hi128_test() {
3187 assert_eq!(Bigint::from_u128(0xABCDEF0123456789ABCDEF0123).hi128(), (0xABCDEF0123456789ABCDEF0123000000, false));
3188 assert_eq!(Bigint::from_u128(0xABCDEF0123456789ABCDEF0123456789).hi128(), (0xABCDEF0123456789ABCDEF0123456789, false));
3189 assert_eq!(Bigint { data: from_u32(&[0x34567890, 0xBCDEF012, 0x3456789A, 0xBCDEF012, 0xA]) }.hi128(), (0xABCDEF0123456789ABCDEF0123456789, false));
3190 assert_eq!(Bigint { data: from_u32(&[0x34567891, 0xBCDEF012, 0x3456789A, 0xBCDEF012, 0xA]) }.hi128(), (0xABCDEF0123456789ABCDEF0123456789, true));
3191 }
3192
3193 #[test]
pad_zero_digits_test()3194 fn pad_zero_digits_test() {
3195 let mut x = Bigint { data: arrvec![0, 0, 0, 1] };
3196 x.pad_zero_digits(3);
3197 assert_eq!(x.data.as_slice(), &[0, 0, 0, 0, 0, 0, 1]);
3198
3199 let mut x = Bigint { data: arrvec![1] };
3200 x.pad_zero_digits(1);
3201 assert_eq!(x.data.as_slice(), &[0, 1]);
3202 }
3203
3204 #[test]
shl_test()3205 fn shl_test() {
3206 // Pattern generated via `''.join(["1" +"0"*i for i in range(20)])`
3207 let mut big = Bigint { data: from_u32(&[0xD2210408]) };
3208 big.ishl(5);
3209 assert_eq!(big.data, from_u32(&[0x44208100, 0x1A]));
3210 big.ishl(32);
3211 assert_eq!(big.data, from_u32(&[0, 0x44208100, 0x1A]));
3212 big.ishl(27);
3213 assert_eq!(big.data, from_u32(&[0, 0, 0xD2210408]));
3214
3215 // 96-bits of previous pattern
3216 let mut big = Bigint { data: from_u32(&[0x20020010, 0x8040100, 0xD2210408]) };
3217 big.ishl(5);
3218 assert_eq!(big.data, from_u32(&[0x400200, 0x802004, 0x44208101, 0x1A]));
3219 big.ishl(32);
3220 assert_eq!(big.data, from_u32(&[0, 0x400200, 0x802004, 0x44208101, 0x1A]));
3221 big.ishl(27);
3222 assert_eq!(big.data, from_u32(&[0, 0, 0x20020010, 0x8040100, 0xD2210408]));
3223 }
3224
3225 #[test]
shr_test()3226 fn shr_test() {
3227 // Simple case.
3228 let mut big = Bigint { data: from_u32(&[0xD2210408]) };
3229 big.ishr(5, false);
3230 assert_eq!(big.data, from_u32(&[0x6910820]));
3231 big.ishr(27, false);
3232 assert_eq!(big.data, from_u32(&[]));
3233
3234 // Pattern generated via `''.join(["1" +"0"*i for i in range(20)])`
3235 let mut big = Bigint { data: from_u32(&[0x20020010, 0x8040100, 0xD2210408]) };
3236 big.ishr(5, false);
3237 assert_eq!(big.data, from_u32(&[0x1001000, 0x40402008, 0x6910820]));
3238 big.ishr(32, false);
3239 assert_eq!(big.data, from_u32(&[0x40402008, 0x6910820]));
3240 big.ishr(27, false);
3241 assert_eq!(big.data, from_u32(&[0xD2210408]));
3242
3243 // Check no-roundup with halfway and even
3244 let mut big = Bigint { data: from_u32(&[0xD2210408]) };
3245 big.ishr(3, true);
3246 assert_eq!(big.data, from_u32(&[0x1A442081]));
3247 big.ishr(1, true);
3248 assert_eq!(big.data, from_u32(&[0xD221040]));
3249
3250 let mut big = Bigint { data: from_u32(&[0xD2210408]) };
3251 big.ishr(4, true);
3252 assert_eq!(big.data, from_u32(&[0xD221040]));
3253
3254 // Check roundup with halfway and odd
3255 let mut big = Bigint { data: from_u32(&[0xD2210438]) };
3256 big.ishr(3, true);
3257 assert_eq!(big.data, from_u32(&[0x1A442087]));
3258 big.ishr(1, true);
3259 assert_eq!(big.data, from_u32(&[0xD221044]));
3260
3261 let mut big = Bigint { data: from_u32(&[0xD2210438]) };
3262 big.ishr(5, true);
3263 assert_eq!(big.data, from_u32(&[0x6910822]));
3264 }
3265
3266 #[test]
bit_length_test()3267 fn bit_length_test() {
3268 let x = Bigint { data: from_u32(&[0, 0, 0, 1]) };
3269 assert_eq!(x.bit_length(), 97);
3270
3271 let x = Bigint { data: from_u32(&[0, 0, 0, 3]) };
3272 assert_eq!(x.bit_length(), 98);
3273
3274 let x = Bigint { data: from_u32(&[1<<31]) };
3275 assert_eq!(x.bit_length(), 32);
3276 }
3277
3278 // SMALL OPS
3279
3280 #[test]
iadd_small_test()3281 fn iadd_small_test() {
3282 // Overflow check (single)
3283 // This should set all the internal data values to 0, the top
3284 // value to (1<<31), and the bottom value to (4>>1).
3285 // This is because the max_value + 1 leads to all 0s, we set the
3286 // topmost bit to 1.
3287 let mut x = Bigint { data: from_u32(&[4294967295]) };
3288 x.iadd_small(5);
3289 assert_eq!(x.data, from_u32(&[4, 1]));
3290
3291 // No overflow, single value
3292 let mut x = Bigint { data: from_u32(&[5]) };
3293 x.iadd_small(7);
3294 assert_eq!(x.data, from_u32(&[12]));
3295
3296 // Single carry, internal overflow
3297 let mut x = Bigint::from_u64(0x80000000FFFFFFFF);
3298 x.iadd_small(7);
3299 assert_eq!(x.data, from_u32(&[6, 0x80000001]));
3300
3301 // Double carry, overflow
3302 let mut x = Bigint::from_u64(0xFFFFFFFFFFFFFFFF);
3303 x.iadd_small(7);
3304 assert_eq!(x.data, from_u32(&[6, 0, 1]));
3305 }
3306
3307 #[test]
isub_small_test()3308 fn isub_small_test() {
3309 // Overflow check (single)
3310 let mut x = Bigint { data: from_u32(&[4, 1]) };
3311 x.isub_small(5);
3312 assert_eq!(x.data, from_u32(&[4294967295]));
3313
3314 // No overflow, single value
3315 let mut x = Bigint { data: from_u32(&[12]) };
3316 x.isub_small(7);
3317 assert_eq!(x.data, from_u32(&[5]));
3318
3319 // Single carry, internal overflow
3320 let mut x = Bigint { data: from_u32(&[6, 0x80000001]) };
3321 x.isub_small(7);
3322 assert_eq!(x.data, from_u32(&[0xFFFFFFFF, 0x80000000]));
3323
3324 // Double carry, overflow
3325 let mut x = Bigint { data: from_u32(&[6, 0, 1]) };
3326 x.isub_small(7);
3327 assert_eq!(x.data, from_u32(&[0xFFFFFFFF, 0xFFFFFFFF]));
3328 }
3329
3330 #[test]
imul_small_test()3331 fn imul_small_test() {
3332 // No overflow check, 1-int.
3333 let mut x = Bigint { data: from_u32(&[5]) };
3334 x.imul_small(7);
3335 assert_eq!(x.data, from_u32(&[35]));
3336
3337 // No overflow check, 2-ints.
3338 let mut x = Bigint::from_u64(0x4000000040000);
3339 x.imul_small(5);
3340 assert_eq!(x.data, from_u32(&[0x00140000, 0x140000]));
3341
3342 // Overflow, 1 carry.
3343 let mut x = Bigint { data: from_u32(&[0x33333334]) };
3344 x.imul_small(5);
3345 assert_eq!(x.data, from_u32(&[4, 1]));
3346
3347 // Overflow, 1 carry, internal.
3348 let mut x = Bigint::from_u64(0x133333334);
3349 x.imul_small(5);
3350 assert_eq!(x.data, from_u32(&[4, 6]));
3351
3352 // Overflow, 2 carries.
3353 let mut x = Bigint::from_u64(0x3333333333333334);
3354 x.imul_small(5);
3355 assert_eq!(x.data, from_u32(&[4, 0, 1]));
3356 }
3357
3358 #[test]
idiv_small_test()3359 fn idiv_small_test() {
3360 let mut x = Bigint { data: from_u32(&[4]) };
3361 assert_eq!(x.idiv_small(7), 4);
3362 assert_eq!(x.data, from_u32(&[]));
3363
3364 let mut x = Bigint { data: from_u32(&[3]) };
3365 assert_eq!(x.idiv_small(7), 3);
3366 assert_eq!(x.data, from_u32(&[]));
3367
3368 // Check roundup, odd, halfway
3369 let mut x = Bigint { data: from_u32(&[15]) };
3370 assert_eq!(x.idiv_small(10), 5);
3371 assert_eq!(x.data, from_u32(&[1]));
3372
3373 // Check 1 carry.
3374 let mut x = Bigint::from_u64(0x133333334);
3375 assert_eq!(x.idiv_small(5), 1);
3376 assert_eq!(x.data, from_u32(&[0x3D70A3D7]));
3377
3378 // Check 2 carries.
3379 let mut x = Bigint::from_u64(0x3333333333333334);
3380 assert_eq!(x.idiv_small(5), 4);
3381 assert_eq!(x.data, from_u32(&[0xD70A3D70, 0xA3D70A3]));
3382 }
3383
3384 #[test]
ipow_test()3385 fn ipow_test() {
3386 let x = Bigint { data: from_u32(&[5]) };
3387 assert_eq!(x.pow(2).data, from_u32(&[25]));
3388 assert_eq!(x.pow(15).data, from_u32(&[452807053, 7]));
3389 assert_eq!(x.pow(16).data, from_u32(&[2264035265, 35]));
3390 assert_eq!(x.pow(17).data, from_u32(&[2730241733, 177]));
3391 assert_eq!(x.pow(302).data, from_u32(&[2443090281, 2149694430, 2297493928, 1584384001, 1279504719, 1930002239, 3312868939, 3735173465, 3523274756, 2025818732, 1641675015, 2431239749, 4292780461, 3719612855, 4174476133, 3296847770, 2677357556, 638848153, 2198928114, 3285049351, 2159526706, 626302612]));
3392 }
3393
3394 // LARGE OPS
3395
3396 #[test]
iadd_large_test()3397 fn iadd_large_test() {
3398 // Overflow, both single values
3399 let mut x = Bigint { data: from_u32(&[4294967295]) };
3400 let y = Bigint { data: from_u32(&[5]) };
3401 x.iadd_large(&y);
3402 assert_eq!(x.data, from_u32(&[4, 1]));
3403
3404 // No overflow, single value
3405 let mut x = Bigint { data: from_u32(&[5]) };
3406 let y = Bigint { data: from_u32(&[7]) };
3407 x.iadd_large(&y);
3408 assert_eq!(x.data, from_u32(&[12]));
3409
3410 // Single carry, internal overflow
3411 let mut x = Bigint::from_u64(0x80000000FFFFFFFF);
3412 let y = Bigint { data: from_u32(&[7]) };
3413 x.iadd_large(&y);
3414 assert_eq!(x.data, from_u32(&[6, 0x80000001]));
3415
3416 // 1st overflows, 2nd doesn't.
3417 let mut x = Bigint::from_u64(0x7FFFFFFFFFFFFFFF);
3418 let y = Bigint::from_u64(0x7FFFFFFFFFFFFFFF);
3419 x.iadd_large(&y);
3420 assert_eq!(x.data, from_u32(&[0xFFFFFFFE, 0xFFFFFFFF]));
3421
3422 // Both overflow.
3423 let mut x = Bigint::from_u64(0x8FFFFFFFFFFFFFFF);
3424 let y = Bigint::from_u64(0x7FFFFFFFFFFFFFFF);
3425 x.iadd_large(&y);
3426 assert_eq!(x.data, from_u32(&[0xFFFFFFFE, 0x0FFFFFFF, 1]));
3427 }
3428
3429 #[test]
isub_large_test()3430 fn isub_large_test() {
3431 // Overflow, both single values
3432 let mut x = Bigint { data: from_u32(&[4, 1]) };
3433 let y = Bigint { data: from_u32(&[5]) };
3434 x.isub_large(&y);
3435 assert_eq!(x.data, from_u32(&[4294967295]));
3436
3437 // No overflow, single value
3438 let mut x = Bigint { data: from_u32(&[12]) };
3439 let y = Bigint { data: from_u32(&[7]) };
3440 x.isub_large(&y);
3441 assert_eq!(x.data, from_u32(&[5]));
3442
3443 // Single carry, internal overflow
3444 let mut x = Bigint { data: from_u32(&[6, 0x80000001]) };
3445 let y = Bigint { data: from_u32(&[7]) };
3446 x.isub_large(&y);
3447 assert_eq!(x.data, from_u32(&[0xFFFFFFFF, 0x80000000]));
3448
3449 // Zeros out.
3450 let mut x = Bigint { data: from_u32(&[0xFFFFFFFF, 0x7FFFFFFF]) };
3451 let y = Bigint { data: from_u32(&[0xFFFFFFFF, 0x7FFFFFFF]) };
3452 x.isub_large(&y);
3453 assert_eq!(x.data, from_u32(&[]));
3454
3455 // 1st overflows, 2nd doesn't.
3456 let mut x = Bigint { data: from_u32(&[0xFFFFFFFE, 0x80000000]) };
3457 let y = Bigint { data: from_u32(&[0xFFFFFFFF, 0x7FFFFFFF]) };
3458 x.isub_large(&y);
3459 assert_eq!(x.data, from_u32(&[0xFFFFFFFF]));
3460 }
3461
3462 #[test]
imul_large_test()3463 fn imul_large_test() {
3464 // Test by empty
3465 let mut x = Bigint { data: from_u32(&[0xFFFFFFFF]) };
3466 let y = Bigint { data: from_u32(&[]) };
3467 x.imul_large(&y);
3468 assert_eq!(x.data, from_u32(&[]));
3469
3470 // Simple case
3471 let mut x = Bigint { data: from_u32(&[0xFFFFFFFF]) };
3472 let y = Bigint { data: from_u32(&[5]) };
3473 x.imul_large(&y);
3474 assert_eq!(x.data, from_u32(&[0xFFFFFFFB, 0x4]));
3475
3476 // Large u32, but still just as easy.
3477 let mut x = Bigint { data: from_u32(&[0xFFFFFFFF]) };
3478 let y = Bigint { data: from_u32(&[0xFFFFFFFE]) };
3479 x.imul_large(&y);
3480 assert_eq!(x.data, from_u32(&[0x2, 0xFFFFFFFD]));
3481
3482 // Let's multiply two large values together
3483 let mut x = Bigint { data: from_u32(&[0xFFFFFFFE, 0x0FFFFFFF, 1]) };
3484 let y = Bigint { data: from_u32(&[0x99999999, 0x99999999, 0xCCCD9999, 0xCCCC]) };
3485 x.imul_large(&y);
3486 assert_eq!(x.data, from_u32(&[0xCCCCCCCE, 0x5CCCCCCC, 0x9997FFFF, 0x33319999, 0x999A7333, 0xD999]));
3487 }
3488
3489 #[test]
imul_karatsuba_mul_test()3490 fn imul_karatsuba_mul_test() {
3491 // Test cases triggered to use `karatsuba_mul`.
3492 let mut x = Bigint { data: from_u32(&[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16]) };
3493 let y = Bigint { data: from_u32(&[4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19]) };
3494 x.imul_large(&y);
3495 assert_eq!(x.data, from_u32(&[4, 13, 28, 50, 80, 119, 168, 228, 300, 385, 484, 598, 728, 875, 1040, 1224, 1340, 1435, 1508, 1558, 1584, 1585, 1560, 1508, 1428, 1319, 1180, 1010, 808, 573, 304]));
3496
3497 // Test cases to use karatsuba_uneven_mul
3498 let mut x = Bigint { data: from_u32(&[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16]) };
3499 let y = Bigint { data: from_u32(&[4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37]) };
3500 x.imul_large(&y);
3501 assert_eq!(x.data, from_u32(&[4, 13, 28, 50, 80, 119, 168, 228, 300, 385, 484, 598, 728, 875, 1040, 1224, 1360, 1496, 1632, 1768, 1904, 2040, 2176, 2312, 2448, 2584, 2720, 2856, 2992, 3128, 3264, 3400, 3536, 3672, 3770, 3829, 3848, 3826, 3762, 3655, 3504, 3308, 3066, 2777, 2440, 2054, 1618, 1131, 592]));
3502 }
3503
3504 #[test]
idiv_large_test()3505 fn idiv_large_test() {
3506 // Simple case.
3507 let mut x = Bigint { data: from_u32(&[0xFFFFFFFF]) };
3508 let y = Bigint { data: from_u32(&[5]) };
3509 let rem = x.idiv_large(&y);
3510 assert_eq!(x.data, from_u32(&[0x33333333]));
3511 assert_eq!(rem.data, from_u32(&[0]));
3512
3513 // Two integer case
3514 let mut x = Bigint { data: from_u32(&[0x2, 0xFFFFFFFF]) };
3515 let y = Bigint { data: from_u32(&[0xFFFFFFFE]) };
3516 let rem = x.idiv_large(&y);
3517 assert_eq!(x.data, from_u32(&[1, 1]));
3518 assert_eq!(rem.data, from_u32(&[4]));
3519
3520 // Larger large case
3521 let mut x = Bigint { data: from_u32(&[0xCCCCCCCF, 0x5CCCCCCC, 0x9997FFFF, 0x33319999, 0x999A7333, 0xD999]) };
3522 let y = Bigint { data: from_u32(&[0x99999999, 0x99999999, 0xCCCD9999, 0xCCCC]) };
3523 let rem = x.idiv_large(&y);
3524 assert_eq!(x.data, from_u32(&[0xFFFFFFFE, 0x0FFFFFFF, 1]));
3525 assert_eq!(rem.data, from_u32(&[1]));
3526
3527 // Extremely large cases, examples from Karatsuba multiplication.
3528 let mut x = Bigint { data: from_u32(&[4, 13, 29, 50, 80, 119, 168, 228, 300, 385, 484, 598, 728, 875, 1040, 1224, 1340, 1435, 1508, 1558, 1584, 1585, 1560, 1508, 1428, 1319, 1180, 1010, 808, 573, 304]) };
3529 let y = Bigint { data: from_u32(&[4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19]) };
3530 let rem = x.idiv_large(&y);
3531 assert_eq!(x.data, from_u32(&[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16]));
3532 assert_eq!(rem.data, from_u32(&[0, 0, 1]));
3533 }
3534
3535 #[test]
quorem_test()3536 fn quorem_test() {
3537 let mut x = Bigint::from_u128(42535295865117307932921825928971026432);
3538 let y = Bigint::from_u128(17218479456385750618067377696052635483);
3539 assert_eq!(x.quorem(&y), 2);
3540 assert_eq!(x.data, from_u32(&[1873752394, 3049207402, 3024501058, 102215382]));
3541 }
3542 }
3543