1
2
3
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
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
//! This crate provides efficient Modular arithmetic operations for various integer types,
//! including primitive integers and `num-bigint`. The latter option is enabled optionally.
//!
//! To achieve fast modular arithmetics, convert integers to any [ModularInteger] implementation
//! use static `new()` or associated [ModularInteger::convert()] functions. Some builtin implementations
//! of [ModularInteger] includes [MontgomeryInt] and [FixedMersenneInt].
//!
//! Example code:
//! ```rust
//! use num_modular::{ModularCoreOps, ModularInteger, MontgomeryInt};
//!
//! // directly using methods in ModularCoreOps
//! let (x, y, m) = (12u8, 13u8, 5u8);
//! assert_eq!(x.mulm(y, &m), x * y % m);
//!
//! // convert integers into ModularInteger
//! let mx = MontgomeryInt::new(x, &m);
//! let my = mx.convert(y); // faster than static MontgomeryInt::new(y, m)
//! assert_eq!((mx * my).residue(), x * y % m);
//! ```
//!
//! # Comparison of fast division / modular arithmetics
//! Several fast division / modulo tricks are provided in these crate, the difference of them are listed below:
//! - [PreModInv]: pre-compute modular inverse of the divisor, only applicable to exact division
//! - Barret (to be implemented): pre-compute (rational approximation of) the reciprocal of the divisor,
//!     applicable to fast division and modulo
//! - [Montgomery]: Convert the dividend into a special form by shifting and pre-compute a modular inverse,
//!     only applicable to fast modulo, but faster than Barret reduction
//! - [FixedMersenne]: Specialization of modulo in form `2^P-K` under 2^127.
//!

// XXX: Other fast modular arithmetic tricks
// REF: https://github.com/lemire/fastmod & https://arxiv.org/pdf/1902.01961.pdf
// REF: https://eprint.iacr.org/2014/040.pdf
// REF: https://github.com/ridiculousfish/libdivide/
// REF: Faster Interleaved Modular Multiplication Based on Barrett and Montgomery Reduction Methods (work for modulus in certain form)

#![no_std]
#[cfg(any(feature = "std", test))]
extern crate std;

use core::ops::{Add, Mul, Neg, Sub};

/// Core modular arithmetic operations.
///
/// Note that all functions will panic if the modulus is zero.
pub trait ModularCoreOps<Rhs = Self, Modulus = Self> {
    type Output;

    /// Return (self + rhs) % m
    fn addm(self, rhs: Rhs, m: Modulus) -> Self::Output;

    /// Return (self - rhs) % m
    fn subm(self, rhs: Rhs, m: Modulus) -> Self::Output;

    /// Return (self * rhs) % m
    fn mulm(self, rhs: Rhs, m: Modulus) -> Self::Output;
}

/// Core unary modular arithmetics
///
/// Note that all functions will panic if the modulus is zero.
pub trait ModularUnaryOps<Modulus = Self> {
    type Output;

    /// Return (-self) % m and make sure the result is normalized in range [0,m)
    fn negm(self, m: Modulus) -> Self::Output;

    /// Calculate modular inverse (x such that self*x = 1 mod m).
    ///
    /// This operation is only available for integer that is coprime to `m`. If not,
    /// the result will be [None].
    fn invm(self, m: Modulus) -> Option<Self::Output>;

    /// Calculate modular double ( x+x mod m)
    fn dblm(self, m: Modulus) -> Self::Output;

    /// Calculate modular square ( x*x mod m )
    fn sqm(self, m: Modulus) -> Self::Output;

    // TODO: Modular sqrt aka Quadratic residue, follow the behavior of FLINT `n_sqrtmod`
    // fn sqrtm(self, m: Modulus) -> Option<Self::Output>;
    // REF: https://stackoverflow.com/questions/6752374/cube-root-modulo-p-how-do-i-do-this
}

/// Modular power functions
pub trait ModularPow<Exp = Self, Modulus = Self> {
    type Output;

    /// Return (self ^ exp) % m
    fn powm(self, exp: Exp, m: Modulus) -> Self::Output;
}

/// Math symbols related to modular arithmetics
pub trait ModularSymbols<Modulus = Self> {
    /// Calculate Legendre Symbol (a|n), where a is `self`.
    ///
    /// Note that this function doesn't perform a full primality check, since
    /// is costly. So if n is not a prime, the result can be not reasonable.
    ///
    /// # Panics
    /// Only if n is not prime
    #[inline]
    fn legendre(&self, n: Modulus) -> i8 {
        self.checked_legendre(n).expect("n shoud be a prime")
    }

    /// Calculate Legendre Symbol (a|n), where a is `self`. Returns [None] only if n is
    /// not a prime.
    ///
    /// Note that this function doesn't perform a full primality check, since
    /// is costly. So if n is not a prime, the result can be not reasonable.
    ///
    /// # Panics
    /// Only if n is not prime
    fn checked_legendre(&self, n: Modulus) -> Option<i8>;

    /// Calculate Jacobi Symbol (a|n), where a is `self`
    ///
    /// # Panics
    /// if n is negative or even
    #[inline]
    fn jacobi(&self, n: Modulus) -> i8 {
        self.checked_jacobi(n)
            .expect("the Jacobi symbol is only defined for non-negative odd integers")
    }

    /// Calculate Jacobi Symbol (a|n), where a is `self`. Returns [None] if n is negative or even.
    fn checked_jacobi(&self, n: Modulus) -> Option<i8>;

    /// Calculate Kronecker Symbol (a|n), where a is `self`
    fn kronecker(&self, n: Modulus) -> i8;
}

// TODO: Discrete log aka index, follow the behavior of FLINT `n_discrete_log_bsgs`
// REF: https://github.com/vks/discrete-log
// fn logm(self, base: Modulus, m: Modulus);

/// Collection of common modular arithmetic operations
pub trait ModularOps<Rhs = Self, Modulus = Self, Output = Self>:
    ModularCoreOps<Rhs, Modulus, Output = Output>
    + ModularUnaryOps<Modulus, Output = Output>
    + ModularPow<Rhs, Modulus, Output = Output>
    + ModularSymbols<Modulus>
{
}
impl<T, Rhs, Modulus> ModularOps<Rhs, Modulus> for T where
    T: ModularCoreOps<Rhs, Modulus, Output = T>
        + ModularUnaryOps<Modulus, Output = T>
        + ModularPow<Rhs, Modulus, Output = T>
        + ModularSymbols<Modulus>
{
}

/// Collection of operations similar to [ModularOps], but takes operands with references
pub trait ModularRefOps: for<'r> ModularOps<&'r Self, &'r Self> + Sized {}
impl<T> ModularRefOps for T where T: for<'r> ModularOps<&'r T, &'r T> {}

/// Provides a utility function to convert signed integers into unsigned modular form
pub trait ModularAbs<Modulus> {
    /// Return self % m, but accepting signed integers
    fn absm(self, m: &Modulus) -> Modulus;
}

/// Represents an number defined in a modulo ring ℤ/nℤ
///
/// The operators should panic if the modulus of two number
/// are not the same.
pub trait ModularInteger:
    Sized
    + PartialEq
    + Add<Self, Output = Self>
    + Sub<Self, Output = Self>
    + Neg<Output = Self>
    + Mul<Self, Output = Self>
{
    /// The underlying representation type of the integer
    type Base;

    /// Return the modulus of the ring
    fn modulus(&self) -> Self::Base;

    /// Return the normalized residue of this integer in the ring
    fn residue(&self) -> Self::Base;

    /// Check if the integer is zero
    fn is_zero(&self) -> bool;

    /// Convert an normal integer into the same ring.
    ///
    /// This method should be perferred over the static
    /// constructor to prevent unnecessary overhead of pre-computation.
    fn convert(&self, n: Self::Base) -> Self;

    /// Calculate the value of self + self
    fn double(self) -> Self;

    /// Calculate the value of self * self
    fn square(self) -> Self;
}

// XXX: implement ModularInteger for ff::PrimeField?
// TODO: implement invm_range (Modular inverse in certain range) and crt (Chinese Remainder Theorem), REF: bubblemath crate

/// Utility function for exact division, with precomputed helper values
///
/// # Available Pre-computation types:
/// - `()`: No pre-computation, the implementation relies on native integer division
/// - [PreModInv]: With Pre-computed modular inverse
pub trait DivExact<Rhs, Precompute>: Sized {
    type Output;

    /// Check if d divides self with the help of the precomputation. If d divides self,
    /// then the quotient is returned.
    fn div_exact(self, d: Rhs, pre: &Precompute) -> Option<Self::Output>;
}

/// A modular reducer that can ensure that the operations on integers are all performed
/// in a modular ring.
///
/// Essential information for performing the modulo operation will be stored in the reducer.
pub trait Reducer<T> {
    /// Create a reducer for a modulus m
    fn new(m: &T) -> Self;

    /// Transform a normal integer into reduced form
    fn transform(&self, target: T) -> T;

    /// Check whether target is a valid reduced form
    fn check(&self, target: &T) -> bool;

    /// Get the modulus in original integer type
    fn modulus(&self) -> T;

    /// Transform a reduced form back to normal integer
    fn residue(&self, target: T) -> T;

    /// Test if the residue() == 0
    fn is_zero(&self, target: &T) -> bool;

    /// Calculate (lhs + rhs) mod m in reduced form
    fn add(&self, lhs: &T, rhs: &T) -> T;

    #[inline]
    fn add_in_place(&self, lhs: &mut T, rhs: &T) {
        *lhs = self.add(lhs, rhs)
    }

    /// Calculate 2*target mod m
    fn dbl(&self, target: T) -> T;

    /// Calculate (lhs - rhs) mod m in reduced form
    fn sub(&self, lhs: &T, rhs: &T) -> T;

    #[inline]
    fn sub_in_place(&self, lhs: &mut T, rhs: &T) {
        *lhs = self.sub(lhs, rhs);
    }

    /// Calculate -monty mod m in reduced form
    fn neg(&self, target: T) -> T;

    /// Calculate (lhs * rhs) mod m in reduced form
    fn mul(&self, lhs: &T, rhs: &T) -> T;

    #[inline]
    fn mul_in_place(&self, lhs: &mut T, rhs: &T) {
        *lhs = self.mul(lhs, rhs);
    }

    /// Calculate target^-1 mod m in reduced form,
    /// it may return None when there is no modular inverse.
    fn inv(&self, target: T) -> Option<T>;

    /// Calculate target^2 mod m in reduced form
    fn sqr(&self, target: T) -> T;

    /// Calculate base ^ exp mod m in reduced form
    fn pow(&self, base: T, exp: &T) -> T;
}

mod barret;
mod double;
mod mersenne;
mod monty;
mod preinv;
mod prim;
mod reduced;
mod word;

pub use barret::{
    Normalized2by1Divisor, Normalized3by2Divisor, PreMulInv1by1, PreMulInv2by1, PreMulInv3by2,
};
pub use double::{udouble, umax};
pub use mersenne::FixedMersenne;
pub use monty::Montgomery;
pub use preinv::PreModInv;
pub use reduced::{ReducedInt, Vanilla, VanillaInt};

/// An integer in modulo ring based on [Montgomery form](https://en.wikipedia.org/wiki/Montgomery_modular_multiplication#Montgomery_form)
pub type MontgomeryInt<T> = ReducedInt<T, Montgomery<T>>;

/// An integer in modulo ring with a fixed (pseudo) Mersenne number as modulus
pub type FixedMersenneInt<const P: u8, const K: umax> = ReducedInt<umax, FixedMersenne<P, K>>;

// pub type BarretInt<T> = ReducedInt<T, BarretInt<T>>;

#[cfg(feature = "num-bigint")]
mod bigint;