1 use integer::Integer;
2 use traits::Zero;
3
4 use biguint::BigUint;
5
6 struct MontyReducer<'a> {
7 n: &'a BigUint,
8 n0inv: u32,
9 }
10
11 // Calculate the modular inverse of `num`, using Extended GCD.
12 //
13 // Reference:
14 // Brent & Zimmermann, Modern Computer Arithmetic, v0.5.9, Algorithm 1.20
inv_mod_u32(num: u32) -> u3215 fn inv_mod_u32(num: u32) -> u32 {
16 // num needs to be relatively prime to 2**32 -- i.e. it must be odd.
17 assert!(num % 2 != 0);
18
19 let mut a: i64 = i64::from(num);
20 let mut b: i64 = i64::from(u32::max_value()) + 1;
21
22 // ExtendedGcd
23 // Input: positive integers a and b
24 // Output: integers (g, u, v) such that g = gcd(a, b) = ua + vb
25 // As we don't need v for modular inverse, we don't calculate it.
26
27 // 1: (u, w) <- (1, 0)
28 let mut u = 1;
29 let mut w = 0;
30 // 3: while b != 0
31 while b != 0 {
32 // 4: (q, r) <- DivRem(a, b)
33 let q = a / b;
34 let r = a % b;
35 // 5: (a, b) <- (b, r)
36 a = b;
37 b = r;
38 // 6: (u, w) <- (w, u - qw)
39 let m = u - w * q;
40 u = w;
41 w = m;
42 }
43
44 assert!(a == 1);
45 // Downcasting acts like a mod 2^32 too.
46 u as u32
47 }
48
49 impl<'a> MontyReducer<'a> {
new(n: &'a BigUint) -> Self50 fn new(n: &'a BigUint) -> Self {
51 let n0inv = inv_mod_u32(n.data[0]);
52 MontyReducer { n: n, n0inv: n0inv }
53 }
54 }
55
56 // Montgomery Reduction
57 //
58 // Reference:
59 // Brent & Zimmermann, Modern Computer Arithmetic, v0.5.9, Algorithm 2.6
monty_redc(a: BigUint, mr: &MontyReducer) -> BigUint60 fn monty_redc(a: BigUint, mr: &MontyReducer) -> BigUint {
61 let mut c = a.data;
62 let n = &mr.n.data;
63 let n_size = n.len();
64
65 // Allocate sufficient work space
66 c.resize(2 * n_size + 2, 0);
67
68 // β is the size of a word, in this case 32 bits. So "a mod β" is
69 // equivalent to masking a to 32 bits.
70 // mu <- -N^(-1) mod β
71 let mu = 0u32.wrapping_sub(mr.n0inv);
72
73 // 1: for i = 0 to (n-1)
74 for i in 0..n_size {
75 // 2: q_i <- mu*c_i mod β
76 let q_i = c[i].wrapping_mul(mu);
77
78 // 3: C <- C + q_i * N * β^i
79 super::algorithms::mac_digit(&mut c[i..], n, q_i);
80 }
81
82 // 4: R <- C * β^(-n)
83 // This is an n-word bitshift, equivalent to skipping n words.
84 let ret = BigUint::new(c[n_size..].to_vec());
85
86 // 5: if R >= β^n then return R-N else return R.
87 if ret < *mr.n {
88 ret
89 } else {
90 ret - mr.n
91 }
92 }
93
94 // Montgomery Multiplication
monty_mult(a: BigUint, b: &BigUint, mr: &MontyReducer) -> BigUint95 fn monty_mult(a: BigUint, b: &BigUint, mr: &MontyReducer) -> BigUint {
96 monty_redc(a * b, mr)
97 }
98
99 // Montgomery Squaring
monty_sqr(a: BigUint, mr: &MontyReducer) -> BigUint100 fn monty_sqr(a: BigUint, mr: &MontyReducer) -> BigUint {
101 // TODO: Replace with an optimised squaring function
102 monty_redc(&a * &a, mr)
103 }
104
monty_modpow(a: &BigUint, exp: &BigUint, modulus: &BigUint) -> BigUint105 pub fn monty_modpow(a: &BigUint, exp: &BigUint, modulus: &BigUint) -> BigUint {
106 let mr = MontyReducer::new(modulus);
107
108 // Calculate the Montgomery parameter
109 let mut v = vec![0; modulus.data.len()];
110 v.push(1);
111 let r = BigUint::new(v);
112
113 // Map the base to the Montgomery domain
114 let mut apri = a * &r % modulus;
115
116 // Binary exponentiation
117 let mut ans = &r % modulus;
118 let mut e = exp.clone();
119 while !e.is_zero() {
120 if e.is_odd() {
121 ans = monty_mult(ans, &apri, &mr);
122 }
123 apri = monty_sqr(apri, &mr);
124 e >>= 1;
125 }
126
127 // Map the result back to the residues domain
128 monty_redc(ans, &mr)
129 }
130