num_modular/
lib.rs

1//! This crate provides efficient Modular arithmetic operations for various integer types,
2//! including primitive integers and `num-bigint`. The latter option is enabled optionally.
3//!
4//! To achieve fast modular arithmetics, convert integers to any [ModularInteger] implementation
5//! use static `new()` or associated [ModularInteger::convert()] functions. Some builtin implementations
6//! of [ModularInteger] includes [MontgomeryInt] and [FixedMersenneInt].
7//!
8//! Example code:
9//! ```rust
10//! use num_modular::{ModularCoreOps, ModularInteger, MontgomeryInt};
11//!
12//! // directly using methods in ModularCoreOps
13//! let (x, y, m) = (12u8, 13u8, 5u8);
14//! assert_eq!(x.mulm(y, &m), x * y % m);
15//!
16//! // convert integers into ModularInteger
17//! let mx = MontgomeryInt::new(x, &m);
18//! let my = mx.convert(y); // faster than static MontgomeryInt::new(y, m)
19//! assert_eq!((mx * my).residue(), x * y % m);
20//! ```
21//!
22//! # Comparison of fast division / modular arithmetics
23//! Several fast division / modulo tricks are provided in these crate, the difference of them are listed below:
24//! - [PreModInv]: pre-compute modular inverse of the divisor, only applicable to exact division
25//! - Barret (to be implemented): pre-compute (rational approximation of) the reciprocal of the divisor,
26//!     applicable to fast division and modulo
27//! - [Montgomery]: Convert the dividend into a special form by shifting and pre-compute a modular inverse,
28//!     only applicable to fast modulo, but faster than Barret reduction
29//! - [FixedMersenne]: Specialization of modulo in form `2^P-K` under 2^127.
30//!
31
32// XXX: Other fast modular arithmetic tricks
33// REF: https://github.com/lemire/fastmod & https://arxiv.org/pdf/1902.01961.pdf
34// REF: https://eprint.iacr.org/2014/040.pdf
35// REF: https://github.com/ridiculousfish/libdivide/
36// REF: Faster Interleaved Modular Multiplication Based on Barrett and Montgomery Reduction Methods (work for modulus in certain form)
37
38#![no_std]
39#[cfg(any(feature = "std", test))]
40extern crate std;
41
42use core::ops::{Add, Mul, Neg, Sub};
43
44/// Core modular arithmetic operations.
45///
46/// Note that all functions will panic if the modulus is zero.
47pub trait ModularCoreOps<Rhs = Self, Modulus = Self> {
48    type Output;
49
50    /// Return (self + rhs) % m
51    fn addm(self, rhs: Rhs, m: Modulus) -> Self::Output;
52
53    /// Return (self - rhs) % m
54    fn subm(self, rhs: Rhs, m: Modulus) -> Self::Output;
55
56    /// Return (self * rhs) % m
57    fn mulm(self, rhs: Rhs, m: Modulus) -> Self::Output;
58}
59
60/// Core unary modular arithmetics
61///
62/// Note that all functions will panic if the modulus is zero.
63pub trait ModularUnaryOps<Modulus = Self> {
64    type Output;
65
66    /// Return (-self) % m and make sure the result is normalized in range [0,m)
67    fn negm(self, m: Modulus) -> Self::Output;
68
69    /// Calculate modular inverse (x such that self*x = 1 mod m).
70    ///
71    /// This operation is only available for integer that is coprime to `m`. If not,
72    /// the result will be [None].
73    fn invm(self, m: Modulus) -> Option<Self::Output>;
74
75    /// Calculate modular double ( x+x mod m)
76    fn dblm(self, m: Modulus) -> Self::Output;
77
78    /// Calculate modular square ( x*x mod m )
79    fn sqm(self, m: Modulus) -> Self::Output;
80
81    // TODO: Modular sqrt aka Quadratic residue, follow the behavior of FLINT `n_sqrtmod`
82    // fn sqrtm(self, m: Modulus) -> Option<Self::Output>;
83    // REF: https://stackoverflow.com/questions/6752374/cube-root-modulo-p-how-do-i-do-this
84}
85
86/// Modular power functions
87pub trait ModularPow<Exp = Self, Modulus = Self> {
88    type Output;
89
90    /// Return (self ^ exp) % m
91    fn powm(self, exp: Exp, m: Modulus) -> Self::Output;
92}
93
94/// Math symbols related to modular arithmetics
95pub trait ModularSymbols<Modulus = Self> {
96    /// Calculate Legendre Symbol (a|n), where a is `self`.
97    ///
98    /// Note that this function doesn't perform a full primality check, since
99    /// is costly. So if n is not a prime, the result can be not reasonable.
100    ///
101    /// # Panics
102    /// Only if n is not prime
103    #[inline]
104    fn legendre(&self, n: Modulus) -> i8 {
105        self.checked_legendre(n).expect("n shoud be a prime")
106    }
107
108    /// Calculate Legendre Symbol (a|n), where a is `self`. Returns [None] only if n is
109    /// not a prime.
110    ///
111    /// Note that this function doesn't perform a full primality check, since
112    /// is costly. So if n is not a prime, the result can be not reasonable.
113    ///
114    /// # Panics
115    /// Only if n is not prime
116    fn checked_legendre(&self, n: Modulus) -> Option<i8>;
117
118    /// Calculate Jacobi Symbol (a|n), where a is `self`
119    ///
120    /// # Panics
121    /// if n is negative or even
122    #[inline]
123    fn jacobi(&self, n: Modulus) -> i8 {
124        self.checked_jacobi(n)
125            .expect("the Jacobi symbol is only defined for non-negative odd integers")
126    }
127
128    /// Calculate Jacobi Symbol (a|n), where a is `self`. Returns [None] if n is negative or even.
129    fn checked_jacobi(&self, n: Modulus) -> Option<i8>;
130
131    /// Calculate Kronecker Symbol (a|n), where a is `self`
132    fn kronecker(&self, n: Modulus) -> i8;
133}
134
135// TODO: Discrete log aka index, follow the behavior of FLINT `n_discrete_log_bsgs`
136// REF: https://github.com/vks/discrete-log
137// fn logm(self, base: Modulus, m: Modulus);
138
139/// Collection of common modular arithmetic operations
140pub trait ModularOps<Rhs = Self, Modulus = Self, Output = Self>:
141    ModularCoreOps<Rhs, Modulus, Output = Output>
142    + ModularUnaryOps<Modulus, Output = Output>
143    + ModularPow<Rhs, Modulus, Output = Output>
144    + ModularSymbols<Modulus>
145{
146}
147impl<T, Rhs, Modulus> ModularOps<Rhs, Modulus> for T where
148    T: ModularCoreOps<Rhs, Modulus, Output = T>
149        + ModularUnaryOps<Modulus, Output = T>
150        + ModularPow<Rhs, Modulus, Output = T>
151        + ModularSymbols<Modulus>
152{
153}
154
155/// Collection of operations similar to [ModularOps], but takes operands with references
156pub trait ModularRefOps: for<'r> ModularOps<&'r Self, &'r Self> + Sized {}
157impl<T> ModularRefOps for T where T: for<'r> ModularOps<&'r T, &'r T> {}
158
159/// Provides a utility function to convert signed integers into unsigned modular form
160pub trait ModularAbs<Modulus> {
161    /// Return self % m, but accepting signed integers
162    fn absm(self, m: &Modulus) -> Modulus;
163}
164
165/// Represents an number defined in a modulo ring ℤ/nℤ
166///
167/// The operators should panic if the modulus of two number
168/// are not the same.
169pub trait ModularInteger:
170    Sized
171    + PartialEq
172    + Add<Self, Output = Self>
173    + Sub<Self, Output = Self>
174    + Neg<Output = Self>
175    + Mul<Self, Output = Self>
176{
177    /// The underlying representation type of the integer
178    type Base;
179
180    /// Return the modulus of the ring
181    fn modulus(&self) -> Self::Base;
182
183    /// Return the normalized residue of this integer in the ring
184    fn residue(&self) -> Self::Base;
185
186    /// Check if the integer is zero
187    fn is_zero(&self) -> bool;
188
189    /// Convert an normal integer into the same ring.
190    ///
191    /// This method should be perferred over the static
192    /// constructor to prevent unnecessary overhead of pre-computation.
193    fn convert(&self, n: Self::Base) -> Self;
194
195    /// Calculate the value of self + self
196    fn double(self) -> Self;
197
198    /// Calculate the value of self * self
199    fn square(self) -> Self;
200}
201
202// XXX: implement ModularInteger for ff::PrimeField?
203// TODO: implement invm_range (Modular inverse in certain range) and crt (Chinese Remainder Theorem), REF: bubblemath crate
204
205/// Utility function for exact division, with precomputed helper values
206///
207/// # Available Pre-computation types:
208/// - `()`: No pre-computation, the implementation relies on native integer division
209/// - [PreModInv]: With Pre-computed modular inverse
210pub trait DivExact<Rhs, Precompute>: Sized {
211    type Output;
212
213    /// Check if d divides self with the help of the precomputation. If d divides self,
214    /// then the quotient is returned.
215    fn div_exact(self, d: Rhs, pre: &Precompute) -> Option<Self::Output>;
216}
217
218/// A modular reducer that can ensure that the operations on integers are all performed
219/// in a modular ring.
220///
221/// Essential information for performing the modulo operation will be stored in the reducer.
222pub trait Reducer<T> {
223    /// Create a reducer for a modulus m
224    fn new(m: &T) -> Self;
225
226    /// Transform a normal integer into reduced form
227    fn transform(&self, target: T) -> T;
228
229    /// Check whether target is a valid reduced form
230    fn check(&self, target: &T) -> bool;
231
232    /// Get the modulus in original integer type
233    fn modulus(&self) -> T;
234
235    /// Transform a reduced form back to normal integer
236    fn residue(&self, target: T) -> T;
237
238    /// Test if the residue() == 0
239    fn is_zero(&self, target: &T) -> bool;
240
241    /// Calculate (lhs + rhs) mod m in reduced form
242    fn add(&self, lhs: &T, rhs: &T) -> T;
243
244    #[inline]
245    fn add_in_place(&self, lhs: &mut T, rhs: &T) {
246        *lhs = self.add(lhs, rhs)
247    }
248
249    /// Calculate 2*target mod m
250    fn dbl(&self, target: T) -> T;
251
252    /// Calculate (lhs - rhs) mod m in reduced form
253    fn sub(&self, lhs: &T, rhs: &T) -> T;
254
255    #[inline]
256    fn sub_in_place(&self, lhs: &mut T, rhs: &T) {
257        *lhs = self.sub(lhs, rhs);
258    }
259
260    /// Calculate -monty mod m in reduced form
261    fn neg(&self, target: T) -> T;
262
263    /// Calculate (lhs * rhs) mod m in reduced form
264    fn mul(&self, lhs: &T, rhs: &T) -> T;
265
266    #[inline]
267    fn mul_in_place(&self, lhs: &mut T, rhs: &T) {
268        *lhs = self.mul(lhs, rhs);
269    }
270
271    /// Calculate target^-1 mod m in reduced form,
272    /// it may return None when there is no modular inverse.
273    fn inv(&self, target: T) -> Option<T>;
274
275    /// Calculate target^2 mod m in reduced form
276    fn sqr(&self, target: T) -> T;
277
278    /// Calculate base ^ exp mod m in reduced form
279    fn pow(&self, base: T, exp: &T) -> T;
280}
281
282mod barret;
283mod double;
284mod mersenne;
285mod monty;
286mod preinv;
287mod prim;
288mod reduced;
289mod word;
290
291pub use barret::{
292    Normalized2by1Divisor, Normalized3by2Divisor, PreMulInv1by1, PreMulInv2by1, PreMulInv3by2,
293};
294pub use double::{udouble, umax};
295pub use mersenne::FixedMersenne;
296pub use monty::Montgomery;
297pub use preinv::PreModInv;
298pub use reduced::{ReducedInt, Vanilla, VanillaInt};
299
300/// An integer in modulo ring based on [Montgomery form](https://en.wikipedia.org/wiki/Montgomery_modular_multiplication#Montgomery_form)
301pub type MontgomeryInt<T> = ReducedInt<T, Montgomery<T>>;
302
303/// An integer in modulo ring with a fixed (pseudo) Mersenne number as modulus
304pub type FixedMersenneInt<const P: u8, const K: umax> = ReducedInt<umax, FixedMersenne<P, K>>;
305
306// pub type BarretInt<T> = ReducedInt<T, BarretInt<T>>;
307
308#[cfg(feature = "num-bigint")]
309mod bigint;