num_modular/
solinas.rs

1use crate::reduced::impl_reduced_binary_pow;
2use crate::{imax, udouble, umax, ModularUnaryOps, Reducer};
3
4// REF: Handbook of Cryptography 14.3.4
5
6macro_rules! impl_fixed_trinomial_solinas {
7    (
8        $TypeName:ident,
9        $T:ty,
10        $K:ty,
11        $D:ty,
12        $half_bits:expr,
13        $max_P1:expr,
14        $kind:ident
15    ) => {
16        impl<const P1: u8, const P2: u8, const K: $K> $TypeName<P1, P2, K> {
17            const BITMASK: $T = match (1 as $T).checked_shl(P1 as u32) {
18                Some(v) => v.wrapping_sub(1),
19                None => <$T>::MAX,
20            };
21            pub const MODULUS: $T = {
22                let p1 = match (1 as $T).checked_shl(P1 as u32) {
23                    Some(v) => v,
24                    None => 0,
25                };
26                let p2 = match (1 as $T).checked_shl(P2 as u32) {
27                    Some(v) => v,
28                    None => panic!("P2 exceeds type width"),
29                };
30                if K >= 0 {
31                    p1.wrapping_sub(p2).wrapping_add(K as $T)
32                } else {
33                    p1.wrapping_sub(p2).wrapping_sub((-K) as $T)
34                }
35            };
36
37            /// Worst-case fold count for `reduce_double`.
38            /// Each fold removes roughly (P1−P2) bits; ⌈P1/(P1−P2)⌉ folds
39            /// shrink from 2·P1 bits to ≤ P1, plus 1 (K>0) or 2 (K<0) for the carry tail.
40            const FOLDS: u32 = {
41                let gap = (P1 - P2) as u32;
42                let folds_ceil = ((P1 as u32) + gap - 1) / gap;
43                if K > 0 {
44                    folds_ceil + 1
45                } else if K < 0 {
46                    folds_ceil + 2
47                } else {
48                    1 // K == 0: trivial reduction, single fold
49                }
50            };
51
52            impl_fixed_trinomial_solinas!(@reduce_single, $kind, $T, $D);
53            impl_fixed_trinomial_solinas!(@reduce_double, $kind, $T, $D);
54        }
55
56        impl<const P1: u8, const P2: u8, const K: $K> Reducer<$T> for $TypeName<P1, P2, K> {
57            #[inline]
58            fn new(m: &$T) -> Self {
59                assert!(
60                    *m == Self::MODULUS,
61                    "the given modulus doesn't match with the generic params"
62                );
63                debug_assert!(P1 <= $max_P1);
64                debug_assert!(P2 > 0 && P1 > P2);
65                debug_assert!(K % 2 != 0); // modulus must be odd
66                // |K| < 2^P2 keeps each reduction step non-negative in Z (required for unsigned arithmetic)
67                debug_assert!((K.unsigned_abs() as u128) < (1u128 << (P2 as u32)));
68                debug_assert!(
69                    (Self::MODULUS == 3 || Self::MODULUS % 3 != 0)
70                        && (Self::MODULUS == 5 || Self::MODULUS % 5 != 0)
71                        && (Self::MODULUS == 7 || Self::MODULUS % 7 != 0)
72                        && (Self::MODULUS == 11 || Self::MODULUS % 11 != 0)
73                        && (Self::MODULUS == 13 || Self::MODULUS % 13 != 0)
74                ); // error on easy composites
75                Self {}
76            }
77            #[inline]
78            fn transform(&self, target: $T) -> $T {
79                Self::reduce_single(target)
80            }
81            #[inline]
82            fn check(&self, target: &$T) -> bool {
83                *target < Self::MODULUS
84            }
85            #[inline]
86            fn residue(&self, target: $T) -> $T {
87                target
88            }
89            #[inline]
90            fn modulus(&self) -> $T {
91                Self::MODULUS
92            }
93            #[inline]
94            fn is_zero(&self, target: &$T) -> bool {
95                target == &0
96            }
97
98            #[inline]
99            fn add(&self, lhs: &$T, rhs: &$T) -> $T {
100                let (sum, overflow) = lhs.overflowing_add(*rhs);
101                if overflow || sum >= Self::MODULUS {
102                    let (sum2, _) = sum.overflowing_sub(Self::MODULUS);
103                    sum2
104                } else {
105                    sum
106                }
107            }
108            #[inline]
109            fn sub(&self, lhs: &$T, rhs: &$T) -> $T {
110                if lhs >= rhs {
111                    lhs - rhs
112                } else {
113                    Self::MODULUS - (rhs - lhs)
114                }
115            }
116            #[inline]
117            fn dbl(&self, target: $T) -> $T {
118                let (sum, overflow) = target.overflowing_add(target);
119                if overflow || sum >= Self::MODULUS {
120                    let (sum2, _) = sum.overflowing_sub(Self::MODULUS);
121                    sum2
122                } else {
123                    sum
124                }
125            }
126            #[inline]
127            fn neg(&self, target: $T) -> $T {
128                if target == 0 {
129                    0
130                } else {
131                    Self::MODULUS - target
132                }
133            }
134            #[inline]
135            fn mul(&self, lhs: &$T, rhs: &$T) -> $T {
136                if (P1 as u32) < $half_bits {
137                    Self::reduce_single(lhs * rhs)
138                } else {
139                    Self::reduce_double(impl_fixed_trinomial_solinas!(@widen_mul, $kind, $T, $D, lhs, rhs))
140                }
141            }
142            #[inline]
143            fn inv(&self, target: $T) -> Option<$T> {
144                // TODO: inv can be specialized
145                // REF: https://xn--2-umb.com/22/goldilocks/
146                if (P1 as u32) < usize::BITS {
147                    (target as usize)
148                        .invm(&(Self::MODULUS as usize))
149                        .map(|v| v as $T)
150                } else {
151                    target.invm(&Self::MODULUS)
152                }
153            }
154            #[inline]
155            fn sqr(&self, target: $T) -> $T {
156                if (P1 as u32) < $half_bits {
157                    Self::reduce_single(target * target)
158                } else {
159                    Self::reduce_double(impl_fixed_trinomial_solinas!(@widen_sqr, $kind, $T, $D, target))
160                }
161            }
162
163            impl_reduced_binary_pow!($T);
164        }
165    };
166
167    // Internal: reduce_single for primitive double-width types (u32→u64, u64→u128)
168    (@reduce_single, primitive, $T:ty, $D:ty) => {
169        const fn reduce_single(v: $T) -> $T {
170            let mut v: $D = v as $D;
171            while v >> P1 > 0 {
172                let lo = (v as $T) & Self::BITMASK;
173                let hi = v >> P1;
174                let mut sum: $D = (hi << (P2 as u32)) + (lo as $D);
175                if K > 0 {
176                    sum -= hi * (K as $D);
177                } else if K < 0 {
178                    sum += hi * ((-K) as $D);
179                }
180                v = sum;
181            }
182            let v = v as $T;
183            if v >= Self::MODULUS {
184                v - Self::MODULUS
185            } else {
186                v
187            }
188        }
189    };
190
191    // Internal: reduce_single for udouble (umax→udouble). Stays in udouble for the same reason
192    // as reduce_double below: `hi << P2` can exceed `umax` during the tail.
193    (@reduce_single, udouble, $T:ty, $D:ty) => {
194        fn reduce_single(v: $T) -> $T {
195            let mut v: $D = udouble { hi: 0, lo: v };
196            while v.hi > 0 || v.lo >> P1 > 0 {
197                let lo = v.lo & Self::BITMASK;
198                let hi = v >> P1;
199                let mut sum = (hi << (P2 as u32)) + lo;
200                if K > 0 {
201                    sum -= hi * (K as umax);
202                } else if K < 0 {
203                    sum += hi * ((-K) as umax);
204                }
205                v = sum;
206            }
207            let v = v.lo;
208            if v >= Self::MODULUS {
209                v - Self::MODULUS
210            } else {
211                v
212            }
213        }
214    };
215
216    // Internal: reduce_double for primitive double-width types (u32→u64, u64→u128)
217    //
218    // When the worst-case fold count is small, replace the while loop with
219    // straight-line unconditional folds. Each fold is a no-op once hi reaches 0.
220    // FOLDS from the expert formula: ⌈P1/(P1−P2)⌉ + 1 (K>0) or +2 (K<0).
221    // Unrolling condition: P2 ≤ ⌊2·P1/3⌋  ⇔  FOLDS ≤ 4.
222    (@reduce_double, primitive, $T:ty, $D:ty) => {
223        fn reduce_double(v: $D) -> $T {
224            let mut lo = (v as $T) & Self::BITMASK;
225            let mut hi = v >> P1;
226            macro_rules! solinas_fold {
227                () => {
228                    let mut sum: $D = (hi << (P2 as u32)) + (lo as $D);
229                    if K > 0 { sum -= hi * (K as $D); }
230                    else if K < 0 { sum += hi * ((-K) as $D); }
231                    lo = (sum as $T) & Self::BITMASK;
232                    hi = sum >> P1;
233                };
234            }
235            if Self::FOLDS <= 3 {
236                #[allow(unused_assignments)] { solinas_fold!(); }
237                #[allow(unused_assignments)] { solinas_fold!(); }
238                #[allow(unused_assignments)] { solinas_fold!(); }
239            } else if Self::FOLDS == 4 {
240                #[allow(unused_assignments)] { solinas_fold!(); }
241                #[allow(unused_assignments)] { solinas_fold!(); }
242                #[allow(unused_assignments)] { solinas_fold!(); }
243                #[allow(unused_assignments)] { solinas_fold!(); }
244            } else {
245                while hi > 0 { solinas_fold!(); }
246            }
247            if lo >= Self::MODULUS {
248                lo - Self::MODULUS
249            } else {
250                lo
251            }
252        }
253    };
254
255    // Internal: reduce_double for udouble (u128→udouble)
256    //
257    // Unlike [Mersenne](crate::FixedMersenne)'s two-phase loop (udouble while `hi.hi > 0`, then
258    // `umax` while `hi.lo > 0`), Solinas keeps `hi` as [udouble] until fully zero. Mersenne's
259    // tail step is `hi * K + lo`, which stays within `umax` when `K < 2^(P-1)`. Solinas uses
260    // `hi << P2`, which can exceed `umax` even when `hi` fits in one word (e.g. `hi * 2^P2`), so
261    // the tail must stay in double-width arithmetic.
262    (@reduce_double, udouble, $T:ty, $D:ty) => {
263        fn reduce_double(v: $D) -> $T {
264            let mut lo = v.lo & Self::BITMASK;
265            let mut hi = v >> P1;
266            macro_rules! udouble_fold {
267                () => {
268                    let mut sum = (hi << (P2 as u32)) + lo;
269                    if K > 0 { sum -= hi * (K as umax); }
270                    else if K < 0 { sum += hi * ((-K) as umax); }
271                    lo = sum.lo & Self::BITMASK;
272                    hi = sum >> P1;
273                };
274            }
275            if Self::FOLDS <= 3 {
276                #[allow(unused_assignments)] { udouble_fold!(); }
277                #[allow(unused_assignments)] { udouble_fold!(); }
278                #[allow(unused_assignments)] { udouble_fold!(); }
279            } else if Self::FOLDS == 4 {
280                #[allow(unused_assignments)] { udouble_fold!(); }
281                #[allow(unused_assignments)] { udouble_fold!(); }
282                #[allow(unused_assignments)] { udouble_fold!(); }
283                #[allow(unused_assignments)] { udouble_fold!(); }
284            } else {
285                while hi.hi > 0 || hi.lo > 0 { udouble_fold!(); }
286            }
287            if lo >= Self::MODULUS {
288                lo - Self::MODULUS
289            } else {
290                lo
291            }
292        }
293    };
294
295    // Internal: widening multiplication for primitive types
296    (@widen_mul, primitive, $T:ty, $D:ty, $lhs:expr, $rhs:expr) => {
297        (*$lhs as $D) * (*$rhs as $D)
298    };
299
300    // Internal: widening multiplication for udouble
301    (@widen_mul, udouble, $T:ty, $D:ty, $lhs:expr, $rhs:expr) => {
302        <$D>::widening_mul(*$lhs, *$rhs)
303    };
304
305    // Internal: widening square for primitive types
306    (@widen_sqr, primitive, $T:ty, $D:ty, $target:expr) => {
307        ($target as $D) * ($target as $D)
308    };
309
310    // Internal: widening square for udouble
311    (@widen_sqr, udouble, $T:ty, $D:ty, $target:expr) => {
312        <$D>::widening_square($target)
313    };
314}
315
316/// A modular reducer for trinomial Solinas numbers `2^P1 - 2^P2 + K` as modulus with 32-bit operands.
317///
318/// Supports `P1` up to 31, `P2 < P1`, and odd signed `K` with `|K| < 2^P2`. All inputs and outputs are `u32`.
319/// The modulus `2^P1 - 2^P2 + K` must be prime for modular inverse and Fermat-based operations to be valid.
320///
321/// # Example
322///
323/// ```rust
324/// use num_modular::{FixedTrinomialSolinas32, Reducer};
325///
326/// const P1: u8 = 4;
327/// const P2: u8 = 2;
328/// const K: i32 = 1;
329/// let modulus = (1u32 << P1) - (1u32 << P2) + (K as u32); // 2^4 - 2^2 + 1 = 13
330/// let reducer = FixedTrinomialSolinas32::<P1, P2, K>::new(&modulus);
331/// let a = reducer.transform(3);
332/// let b = reducer.transform(5);
333/// assert_eq!(reducer.residue(reducer.add(&a, &b)), 8);
334/// ```
335#[derive(Debug, Clone, Copy)]
336pub struct FixedTrinomialSolinas32<const P1: u8, const P2: u8, const K: i32>();
337
338impl_fixed_trinomial_solinas!(FixedTrinomialSolinas32, u32, i32, u64, 16, 31, primitive);
339
340/// A modular reducer for trinomial Solinas numbers `2^P1 - 2^P2 + K` as modulus with 64-bit operands.
341///
342/// Supports `P1` up to 64, `P2 < P1`, and odd signed `K` with `|K| < 2^P2`. All inputs and outputs are `u64`.
343/// Uses `u128` as the double-width intermediate for multiplication and reduction.
344/// The modulus `2^P1 - 2^P2 + K` must be prime for modular inverse and Fermat-based operations to be valid.
345///
346/// # Example
347///
348/// ```rust
349/// use num_modular::{FixedTrinomialSolinas64, Reducer};
350///
351/// const P1: u8 = 6;
352/// const P2: u8 = 2;
353/// const K: i64 = 1;
354/// let modulus = (1u64 << P1) - (1u64 << P2) + (K as u64); // 2^6 - 2^2 + 1 = 61
355/// let reducer = FixedTrinomialSolinas64::<P1, P2, K>::new(&modulus);
356/// let a = reducer.transform(10);
357/// let b = reducer.transform(20);
358/// assert_eq!(reducer.residue(reducer.mul(&a, &b)), (10u64 * 20) % 61);
359/// ```
360#[derive(Debug, Clone, Copy)]
361pub struct FixedTrinomialSolinas64<const P1: u8, const P2: u8, const K: i64>();
362
363impl_fixed_trinomial_solinas!(FixedTrinomialSolinas64, u64, i64, u128, 32, 64, primitive);
364
365/// A modular reducer for trinomial Solinas numbers `2^P1 - 2^P2 + K` as modulus.
366///
367/// Supports `P1` up to 127, `P2 < P1`, and odd signed `K` with `|K| < 2^P2`. All inputs and outputs are [umax] (currently `u128`).
368///
369/// The `P1` is limited to 127 so that overflow checks aren't necessary. This covers all trinomial
370/// Solinas primes within the range of [umax] (i.e. `u128`).
371///
372/// # Example
373///
374/// ```rust
375/// use num_modular::{FixedTrinomialSolinas, Reducer};
376///
377/// const P1: u8 = 31;
378/// const P2: u8 = 13;
379/// const K: i128 = 1;
380/// let modulus = (1u128 << P1) - (1u128 << P2) + (K as u128);
381/// let reducer = FixedTrinomialSolinas::<P1, P2, K>::new(&modulus);
382/// let a = reducer.transform(1000);
383/// let b = reducer.transform(2000);
384/// assert_eq!(reducer.residue(reducer.mul(&a, &b)), (1000u128 * 2000) % modulus);
385/// ```
386#[derive(Debug, Clone, Copy)]
387pub struct FixedTrinomialSolinas<const P1: u8, const P2: u8, const K: imax>();
388
389impl_fixed_trinomial_solinas!(FixedTrinomialSolinas, umax, imax, udouble, 64, 127, udouble);
390
391#[cfg(test)]
392mod tests {
393    use super::*;
394    use crate::{ModularCoreOps, ModularPow};
395    use rand::random;
396
397    // u128 types
398    type S1 = FixedTrinomialSolinas<31, 13, 1>;
399    type S2 = FixedTrinomialSolinas<61, 30, 1>;
400    type S3 = FixedTrinomialSolinas<127, 64, 1>;
401    type S4 = FixedTrinomialSolinas<32, 16, 1>;
402    type S5 = FixedTrinomialSolinas<56, 28, -1>;
403    type S6 = FixedTrinomialSolinas<122, 61, -3>;
404
405    // u64 types
406    type S64_1 = FixedTrinomialSolinas64<31, 13, 1>;
407    type S64_2 = FixedTrinomialSolinas64<61, 30, 1>;
408    type S64_3 = FixedTrinomialSolinas64<32, 16, 1>;
409    type S64_4 = FixedTrinomialSolinas64<64, 32, 1>; // 2^64 - 2^32 + 1
410
411    // u32 types
412    type S32_1 = FixedTrinomialSolinas32<4, 2, 1>;
413    type S32_2 = FixedTrinomialSolinas32<5, 3, -1>;
414    type S32_3 = FixedTrinomialSolinas32<6, 2, 1>;
415
416    const NRANDOM: u32 = 10;
417
418    #[test]
419    fn creation_test_u128() {
420        const P: umax = <S1>::MODULUS;
421        let m = S1::new(&P);
422        assert_eq!(m.residue(m.transform(0)), 0);
423        assert_eq!(m.residue(m.transform(1)), 1);
424        assert_eq!(m.residue(m.transform(P)), 0);
425        assert_eq!(m.residue(m.transform(P - 1)), P - 1);
426        assert_eq!(m.residue(m.transform(P + 1)), 1);
427
428        for _ in 0..NRANDOM {
429            let a = random::<umax>();
430
431            const P1: umax = <S1>::MODULUS;
432            let m1 = S1::new(&P1);
433            assert_eq!(m1.residue(m1.transform(a)), a % P1);
434            const P2: umax = <S2>::MODULUS;
435            let m2 = S2::new(&P2);
436            assert_eq!(m2.residue(m2.transform(a)), a % P2);
437            const P3: umax = <S3>::MODULUS;
438            let m3 = S3::new(&P3);
439            assert_eq!(m3.residue(m3.transform(a)), a % P3);
440            const P4: umax = <S4>::MODULUS;
441            let m4 = S4::new(&P4);
442            assert_eq!(m4.residue(m4.transform(a)), a % P4);
443            const P5: umax = <S5>::MODULUS;
444            let m5 = S5::new(&P5);
445            assert_eq!(m5.residue(m5.transform(a)), a % P5);
446            const P6: umax = <S6>::MODULUS;
447            let m6 = S6::new(&P6);
448            assert_eq!(m6.residue(m6.transform(a)), a % P6);
449        }
450    }
451
452    #[test]
453    fn creation_test_u64() {
454        for _ in 0..NRANDOM {
455            let a = random::<u64>();
456
457            const P1: u64 = <S64_1>::MODULUS;
458            let m1 = S64_1::new(&P1);
459            assert_eq!(m1.residue(m1.transform(a)), a % P1);
460            const P2: u64 = <S64_2>::MODULUS;
461            let m2 = S64_2::new(&P2);
462            assert_eq!(m2.residue(m2.transform(a)), a % P2);
463            const P3: u64 = <S64_3>::MODULUS;
464            let m3 = S64_3::new(&P3);
465            assert_eq!(m3.residue(m3.transform(a)), a % P3);
466            const P4: u64 = <S64_4>::MODULUS;
467            let m4 = S64_4::new(&P4);
468            assert_eq!(m4.residue(m4.transform(a)), a % P4);
469        }
470    }
471
472    #[test]
473    fn creation_test_u32() {
474        for _ in 0..NRANDOM {
475            let a = random::<u32>();
476
477            const P1: u32 = <S32_1>::MODULUS;
478            let m1 = S32_1::new(&P1);
479            assert_eq!(m1.residue(m1.transform(a)), a % P1);
480            const P2: u32 = <S32_2>::MODULUS;
481            let m2 = S32_2::new(&P2);
482            assert_eq!(m2.residue(m2.transform(a)), a % P2);
483            const P3: u32 = <S32_3>::MODULUS;
484            let m3 = S32_3::new(&P3);
485            assert_eq!(m3.residue(m3.transform(a)), a % P3);
486        }
487    }
488
489    #[test]
490    fn test_against_modops_u128() {
491        macro_rules! tests_for {
492            ($a:tt, $b:tt, $e:tt; $($M:ty)*) => ($({
493                const P: umax = <$M>::MODULUS;
494                let r = <$M>::new(&P);
495                let am = r.transform($a);
496                let bm = r.transform($b);
497                assert_eq!(r.add(&am, &bm), $a.addm($b, &P));
498                assert_eq!(r.sub(&am, &bm), $a.subm($b, &P));
499                assert_eq!(r.mul(&am, &bm), $a.mulm($b, &P));
500                assert_eq!(r.neg(am), $a.negm(&P));
501                assert_eq!(r.inv(am), $a.invm(&P));
502                assert_eq!(r.dbl(am), $a.dblm(&P));
503                assert_eq!(r.sqr(am), $a.sqm(&P));
504                assert_eq!(r.pow(am, &$e), $a.powm($e, &P));
505            })*);
506        }
507
508        for _ in 0..NRANDOM {
509            let (a, b) = (random::<u128>(), random::<u128>());
510            let e = random::<u8>() as umax;
511            tests_for!(a, b, e; S1 S2 S3 S4 S5 S6);
512        }
513    }
514
515    #[test]
516    fn test_against_modops_u64() {
517        macro_rules! tests_for {
518            ($a:ident, $b:ident, $e:ident; $($M:ty)*) => ($({
519                const P: u64 = <$M>::MODULUS;
520                let r = <$M>::new(&P);
521                let am = r.transform($a);
522                let bm = r.transform($b);
523                assert_eq!(r.add(&am, &bm), $a.addm($b, &P));
524                assert_eq!(r.sub(&am, &bm), $a.subm($b, &P));
525                assert_eq!(r.mul(&am, &bm), $a.mulm($b, &P));
526                assert_eq!(r.neg(am), $a.negm(&P));
527                assert_eq!(r.inv(am), $a.invm(&P));
528                assert_eq!(r.dbl(am), $a.dblm(&P));
529                assert_eq!(r.sqr(am), $a.sqm(&P));
530                assert_eq!(r.pow(am, &$e), $a.powm($e, &P));
531            })*);
532        }
533
534        for _ in 0..NRANDOM {
535            let a = random::<u64>();
536            let b = random::<u64>();
537            let e = random::<u8>() as u64;
538            tests_for!(a, b, e; S64_1 S64_2 S64_3 S64_4);
539        }
540    }
541
542    #[test]
543    fn test_against_modops_u32() {
544        macro_rules! tests_for {
545            ($a:ident, $b:ident, $e:ident; $($M:ty)*) => ($({
546                const P: u32 = <$M>::MODULUS;
547                let r = <$M>::new(&P);
548                let am = r.transform($a);
549                let bm = r.transform($b);
550                assert_eq!(r.add(&am, &bm), $a.addm($b, &P));
551                assert_eq!(r.sub(&am, &bm), $a.subm($b, &P));
552                assert_eq!(r.mul(&am, &bm), $a.mulm($b, &P));
553                assert_eq!(r.neg(am), $a.negm(&P));
554                assert_eq!(r.inv(am), $a.invm(&P));
555                assert_eq!(r.dbl(am), $a.dblm(&P));
556                assert_eq!(r.sqr(am), $a.sqm(&P));
557                assert_eq!(r.pow(am, &$e), $a.powm($e, &P));
558            })*);
559        }
560
561        for _ in 0..NRANDOM {
562            let a = random::<u32>();
563            let b = random::<u32>();
564            let e = random::<u8>() as u32;
565            tests_for!(a, b, e; S32_1 S32_2 S32_3);
566        }
567    }
568
569    #[test]
570    fn test_add_near_overflow_u64() {
571        // 2^64 - 2^32 + 1 = 0xFFFFFFFF00000001, near u64::MAX
572        type S = FixedTrinomialSolinas64<64, 32, 1>;
573        const P: u64 = <S>::MODULUS;
574        assert_eq!(P, 0xFFFFFFFF00000001);
575        let r = S::new(&P);
576        // Values near P-1; their sum exceeds u64::MAX
577        // (P-1) + (P-2) = 2P-3 ≡ P-3 (mod P)
578        let a = r.transform(P - 1);
579        let b = r.transform(P - 2);
580        assert_eq!(r.residue(r.add(&a, &b)), P - 3);
581        // dbl near overflow: 2*(P-1) = 2P-2 ≡ P-2 (mod P)
582        let c = r.transform(P - 1);
583        assert_eq!(r.residue(r.dbl(c)), P - 2);
584    }
585}
num_modular/solinas.rs

num_modular/
solinas.rs
