ruint/algorithms/div/
mod.rs

Help
1//! ⚠️ Collection of division algorithms.
2//!
3//! **Warning.** Most functions in this module are currently not considered part
4//! of the stable API and may be changed or removed in future minor releases.
5//!
6//! All division algorithms also compute the remainder. There is no benefit
7//! to splitting the division and remainder into separate functions, since
8//! the remainder is always computed as part of the division algorithm.
9//!
10//! These functions are adapted from algorithms in [MG10] and [K97].
11//!
12//! [K97]: https://cs.stanford.edu/~knuth/taocp.html
13//! [MG10]: https://gmplib.org/~tege/division-paper.pdf
14
15#![allow(clippy::similar_names)] // TODO
16
17mod knuth;
18mod reciprocal;
19mod small;
20
21pub use self::{
22    knuth::{div_nxm, div_nxm_normalized},
23    reciprocal::{reciprocal, reciprocal_2, reciprocal_2_mg10, reciprocal_mg10, reciprocal_ref},
24    small::{
25        div_2x1, div_2x1_mg10, div_2x1_ref, div_3x2, div_3x2_mg10, div_3x2_ref, div_nx1,
26        div_nx1_normalized, div_nx2, div_nx2_normalized,
27    },
28};
29use crate::algorithms::DoubleWord;
30
31/// ⚠️ Division with remainder.
32///
33/// **Warning.** This function is not part of the stable API.
34///
35/// The quotient is stored in the `numerator` and the remainder is stored
36/// in the `divisor`.
37///
38/// # Algorithm
39///
40/// It trims zeros from the numerator and divisor then solves the trivial cases
41/// directly, or dispatches to the [`div_nx1`], [`div_nx2`] or [`div_nxm`]
42/// functions.
43///
44/// # Panics
45///
46/// Panics if `divisor` is zero.
47#[inline]
48pub fn div(numerator: &mut [u64], divisor: &mut [u64]) {
49    // Trim most significant zeros from divisor.
50    let i = divisor
51        .iter()
52        .rposition(|&x| x != 0)
53        .expect("Divisor is zero");
54    let divisor = &mut divisor[..=i];
55    debug_assert!(!divisor.is_empty());
56    debug_assert!(divisor.last() != Some(&0));
57
58    // Trim zeros from numerator
59    let numerator = if let Some(i) = numerator.iter().rposition(|&n| n != 0) {
60        &mut numerator[..=i]
61    } else {
62        // Empty numerator (q, r) = (0,0)
63        divisor.fill(0);
64        return;
65    };
66    debug_assert!(!numerator.is_empty());
67    debug_assert!(*numerator.last().unwrap() != 0);
68
69    // If numerator is smaller than divisor (q, r) = (0, numerator)
70    if numerator.len() < divisor.len() {
71        let (remainder, padding) = divisor.split_at_mut(numerator.len());
72        remainder.copy_from_slice(numerator);
73        padding.fill(0);
74        numerator.fill(0);
75        return;
76    }
77    debug_assert!(numerator.len() >= divisor.len());
78
79    // Compute quotient and remainder, branching out to different algorithms.
80    if divisor.len() <= 2 {
81        if divisor.len() == 1 {
82            if numerator.len() == 1 {
83                let q = numerator[0] / divisor[0];
84                let r = numerator[0] % divisor[0];
85                numerator[0] = q;
86                divisor[0] = r;
87            } else {
88                divisor[0] = div_nx1(numerator, divisor[0]);
89            }
90        } else {
91            let d = u128::join(divisor[1], divisor[0]);
92            let remainder = div_nx2(numerator, d);
93            divisor[0] = remainder.low();
94            divisor[1] = remainder.high();
95        }
96    } else {
97        div_nxm(numerator, divisor);
98    }
99}
100
101#[cfg(test)]
102mod tests {
103    use super::*;
104    use crate::aliases::U512;
105
106    // Test vectors from <https://github.com/chfast/intx/blob/8b5f4748a7386a9530769893dae26b3273e0ffe2/test/unittests/test_div.cpp#L58>
107    // [[numerator, divisor, quotient, remainder]; _]
108    const INTX_TESTS: [[U512; 4]; 45] = uint!([
109        [2_U512, 1_U512, 2_U512, 0_U512],
110        [
111            0x10000000000000000_U512,
112            2_U512,
113            0x8000000000000000_U512,
114            0_U512,
115        ],
116        [
117            0x7000000000000000_U512,
118            0x8000000000000000_U512,
119            0_U512,
120            0x7000000000000000_U512,
121        ],
122        [
123            0x8000000000000000_U512,
124            0x8000000000000000_U512,
125            1_U512,
126            0_U512,
127        ],
128        [
129            0x8000000000000001_U512,
130            0x8000000000000000_U512,
131            1_U512,
132            1_U512,
133        ],
134        [
135            0x80000000000000010000000000000000_U512,
136            0x80000000000000000000000000000000_U512,
137            1_U512,
138            0x10000000000000000_U512,
139        ],
140        [
141            0x80000000000000000000000000000000_U512,
142            0x80000000000000000000000000000001_U512,
143            0_U512,
144            0x80000000000000000000000000000000_U512,
145        ],
146        [
147            0x478392145435897052_U512,
148            0x111_U512,
149            0x430f89ebadad0baa_U512,
150            8_U512,
151        ],
152        [
153            0x400000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000_U512,
154            0x800000000000000000000000000000000000000000000000_U512,
155            0x800000000000000000000000000000000000000000000000_U512,
156            0_U512,
157        ],
158        [
159            0x80000000000000000000000000000000000000000000000000000000000000000000000000000000_U512,
160            0x800000000000000000000000000000000000000000000000_U512,
161            0x100000000000000000000000000000000_U512,
162            0_U512,
163        ],
164        [
165            0x1e00000000000000000000090000000000000000000000000000000000000000000000000000000000000000000000000000000009000000000000000000_U512,
166            0xa_U512,
167            0x30000000000000000000000e6666666666666666666666666666666666666666666666666666666666666666666666666666666674ccccccccccccccccc_U512,
168            8_U512,
169        ],
170        [
171            0x767676767676767676000000767676767676_U512,
172            0x2900760076761e00020076760000000076767676000000_U512,
173            0_U512,
174            0x767676767676767676000000767676767676_U512,
175        ],
176        [
177            0x12121212121212121212121212121212_U512,
178            0x232323232323232323_U512,
179            0x83a83a83a83a83_U512,
180            0x171729292929292929_U512,
181        ],
182        [
183            0xabc0abc0abc0abc0abc0abc0abc0abc0abc0abc0abc0abc0abc0abc0abc0abc0abc0abc0abc0abc0abc0abc0abc0abc0abc0abc0abc0abc0abc0abc0abc0_U512,
184            0x1c01c01c01c01c01c01c01c01c_U512,
185            0x621ed21ed21ed21ed21ed21ed224f40bf40bf40bf40bf40bf40bf46e12de12de12de12de12de12de1900000000000000000_U512,
186            0xabc0abc0abc0abc0_U512,
187        ],
188        [
189            0xfffff716b61616160b0b0b2b0b0b0becf4bef50a0df4f48b090b2b0bc60a0a00_U512,
190            0xfffff716b61616160b0b0b2b0b230b000008010d0a2b00_U512,
191            0xffffffffffffffffff_U512,
192            0xfffff7169e17030ac1ff082ed51796090b330cd3143500_U512,
193        ],
194        [
195            0x50beb1c60141a0000dc2b0b0b0b0b0b410a0a0df4f40b090b2b0bc60a0a00_U512,
196            0x2000110000000d0a300e750a000000090a0a_U512,
197            0x285f437064cd09ff8bc5b7857d_U512,
198            0x1fda1c384d86199e14bb4edfc6693042f11e_U512,
199        ],
200        [
201            0x4b00000b41000b0b0b2b0b0b0b0b0b410a0aeff4f40b090b2b0bc60a0a1000_U512,
202            0x4b00000b41000b0b0b2b0b0b0b0b0b410a0aeff4f40b0a0a_U512,
203            0xffffffffffffff_U512,
204            0x4b00000b41000b0b0b2b0b0b0b0b0b400b35fbbafe151a0a_U512,
205        ],
206        [
207            0xeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeeee_U512,
208            7_U512,
209            0x22222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222222_U512,
210            0_U512,
211        ],
212        [
213            0xf6376770abd3a36b20394c5664afef1194c801c3f05e42566f085ed24d002bb0_U512,
214            0xb368d219438b7f3f_U512,
215            0x15f53bce87e9fb63c7c3ab03f6c0ba30d3ecf982fa97cdf0a_U512,
216            0x4bfd94dbec31523a_U512,
217        ],
218        [
219            0x0_U512,
220            0x10900000000000000000000000000000000000000000000000000_U512,
221            0x0_U512,
222            0x0_U512,
223        ],
224        [
225            0x77676767676760000000000000001002e000000000000040000000e000000000000007f0000000000000000000000000000000000000000f7000000000000_U512,
226            0xfffc000000000000767676240000000000002b05760476000000000000000000767676767600000000000000000000000000000000_U512,
227            0x7769450c7b994e65025_U512,
228            0x241cb1aa4f67c22ae65c9920bf3bb7ad8280311a887aee8be4054a3e242a5ea9ab35d800f2000000000000000000f7000000000000_U512,
229        ],
230        [
231            0xdffffffffffffffffffffffffffffffffff00000000000000000000000000000000000000000001000000000000000000000000008100000000001001_U512,
232            0xdffffffffffffffffffffffffffffffffffffffffffff3fffffffffffffffffffffffffffff_U512,
233            0xffffffffffffffffffffffffffffffffffedb6db6db6e9_U512,
234            0x200000000000000000000000000010000f2492492492ec000000000000080ffedb6db6dc6ea_U512,
235        ],
236        [
237            0xff000000000000000000000000000000000000000400000092767600000000000000000000000081000000000000000000000001020000000000eeffffffffff_U512,
238            0xffffffffffffffffffffffffffffffffffffffffffffffffffffffffff000000000000000000000005000000000000000000ffffffffff100000000000000000_U512,
239            0x0_U512,
240            0xff000000000000000000000000000000000000000400000092767600000000000000000000000081000000000000000000000001020000000000eeffffffffff_U512,
241        ],
242        [
243            0xfffffffffffffffffffffffffffffffffffffffbffffff6d8989ffffffffffffffffffffffff7efffffffffffffffffffffffefdffffffffff110000000001_U512,
244            0xfffffffffffffffffffffffaffffffffffffffffff0000000000f00000000000000000_U512,
245            0x1000000000000000000000004fffffffffffffffc00ffff8689890fff_U512,
246            0xffffffec09fffda0afa81efafc00ffff868d481fff71de0d8100efffff110000000001_U512,
247        ],
248        [
249            0x767676767676000000000076000000000000005600000000000000000000_U512,
250            0x767676767676000000000076000000760000_U512,
251            0xffffffffffffffffffffffff_U512,
252            0x767676007676005600000076000000760000_U512,
253        ],
254        [
255            0x8200000000000000000000000000000000000000000000000000000000000000_U512,
256            0x8200000000000000fe000004000000ffff000000fffff700_U512,
257            0xfffffffffffffffe_U512,
258            0x5fffffbffffff01fd00000700000afffe000001ffffee00_U512,
259        ],
260        [
261            0xdac7fff9ffd9e1322626262626262600_U512,
262            0xd021262626262626_U512,
263            0x10d1a094108c5da55_U512,
264            0x6f386ccc73c11f62_U512,
265        ],
266        [
267            0x8000000000000001800000000000000080000000000000008000000000000000_U512,
268            0x800000000000000080000000000000008000000000000000_U512,
269            0x10000000000000001_U512,
270            0x7fffffffffffffff80000000000000000000000000000000_U512,
271        ],
272        [
273            0x00e8e8e8e2000100000009ea02000000000000ff3ffffff800000010002200000000000000000000000000000000000000000000000000000000000000000000_U512,
274            0x00e8e8e8e2000100000009ea02000000000000ff3ffffff800000010002280ff0000000000000000000000000000000000000000000000000000000000000000_U512,
275            0_U512,
276            0x00e8e8e8e2000100000009ea02000000000000ff3ffffff800000010002200000000000000000000000000000000000000000000000000000000000000000000_U512,
277        ],
278        [
279            0x000000c9700000000000000000023f00c00014ff0000000000000000223008050000000000000000000000000000000000000000000000000000000000000000_U512,
280            0x00000000c9700000000000000000023f00c00014ff002c0000000000002231080000000000000000000000000000000000000000000000000000000000000000_U512,
281            0xff_U512,
282            0x00000000c9700000000000000000023f00c00014fed42c00000000000021310d0000000000000000000000000000000000000000000000000000000000000000_U512,
283        ],
284        [
285            0x40000000fd000000db00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000001_U512,
286            0x40000000fd000000db0000000000000000000040000000fd000000db000001_U512,
287            0xfffffffffffffffffffffffffffffffffffffeffffffffffffffff_U512,
288            0x3ffffffffd000000db000040000000fd0000011b000001fd000000db000002_U512,
289        ],
290        [
291            0x40000000fd000000db00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000001_U512,
292            0x40000000fd000000db0000000000000000000040000000fd000000db0000d3_U512,
293            0xfffffffffffffffffffffffffffffffffffffeffffffffffffffff_U512,
294            0x3fffff2dfd000000db000040000000fd0000011b0000d3fd000000db0000d4_U512,
295        ],
296        [
297            0x001f000000000000000000000000000000200000000100000000000000000000_U512,
298            0x0000000000000000000100000000ffffffffffffffff0000000000002e000000_U512,
299            0x1effffffe10000001f_U512,
300            0xfffa6e20000591fffffa6e000000_U512,
301        ],
302        [
303            0x7effffff80000000000000000000000000020000440000000000000000000001_U512,
304            0x7effffff800000007effffff800000008000ff0000010000_U512,
305            0xfffffffffffffffe_U512,
306            0x7effffff800000007e0100ff43ff00010001fe0000020001_U512,
307        ],
308        [
309            0x5fd8fffffffffffffffffffffffffffffc090000ce700004d0c9ffffff000001_U512,
310            0x2ffffffffffffffffffffffffffffffffff000000030000_U512,
311            0x1ff2ffffffffffffff_U512,
312            0x2fffffffffffffffc28f300ce102704d0c8ffffff030001_U512,
313        ],
314        [
315            0x62d8fffffffffffffffffffffffffffffc18000000000000000000ca00000001_U512,
316            0x2ffffffffffffffffffffffffffffffffff200000000000_U512,
317            0x20f2ffffffffffffff_U512,
318            0x2fffffffffffffffc34d49fffffffffffff20ca00000001_U512,
319        ],
320        [
321            0x7effffff8000000000000000000000000000000000000000d900000000000001_U512,
322            0x7effffff8000000000000000000000000000000000008001_U512,
323            0xffffffffffffffff_U512,
324            0x7effffff7fffffffffffffffffff7fffd900000000008002_U512,
325        ],
326        [
327            0x0000000000000006400aff20ff00200004e7fd1eff08ffca0afd1eff08ffca0a_U512,
328            0x00000000000000210000000000000022_U512,
329            0x307c7456554d945ce57749fd52bfdb7f_U512,
330            0x1491254b5a0b84a32c_U512,
331        ],
332        [
333            0x7effffff8000000000000000000000000000000000150000d900000000000001_U512,
334            0x7effffff8000000000000000000000000000000000f9e101_U512,
335            0xffffffffffffffff_U512,
336            0x7effffff7fffffffffffffffff1b1effd900000000f9e102_U512,
337        ],
338        [
339            0xffffffff0100000000000000000000000000ffff0000ffffffff0100000000_U512,
340            0xffffffff010000000000000000000000ffff0000ffffff_U512,
341            0xffffffffffffffff_U512,
342            0xffffffff00ffffff0001fffe00010100fffe0100ffffff_U512,
343        ],
344        [
345            0xabfffff0000ffffffffff36363636363636363636d00500000000ffffffffffffe90000ff00000000000000000000ffff0000000000_U512,
346            0xabfffff0000ffffffffff36363636363636363636d00500000000ffffffffffffe9ff001f_U512,
347            0xffffffffffffffffffffffffffffffffff_U512,
348            0xabfffff0000ffffffffff36363636363537371636d00500000001000000fffeffe9ff001f_U512,
349        ],
350        [
351            0xff00ffffffffffffffcaffffffff0100_U512,
352            0x0100000000000000ff800000000000ff_U512,
353            0xff_U512,
354            0xffffffffff017f4afffffffe02ff_U512,
355        ],
356        [
357            0x9000ffffffffffffffcaffffffff0100_U512,
358            0x800000000000007fc000000000007f80_U512,
359            1_U512,
360            0x1000ffffffffff803fcafffffffe8180_U512,
361        ],
362        [
363            // Very special case for reciprocal_3by2().
364            170141183460488574554024512018559533058_U512,
365            170141183460488574554024512018559533057_U512,
366            1_U512,
367            1_U512,
368        ],
369        [
370            0x6e2d23924d38f0ab643864e9b2a328a54914f48533114fae3475168bfd74a61ae91e676b4a4f33a5b3b6cc189536ccb4afc46d02b061d6daaf0298c993376ab4_U512,
371            170141183460488574554024512018559533057_U512,
372            0xdc5a47249a56560d078334729ffb61da211f5d2ec622c22f88bc3b4ebae1abdac6b03621554ef71070bc1e0dc5c301bc_U512,
373            0x6dc100ea02272bdcf68a4a5b95f468f8_U512,
374        ]
375    ]);
376
377    macro_rules! test_cases {
378        ($n:ty, $d:ty) => {
379            for [numerator, divisor, quotient, remainder] in INTX_TESTS {
380                if numerator.bit_len() > <$n>::BITS || divisor.bit_len() > <$d>::BITS {
381                    continue;
382                }
383                let mut numerator = <$n>::from(numerator).into_limbs();
384                let mut divisor = <$d>::from(divisor).into_limbs();
385                let quotient = <$n>::from(quotient).into_limbs();
386                let remainder = <$d>::from(remainder).into_limbs();
387                div(&mut numerator, &mut divisor);
388                assert_eq!(numerator, quotient);
389                assert_eq!(divisor, remainder);
390            }
391        };
392    }
393
394    #[test]
395    fn test_div_8x8() {
396        use crate::aliases::U512;
397        test_cases!(U512, U512);
398    }
399
400    #[test]
401    fn test_div_6x6() {
402        use crate::aliases::U384;
403        test_cases!(U384, U384);
404    }
405
406    #[test]
407    fn test_div_4x4() {
408        use crate::aliases::U256;
409        test_cases!(U256, U256);
410    }
411
412    #[test]
413    fn test_div_5x4() {
414        use crate::aliases::{U256, U320};
415        test_cases!(U320, U256);
416    }
417
418    #[test]
419    fn test_div_8x4() {
420        use crate::aliases::{U256, U512};
421        test_cases!(U512, U256);
422    }
423
424    #[test]
425    #[allow(clippy::unreadable_literal)]
426    fn test_div_8by4_one() {
427        let mut numerator = [
428            0x9c2bcebfa9cca2c6_u64,
429            0x274e154bb5e24f7a_u64,
430            0xe1442d5d3842be2b_u64,
431            0xf18f5adfd420853f_u64,
432            0x04ed6127eba3b594_u64,
433            0xc5c179973cdb1663_u64,
434            0x7d7f67780bb268ff_u64,
435            0x0000000000000003_u64,
436            0x0000000000000000_u64,
437        ];
438        let mut divisor = [
439            0x0181880b078ab6a1_u64,
440            0x62d67f6b7b0bda6b_u64,
441            0x92b1840f9c792ded_u64,
442            0x0000000000000019_u64,
443        ];
444        let expected_quotient = [
445            0x9128464e61d6b5b3_u64,
446            0xd9eea4fc30c5ac6c_u64,
447            0x944a2d832d5a6a08_u64,
448            0x22f06722e8d883b1_u64,
449            0x0000000000000000_u64,
450        ];
451        let expected_remainder = [
452            0x1dfa5a7ea5191b33_u64,
453            0xb5aeb3f9ad5e294e_u64,
454            0xfc710038c13e4eed_u64,
455            0x000000000000000b_u64,
456        ];
457        div(&mut numerator, &mut divisor);
458        assert_eq!(numerator[..5], expected_quotient);
459        assert_eq!(divisor, expected_remainder);
460    }
461}
ruint/algorithms/div/mod.rs

ruint/algorithms/div/
mod.rs
