num_modular/
barrett.rs

1//! All methods that using pre-computed inverse of the modulus will be contained in this module,
2//! as it shares the idea of barrett reduction.
3
4// Version 1: Vanilla barrett reduction (for x mod n, x < n^2)
5// - Choose k = ceil(log2(n))
6// - Precompute r = floor(2^(k+1)/n)
7// - t = x - floor(x*r/2^(k+1)) * n
8// - if t > n, t -= n
9// - return t
10//
11// Version 2: Full width barrett reduction
12// - Similar to version 1 but support n up to full width
13// - Ref (u128): <https://math.stackexchange.com/a/3455956/815652>
14//
15// Version 3: Floating point barrett reduction
16// - Using floating point to store r
17// - Ref: <http://flintlib.org/doc/ulong_extras.html#c.n_mulmod_precomp>
18//
19// Version 4: "Improved division by invariant integers" by Granlund
20// - Ref: <https://gmplib.org/~tege/division-paper.pdf>
21//        <https://gmplib.org/~tege/divcnst-pldi94.pdf>
22//
23// Comparison between vanilla Barrett reduction and Montgomery reduction:
24// - Barrett reduction requires one 2k-by-k bits and one k-by-k bits multiplication while Montgomery only involves two k-by-k multiplications
25// - Extra conversion step is required for Montgomery form to get a normal integer
26// (Referece: <https://www.nayuki.io/page/barrett-reduction-algorithm>)
27//
28// The latter two versions are efficient and practical for use.
29
30use crate::reduced::{impl_reduced_binary_pow, Vanilla};
31use crate::{DivExact, DivExactAssign, ModularUnaryOps, Reducer};
32
33/// Divide a Word by a prearranged divisor.
34///
35/// Granlund, Montgomerry "Division by Invariant Integers using Multiplication"
36/// Algorithm 4.1.
37#[must_use]
38#[derive(Debug, Clone, Copy, PartialEq, Eq)]
39pub struct PreMulInv1by1<T> {
40    // Let n = ceil(log_2(divisor))
41    // 2^(n-1) < divisor <= 2^n
42    // m = floor(B * 2^n / divisor) + 1 - B, where B = 2^N
43    m: T,
44
45    // shift = n - 1
46    shift: u32,
47}
48
49macro_rules! impl_premulinv_1by1_for {
50    ($T:ty) => {
51        impl PreMulInv1by1<$T> {
52            pub const fn new(divisor: $T) -> Self {
53                debug_assert!(divisor > 1);
54
55                // n = ceil(log2(divisor))
56                let n = <$T>::BITS - (divisor - 1).leading_zeros();
57
58                /* Calculate:
59                 * m = floor(B * 2^n / divisor) + 1 - B
60                 * m >= B + 1 - B >= 1
61                 * m <= B * 2^n / (2^(n-1) + 1) + 1 - B
62                 *    = (B * 2^n + 2^(n-1) + 1) / (2^(n-1) + 1) - B
63                 *    = B * (2^n + 2^(n-1-N) + 2^-N) / (2^(n-1)+1) - B
64                 *    < B * (2^n + 2^1) / (2^(n-1)+1) - B
65                 *    = B
66                 * So m fits in a Word.
67                 *
68                 * Note:
69                 * divisor * (B + m) = divisor * floor(B * 2^n / divisor + 1)
70                 * = B * 2^n + k, 1 <= k <= divisor
71                 */
72
73                // m = floor(B * (2^n-1 - (divisor-1)) / divisor) + 1
74                let (lo, _hi) = split(merge(0, ones(n) - (divisor - 1)) / extend(divisor));
75                debug_assert!(_hi == 0);
76                Self {
77                    shift: n - 1,
78                    m: lo + 1,
79                }
80            }
81
82            /// (a / divisor, a % divisor)
83            #[inline]
84            pub const fn div_rem(&self, a: $T, d: $T) -> ($T, $T) {
85                // q = floor( (B + m) * a / (B * 2^n) )
86                /*
87                 * Remember that divisor * (B + m) = B * 2^n + k, 1 <= k <= 2^n
88                 *
89                 * (B + m) * a / (B * 2^n)
90                 * = a / divisor * (B * 2^n + k) / (B * 2^n)
91                 * = a / divisor + k * a / (divisor * B * 2^n)
92                 * On one hand, this is >= a / divisor
93                 * On the other hand, this is:
94                 * <= a / divisor + 2^n * (B-1) / (2^n * B) / divisor
95                 * < (a + 1) / divisor
96                 *
97                 * Therefore the floor is always the exact quotient.
98                 */
99
100                // t = m * n / B
101                let (_, t) = split(wmul(self.m, a));
102                // q = (t + a) / 2^n = (t + (a - t)/2) / 2^(n-1)
103                let q = (t + ((a - t) >> 1)) >> self.shift;
104                let r = a - q * d;
105                (q, r)
106            }
107        }
108
109        impl DivExact<$T, PreMulInv1by1<$T>> for $T {
110            type Output = $T;
111
112            #[inline]
113            fn div_exact(self, d: $T, pre: &PreMulInv1by1<$T>) -> Option<Self::Output> {
114                let (q, r) = pre.div_rem(self, d);
115                if r == 0 {
116                    Some(q)
117                } else {
118                    None
119                }
120            }
121        }
122
123        impl DivExactAssign<$T, PreMulInv1by1<$T>> for $T {
124            #[inline]
125            fn div_exact_assign(&mut self, d: $T, pre: &PreMulInv1by1<$T>) -> bool {
126                match DivExact::div_exact(*self, d, pre) {
127                    Some(q) => {
128                        *self = q;
129                        true
130                    }
131                    None => false,
132                }
133            }
134        }
135    };
136}
137
138/// Divide a DoubleWord by a prearranged divisor.
139///
140/// Assumes quotient fits in a Word.
141///
142/// Möller, Granlund, "Improved division by invariant integers", Algorithm 4.
143#[must_use]
144#[derive(Debug, Clone, Copy, PartialEq, Eq)]
145pub struct Normalized2by1Divisor<T> {
146    // Normalized (top bit must be set).
147    divisor: T,
148
149    // floor((B^2 - 1) / divisor) - B, where B = 2^T::BITS
150    m: T,
151}
152
153macro_rules! impl_normdiv_2by1_for {
154    ($T:ty, $D:ty) => {
155        impl Normalized2by1Divisor<$T> {
156            /// Calculate the inverse m > 0 of a normalized divisor (fit in a word), such that
157            ///
158            /// (m + B) * divisor = B^2 - k for some 1 <= k <= divisor
159            ///
160            #[inline]
161            pub const fn invert_word(divisor: $T) -> $T {
162                let (m, _hi) = split(<$D>::MAX / extend(divisor));
163                debug_assert!(_hi == 1);
164                m
165            }
166
167            /// Initialize from a given normalized divisor.
168            ///
169            /// The divisor must have top bit of 1
170            #[inline]
171            pub const fn new(divisor: $T) -> Self {
172                assert!(divisor.leading_zeros() == 0);
173                Self {
174                    divisor,
175                    m: Self::invert_word(divisor),
176                }
177            }
178
179            /// Returns (a / divisor, a % divisor)
180            #[inline]
181            pub const fn div_rem_1by1(&self, a: $T) -> ($T, $T) {
182                if a < self.divisor {
183                    (0, a)
184                } else {
185                    (1, a - self.divisor) // because self.divisor is normalized
186                }
187            }
188
189            /// Returns (a / divisor, a % divisor)
190            /// The result must fit in a single word.
191            #[inline]
192            pub const fn div_rem_2by1(&self, a: $D) -> ($T, $T) {
193                let (a_lo, a_hi) = split(a);
194                debug_assert!(a_hi < self.divisor);
195
196                // Approximate quotient is (m + B) * a / B^2 ~= (m * a/B + a)/B.
197                // This is q1 below.
198                // This doesn't overflow because a_hi < self.divisor <= Word::MAX.
199                let (q0, q1) = split(wmul(self.m, a_hi) + a);
200
201                // q = q1 + 1 is our first approximation, but calculate mod B.
202                // r = a - q * d
203                let q = q1.wrapping_add(1);
204                let r = a_lo.wrapping_sub(q.wrapping_mul(self.divisor));
205
206                /* Theorem: max(-d, q0+1-B) <= r < max(B-d, q0)
207                 * Proof:
208                 * r = a - q * d = a - q1 * d - d
209                 * = a - (q1 * B + q0 - q0) * d/B - d
210                 * = a - (m * a_hi + a - q0) * d/B - d
211                 * = a - ((m+B) * a_hi + a_lo - q0) * d/B - d
212                 * = a - ((B^2-k)/d * a_hi + a_lo - q0) * d/B - d
213                 * = a - B * a_hi + (a_hi * k - a_lo * d + q0 * d) / B - d
214                 * = (a_hi * k + a_lo * (B - d) + q0 * d) / B - d
215                 *
216                 * r >= q0 * d / B - d
217                 * r >= -d
218                 * r >= d/B (q0 - B) > q0-B
219                 * r >= max(-d, q0+1-B)
220                 *
221                 * r < (d * d + B * (B-d) + q0 * d) / B - d
222                 * = (B-d)^2 / B + q0 * d / B
223                 * = (1 - d/B) * (B-d) + (d/B) * q0
224                 * <= max(B-d, q0)
225                 * QED
226                 */
227
228                // if r mod B > q0 { q -= 1; r += d; }
229                //
230                // Consider two cases:
231                // a) r >= 0:
232                // Then r = r mod B > q0, hence r < B-d. Adding d will not overflow r.
233                // b) r < 0:
234                // Then r mod B = r-B > q0, and r >= -d, so adding d will make r non-negative.
235                // In either case, this will result in 0 <= r < B.
236
237                // In a branch-free way:
238                // decrease = 0xffff.fff = -1 if r mod B > q0, 0 otherwise.
239                let (_, decrease) = split(extend(q0).wrapping_sub(extend(r)));
240                let mut q = q.wrapping_add(decrease);
241                let mut r = r.wrapping_add(decrease & self.divisor);
242
243                // At this point 0 <= r < B, i.e. 0 <= r < 2d.
244                // the following fix step is unlikely to happen
245                if r >= self.divisor {
246                    q += 1;
247                    r -= self.divisor;
248                }
249
250                (q, r)
251            }
252        }
253    };
254}
255
256/// A wrapper of [Normalized2by1Divisor] that can be used as a [Reducer]
257#[must_use]
258#[derive(Debug, Clone, Copy, PartialEq, Eq)]
259pub struct PreMulInv2by1<T> {
260    div: Normalized2by1Divisor<T>,
261    shift: u32,
262}
263
264impl<T> PreMulInv2by1<T> {
265    #[inline]
266    pub const fn divider(&self) -> &Normalized2by1Divisor<T> {
267        &self.div
268    }
269    #[inline]
270    pub const fn shift(&self) -> u32 {
271        self.shift
272    }
273}
274
275macro_rules! impl_premulinv_2by1_reducer_for {
276    ($T:ty) => {
277        impl PreMulInv2by1<$T> {
278            #[inline]
279            pub const fn new(divisor: $T) -> Self {
280                let shift = divisor.leading_zeros();
281                let div = Normalized2by1Divisor::<$T>::new(divisor << shift);
282                Self { div, shift }
283            }
284
285            /// Get the **normalized** divisor.
286            #[inline]
287            pub const fn divisor(&self) -> $T {
288                self.div.divisor
289            }
290        }
291
292        impl Reducer<$T> for PreMulInv2by1<$T> {
293            #[inline]
294            fn new(m: &$T) -> Self {
295                PreMulInv2by1::<$T>::new(*m)
296            }
297            #[inline]
298            fn transform(&self, target: $T) -> $T {
299                if self.shift == 0 {
300                    self.div.div_rem_1by1(target).1
301                } else {
302                    self.div.div_rem_2by1(extend(target) << self.shift).1
303                }
304            }
305            #[inline]
306            fn check(&self, target: &$T) -> bool {
307                *target < self.div.divisor && target & ones(self.shift) == 0
308            }
309            #[inline]
310            fn residue(&self, target: $T) -> $T {
311                target >> self.shift
312            }
313            #[inline]
314            fn modulus(&self) -> $T {
315                self.div.divisor >> self.shift
316            }
317            #[inline]
318            fn is_zero(&self, target: &$T) -> bool {
319                *target == 0
320            }
321
322            #[inline(always)]
323            fn add(&self, lhs: &$T, rhs: &$T) -> $T {
324                Vanilla::<$T>::add(&self.div.divisor, *lhs, *rhs)
325            }
326            #[inline(always)]
327            fn dbl(&self, target: $T) -> $T {
328                Vanilla::<$T>::dbl(&self.div.divisor, target)
329            }
330            #[inline(always)]
331            fn sub(&self, lhs: &$T, rhs: &$T) -> $T {
332                Vanilla::<$T>::sub(&self.div.divisor, *lhs, *rhs)
333            }
334            #[inline(always)]
335            fn neg(&self, target: $T) -> $T {
336                Vanilla::<$T>::neg(&self.div.divisor, target)
337            }
338
339            #[inline(always)]
340            fn inv(&self, target: $T) -> Option<$T> {
341                self.residue(target)
342                    .invm(&self.modulus())
343                    .map(|v| v << self.shift)
344            }
345            #[inline]
346            fn mul(&self, lhs: &$T, rhs: &$T) -> $T {
347                self.div.div_rem_2by1(wmul(lhs >> self.shift, *rhs)).1
348            }
349            #[inline]
350            fn sqr(&self, target: $T) -> $T {
351                self.div.div_rem_2by1(wsqr(target) >> self.shift).1
352            }
353
354            impl_reduced_binary_pow!($T);
355        }
356    };
357}
358
359/// Divide a 3-Word by a prearranged DoubleWord divisor.
360///
361/// Assumes quotient fits in a Word.
362///
363/// Möller, Granlund, "Improved division by invariant integers"
364/// Algorithm 5.
365#[must_use]
366#[derive(Debug, Clone, Copy, PartialEq, Eq)]
367pub struct Normalized3by2Divisor<T, D> {
368    // Top bit must be 1.
369    divisor: D,
370
371    // floor ((B^3 - 1) / divisor) - B, where B = 2^WORD_BITS
372    m: T,
373}
374
375macro_rules! impl_normdiv_3by2_for {
376    ($T:ty, $D:ty) => {
377        impl Normalized3by2Divisor<$T, $D> {
378            /// Calculate the inverse m > 0 of a normalized divisor (fit in a DoubleWord), such that
379            ///
380            /// (m + B) * divisor = B^3 - k for some 1 <= k <= divisor
381            ///
382            /// Möller, Granlund, "Improved division by invariant integers", Algorithm 6.
383            #[inline]
384            pub const fn invert_double_word(divisor: $D) -> $T {
385                let (d0, d1) = split(divisor);
386                let mut v = Normalized2by1Divisor::<$T>::invert_word(d1);
387                // then B^2 - d1 <= (B + v)d1 < B^2
388
389                let (mut p, c) = d1.wrapping_mul(v).overflowing_add(d0);
390                if c {
391                    v -= 1;
392                    if p >= d1 {
393                        v -= 1;
394                        p -= d1;
395                    }
396                    p = p.wrapping_sub(d1);
397                }
398                // then B^2 - d1 <= (B + v)d1 + d0 < B^2
399
400                let (t0, t1) = split(extend(v) * extend(d0));
401                let (p, c) = p.overflowing_add(t1);
402                if c {
403                    v -= 1;
404                    if merge(t0, p) >= divisor {
405                        v -= 1;
406                    }
407                }
408
409                v
410            }
411
412            /// Initialize from a given normalized divisor.
413            ///
414            /// divisor must have top bit of 1
415            #[inline]
416            pub const fn new(divisor: $D) -> Self {
417                assert!(divisor.leading_zeros() == 0);
418                Self {
419                    divisor,
420                    m: Self::invert_double_word(divisor),
421                }
422            }
423
424            #[inline]
425            pub const fn div_rem_2by2(&self, a: $D) -> ($D, $D) {
426                if a < self.divisor {
427                    (0, a)
428                } else {
429                    (1, a - self.divisor) // because self.divisor is normalized
430                }
431            }
432
433            /// The input a is arranged as (lo, mi & hi)
434            /// The output is (a / divisor, a % divisor)
435            pub const fn div_rem_3by2(&self, a_lo: $T, a_hi: $D) -> ($T, $D) {
436                debug_assert!(a_hi < self.divisor);
437                let (a1, a2) = split(a_hi);
438                let (d0, d1) = split(self.divisor);
439
440                // This doesn't overflow because a2 <= self.divisor / B <= Word::MAX.
441                let (q0, q1) = split(wmul(self.m, a2) + a_hi);
442                let r1 = a1.wrapping_sub(q1.wrapping_mul(d1));
443                let t = wmul(d0, q1);
444                let r = merge(a_lo, r1).wrapping_sub(t).wrapping_sub(self.divisor);
445
446                // The first guess of quotient is q1 + 1
447                // if r1 >= q0 { r += d; } else { q1 += 1; }
448                // In a branch-free way:
449                // decrease = 0 if r1 >= q0, = 0xffff.fff = -1 otherwise
450                let (_, r1) = split(r);
451                let (_, decrease) = split(extend(r1).wrapping_sub(extend(q0)));
452                let mut q1 = q1.wrapping_sub(decrease);
453                let mut r = r.wrapping_add(merge(!decrease, !decrease) & self.divisor);
454
455                // the following fix step is unlikely to happen
456                if r >= self.divisor {
457                    q1 += 1;
458                    r -= self.divisor;
459                }
460
461                (q1, r)
462            }
463
464            /// Divide a 4-word number with double word divisor
465            ///
466            /// The output is (a / divisor, a % divisor)
467            pub const fn div_rem_4by2(&self, a_lo: $D, a_hi: $D) -> ($D, $D) {
468                let (a0, a1) = split(a_lo);
469                let (q1, r1) = self.div_rem_3by2(a1, a_hi);
470                let (q0, r0) = self.div_rem_3by2(a0, r1);
471                (merge(q0, q1), r0)
472            }
473        }
474    };
475}
476
477/// A wrapper of [Normalized3by2Divisor] that can be used as a [Reducer]
478#[must_use]
479#[derive(Debug, Clone, Copy, PartialEq, Eq)]
480pub struct PreMulInv3by2<T, D> {
481    div: Normalized3by2Divisor<T, D>,
482    shift: u32,
483}
484
485impl<T, D> PreMulInv3by2<T, D> {
486    #[inline]
487    pub const fn divider(&self) -> &Normalized3by2Divisor<T, D> {
488        &self.div
489    }
490    #[inline]
491    pub const fn shift(&self) -> u32 {
492        self.shift
493    }
494}
495
496macro_rules! impl_premulinv_3by2_reducer_for {
497    ($T:ty, $D:ty) => {
498        impl PreMulInv3by2<$T, $D> {
499            #[inline]
500            pub const fn new(divisor: $D) -> Self {
501                let shift = divisor.leading_zeros();
502                let div = Normalized3by2Divisor::<$T, $D>::new(divisor << shift);
503                Self { div, shift }
504            }
505
506            /// Get the **normalized** divisor.
507            #[inline]
508            pub const fn divisor(&self) -> $D {
509                self.div.divisor
510            }
511        }
512
513        impl Reducer<$D> for PreMulInv3by2<$T, $D> {
514            #[inline]
515            fn new(m: &$D) -> Self {
516                assert!(*m > <$T>::MAX as $D);
517                let shift = m.leading_zeros();
518                let div = Normalized3by2Divisor::<$T, $D>::new(m << shift);
519                Self { div, shift }
520            }
521            #[inline]
522            fn transform(&self, target: $D) -> $D {
523                if self.shift == 0 {
524                    self.div.div_rem_2by2(target).1
525                } else {
526                    let (lo, hi) = split(target);
527                    let (n0, carry) = split(extend(lo) << self.shift);
528                    let n12 = (extend(hi) << self.shift) | extend(carry);
529                    self.div.div_rem_3by2(n0, n12).1
530                }
531            }
532            #[inline]
533            fn check(&self, target: &$D) -> bool {
534                *target < self.div.divisor && split(*target).0 & ones(self.shift) == 0
535            }
536            #[inline]
537            fn residue(&self, target: $D) -> $D {
538                target >> self.shift
539            }
540            #[inline]
541            fn modulus(&self) -> $D {
542                self.div.divisor >> self.shift
543            }
544            #[inline]
545            fn is_zero(&self, target: &$D) -> bool {
546                *target == 0
547            }
548
549            #[inline(always)]
550            fn add(&self, lhs: &$D, rhs: &$D) -> $D {
551                Vanilla::<$D>::add(&self.div.divisor, *lhs, *rhs)
552            }
553            #[inline(always)]
554            fn dbl(&self, target: $D) -> $D {
555                Vanilla::<$D>::dbl(&self.div.divisor, target)
556            }
557            #[inline(always)]
558            fn sub(&self, lhs: &$D, rhs: &$D) -> $D {
559                Vanilla::<$D>::sub(&self.div.divisor, *lhs, *rhs)
560            }
561            #[inline(always)]
562            fn neg(&self, target: $D) -> $D {
563                Vanilla::<$D>::neg(&self.div.divisor, target)
564            }
565
566            #[inline(always)]
567            fn inv(&self, target: $D) -> Option<$D> {
568                self.residue(target)
569                    .invm(&self.modulus())
570                    .map(|v| v << self.shift)
571            }
572            #[inline]
573            fn mul(&self, lhs: &$D, rhs: &$D) -> $D {
574                let prod = DoubleWordModule::wmul(lhs >> self.shift, *rhs);
575                let (lo, hi) = DoubleWordModule::split(prod);
576                self.div.div_rem_4by2(lo, hi).1
577            }
578            #[inline]
579            fn sqr(&self, target: $D) -> $D {
580                let prod = DoubleWordModule::wsqr(target) >> self.shift;
581                let (lo, hi) = DoubleWordModule::split(prod);
582                self.div.div_rem_4by2(lo, hi).1
583            }
584
585            impl_reduced_binary_pow!($D);
586        }
587    };
588}
589
590macro_rules! collect_impls {
591    ($T:ident, $ns:ident) => {
592        mod $ns {
593            use super::*;
594            use crate::word::$T::*;
595
596            impl_premulinv_1by1_for!(Word);
597            impl_normdiv_2by1_for!(Word, DoubleWord);
598            impl_premulinv_2by1_reducer_for!(Word);
599            impl_normdiv_3by2_for!(Word, DoubleWord);
600            impl_premulinv_3by2_reducer_for!(Word, DoubleWord);
601        }
602    };
603}
604collect_impls!(u8, u8_impl);
605collect_impls!(u16, u16_impl);
606collect_impls!(u32, u32_impl);
607collect_impls!(u64, u64_impl);
608collect_impls!(usize, usize_impl);
609
610#[cfg(test)]
611mod tests {
612    use super::*;
613    use crate::reduced::tests::ReducedTester;
614    use rand::prelude::*;
615
616    #[test]
617    #[allow(unstable_name_collisions)]
618    fn test_mul_inv_1by1() {
619        type Word = u64;
620        let mut rng = StdRng::seed_from_u64(1);
621        for _ in 0..400000 {
622            let d_bits = rng.gen_range(2..=Word::BITS);
623            let max_d = Word::MAX >> (Word::BITS - d_bits);
624            let d = rng.gen_range(max_d / 2 + 1..=max_d);
625            let fast_div = PreMulInv1by1::<Word>::new(d);
626            let n = rng.gen();
627            let (q, r) = fast_div.div_rem(n, d);
628            assert_eq!(q, n / d);
629            assert_eq!(r, n % d);
630
631            if r == 0 {
632                assert_eq!(n.div_exact(d, &fast_div), Some(q));
633            } else {
634                assert_eq!(n.div_exact(d, &fast_div), None);
635            }
636
637            let mut n2 = n;
638            if r == 0 {
639                assert!(n2.div_exact_assign(d, &fast_div));
640                assert_eq!(n2, q);
641            } else {
642                assert!(!n2.div_exact_assign(d, &fast_div));
643                assert_eq!(n2, n);
644            }
645        }
646    }
647
648    #[test]
649    fn test_mul_inv_2by1() {
650        type Word = u64;
651        type Divider = Normalized2by1Divisor<Word>;
652        use crate::word::u64::*;
653
654        let fast_div = Divider::new(Word::MAX);
655        assert_eq!(fast_div.div_rem_2by1(0), (0, 0));
656
657        let mut rng = StdRng::seed_from_u64(1);
658        for _ in 0..200000 {
659            let d = rng.gen_range(Word::MAX / 2 + 1..=Word::MAX);
660            let q = rng.gen();
661            let r = rng.gen_range(0..d);
662            let (a0, a1) = split(wmul(q, d) + extend(r));
663            let fast_div = Divider::new(d);
664            assert_eq!(fast_div.div_rem_2by1(merge(a0, a1)), (q, r));
665        }
666    }
667
668    #[test]
669    fn test_mul_inv_3by2() {
670        type Word = u64;
671        type DoubleWord = u128;
672        type Divider = Normalized3by2Divisor<Word, DoubleWord>;
673        use crate::word::u64::*;
674
675        let d = DoubleWord::MAX;
676        let fast_div = Divider::new(d);
677        assert_eq!(fast_div.div_rem_3by2(0, 0), (0, 0));
678
679        let mut rng = StdRng::seed_from_u64(1);
680        for _ in 0..100000 {
681            let d = rng.gen_range(DoubleWord::MAX / 2 + 1..=DoubleWord::MAX);
682            let r = rng.gen_range(0..d);
683            let q = rng.gen();
684
685            let (d0, d1) = split(d);
686            let (r0, r1) = split(r);
687            let (a0, c) = split(wmul(q, d0) + extend(r0));
688            let (a1, a2) = split(wmul(q, d1) + extend(r1) + extend(c));
689            let a12 = merge(a1, a2);
690
691            let fast_div = Divider::new(d);
692            assert_eq!(
693                fast_div.div_rem_3by2(a0, a12),
694                (q, r),
695                "failed at {:?} / {}",
696                (a0, a12),
697                d
698            );
699        }
700    }
701
702    #[test]
703    fn test_mul_inv_4by2() {
704        type Word = u64;
705        type DoubleWord = u128;
706        type Divider = Normalized3by2Divisor<Word, DoubleWord>;
707        use crate::word::u128::*;
708
709        let mut rng = StdRng::seed_from_u64(1);
710        for _ in 0..20000 {
711            let d = rng.gen_range(DoubleWord::MAX / 2 + 1..=DoubleWord::MAX);
712            let q = rng.gen();
713            let r = rng.gen_range(0..d);
714            let (a_lo, a_hi) = split(wmul(q, d) + r as DoubleWord);
715            let fast_div = Divider::new(d);
716            assert_eq!(fast_div.div_rem_4by2(a_lo, a_hi), (q, r));
717        }
718    }
719
720    #[test]
721    fn test_2by1_against_modops() {
722        for _ in 0..10 {
723            ReducedTester::<u8>::test_against_modops::<PreMulInv2by1<u8>>(0);
724            ReducedTester::<u16>::test_against_modops::<PreMulInv2by1<u16>>(0);
725            ReducedTester::<u32>::test_against_modops::<PreMulInv2by1<u32>>(0);
726            ReducedTester::<u64>::test_against_modops::<PreMulInv2by1<u64>>(0);
727            // ReducedTester::<u128>::test_against_modops::<PreMulInv2by1<u128>>();
728            ReducedTester::<usize>::test_against_modops::<PreMulInv2by1<usize>>(0);
729        }
730    }
731
732    #[test]
733    fn test_3by2_against_modops() {
734        for _ in 0..10 {
735            ReducedTester::<u16>::test_against_modops::<PreMulInv3by2<u8, u16>>(2);
736            ReducedTester::<u32>::test_against_modops::<PreMulInv3by2<u16, u32>>(2);
737            ReducedTester::<u64>::test_against_modops::<PreMulInv3by2<u32, u64>>(2);
738            ReducedTester::<u128>::test_against_modops::<PreMulInv3by2<u64, u128>>(2);
739        }
740    }
741}
num_modular/barrett.rs

num_modular/
barrett.rs
