1extern crate num;
  2extern crate num_bigint;
  3extern crate num_traits;
  4
  5use num_bigint::{BigInt, ToBigInt};
  6use num_traits::{One, Zero};
  7
  8pub fn modulus(a: &BigInt, m: &BigInt) -> BigInt {
  9    ((a % m) + m) % m
 10}
 11
 12pub fn modinv(a: &BigInt, q: &BigInt) -> Result<BigInt, String> {
 13    let big_zero: BigInt = Zero::zero();
 14    if a == &big_zero {
 15        return Err("no mod inv of Zero".to_string());
 16    }
 17
 18    let mut mn = (q.clone(), a.clone());
 19    let mut xy: (BigInt, BigInt) = (Zero::zero(), One::one());
 20
 21    while mn.1 != big_zero {
 22        xy = (xy.1.clone(), xy.0 - (mn.0.clone() / mn.1.clone()) * xy.1);
 23        mn = (mn.1.clone(), modulus(&mn.0, &mn.1));
 24    }
 25
 26    while xy.0 < Zero::zero() {
 27        xy.0 = modulus(&xy.0, q);
 28    }
 29    Ok(xy.0)
 30}
 31
 32/*
 33pub fn modinv_v2(a0: &BigInt, m0: &BigInt) -> BigInt {
 34    if m0 == &One::one() {
 35        return One::one();
 36    }
 37
 38    let (mut a, mut m, mut x0, mut inv): (BigInt, BigInt, BigInt, BigInt) =
 39        (a0.clone(), m0.clone(), Zero::zero(), One::one());
 40
 41    while a > One::one() {
 42        inv = inv - (&a / m.clone()) * x0.clone();
 43        a = a % m.clone();
 44        std::mem::swap(&mut a, &mut m);
 45        std::mem::swap(&mut x0, &mut inv);
 46    }
 47
 48    if inv < Zero::zero() {
 49        inv += m0.clone()
 50    }
 51    inv
 52}
 53
 54pub fn modinv_v3(a: &BigInt, q: &BigInt) -> BigInt {
 55    let mut aa: BigInt = a.clone();
 56    let mut qq: BigInt = q.clone();
 57    if qq < Zero::zero() {
 58        qq = -qq;
 59    }
 60    if aa < Zero::zero() {
 61        aa = -aa;
 62    }
 63    let d = num::Integer::gcd(&aa, &qq);
 64    if d != One::one() {
 65        println!("ERR no mod_inv");
 66    }
 67    let res: BigInt;
 68    if d < Zero::zero() {
 69        res = d + qq;
 70    } else {
 71        res = d;
 72    }
 73    res
 74}
 75pub fn modinv_v4(x: &BigInt, q: &BigInt) -> BigInt {
 76    let (gcd, inverse, _) = extended_gcd(x.clone(), q.clone());
 77    let one: BigInt = One::one();
 78    if gcd == one {
 79        modulus(&inverse, q)
 80    } else {
 81        panic!("error: gcd!=one")
 82    }
 83}
 84pub fn extended_gcd(a: BigInt, b: BigInt) -> (BigInt, BigInt, BigInt) {
 85    let (mut s, mut old_s) = (BigInt::zero(), BigInt::one());
 86    let (mut t, mut old_t) = (BigInt::one(), BigInt::zero());
 87    let (mut r, mut old_r) = (b, a);
 88
 89    while r != BigInt::zero() {
 90        let quotient = &old_r / &r;
 91        old_r -= &quotient * &r;
 92        std::mem::swap(&mut old_r, &mut r);
 93        old_s -= &quotient * &s;
 94        std::mem::swap(&mut old_s, &mut s);
 95        old_t -= quotient * &t;
 96        std::mem::swap(&mut old_t, &mut t);
 97    }
 98
 99    let _quotients = (t, s); // == (a, b) / gcd
100
101    (old_r, old_s, old_t)
102}
103*/
104
105pub fn concatenate_arrays<T: Clone>(x: &[T], y: &[T]) -> Vec<T> {
106    x.iter().chain(y).cloned().collect()
107}
108
109pub fn modsqrt(a: &BigInt, q: &BigInt) -> Result<BigInt, String> {
110    // Tonelli-Shanks Algorithm (https://en.wikipedia.org/wiki/Tonelli%E2%80%93Shanks_algorithm)
111    //
112    // This implementation is following the Go lang core implementation https://golang.org/src/math/big/int.go?s=23173:23210#L859
113    // Also described in https://www.maa.org/sites/default/files/pdf/upload_library/22/Polya/07468342.di020786.02p0470a.pdf
114    // -> section 6
115
116    let zero: BigInt = Zero::zero();
117    let one: BigInt = One::one();
118    if legendre_symbol(&a, q) != 1 {
119        return Err("not a mod p square".to_string());
120    } else if a == &zero {
121        return Err("not a mod p square".to_string());
122    } else if q == &2.to_bigint().unwrap() {
123        return Err("not a mod p square".to_string());
124    } else if q % 4.to_bigint().unwrap() == 3.to_bigint().unwrap() {
125        let r = a.modpow(&((q + one) / 4), &q);
126        return Ok(r);
127    }
128
129    let mut s = q - &one;
130    let mut e: BigInt = Zero::zero();
131    while &s % 2 == zero {
132        s = s >> 1;
133        e = e + &one;
134    }
135
136    let mut n: BigInt = 2.to_bigint().unwrap();
137    while legendre_symbol(&n, q) != -1 {
138        n = &n + &one;
139    }
140
141    let mut y = a.modpow(&((&s + &one) >> 1), q);
142    let mut b = a.modpow(&s, q);
143    let mut g = n.modpow(&s, q);
144    let mut r = e;
145
146    loop {
147        let mut t = b.clone();
148        let mut m: BigInt = Zero::zero();
149        while &t != &one {
150            t = modulus(&(&t * &t), q);
151            m = m + &one;
152        }
153
154        if m == zero {
155            return Ok(y.clone());
156        }
157
158        t = g.modpow(&(2.to_bigint().unwrap().modpow(&(&r - &m - 1), q)), q);
159        g = g.modpow(&(2.to_bigint().unwrap().modpow(&(r - &m), q)), q);
160        y = modulus(&(y * t), q);
161        b = modulus(&(b * &g), q);
162        r = m.clone();
163    }
164}
165
166#[allow(dead_code)]
167pub fn modsqrt_v2(a: &BigInt, q: &BigInt) -> Result<BigInt, String> {
168    // Tonelli-Shanks Algorithm (https://en.wikipedia.org/wiki/Tonelli%E2%80%93Shanks_algorithm)
169    //
170    // This implementation is following this Python implementation by Dusk https://github.com/dusk-network/dusk-zerocaf/blob/master/tools/tonelli.py
171
172    let zero: BigInt = Zero::zero();
173    let one: BigInt = One::one();
174    if legendre_symbol(&a, q) != 1 {
175        return Err("not a mod p square".to_string());
176    } else if a == &zero {
177        return Err("not a mod p square".to_string());
178    } else if q == &2.to_bigint().unwrap() {
179        return Err("not a mod p square".to_string());
180    } else if q % 4.to_bigint().unwrap() == 3.to_bigint().unwrap() {
181        let r = a.modpow(&((q + one) / 4), &q);
182        return Ok(r);
183    }
184
185    let mut p = q - &one;
186    let mut s: BigInt = Zero::zero();
187    while &p % 2.to_bigint().unwrap() == zero {
188        s = s + &one;
189        p = p >> 1;
190    }
191
192    let mut z: BigInt = One::one();
193    while legendre_symbol(&z, q) != -1 {
194        z = &z + &one;
195    }
196    let mut c = z.modpow(&p, q);
197
198    let mut x = a.modpow(&((&p + &one) >> 1), q);
199    let mut t = a.modpow(&p, q);
200    let mut m = s;
201
202    while &t != &one {
203        let mut i: BigInt = One::one();
204        let mut e: BigInt = 2.to_bigint().unwrap();
205        while i < m {
206            if t.modpow(&e, q) == one {
207                break;
208            }
209            e = e * 2.to_bigint().unwrap();
210            i = i + &one;
211        }
212
213        let b = c.modpow(&(2.to_bigint().unwrap().modpow(&(&m - &i - 1), q)), q);
214        x = modulus(&(x * &b), q);
215        t = modulus(&(t * &b * &b), q);
216        c = modulus(&(&b * &b), q);
217        m = i.clone();
218    }
219    return Ok(x);
220}
221
222pub fn legendre_symbol(a: &BigInt, q: &BigInt) -> i32 {
223    // returns 1 if has a square root modulo q
224    let one: BigInt = One::one();
225    let ls: BigInt = a.modpow(&((q - &one) >> 1), &q);
226    if &(ls) == &(q - one) {
227        return -1;
228    }
229    1
230}
231
232#[cfg(test)]
233mod tests {
234    use super::*;
235
236    #[test]
237    fn test_mod_inverse() {
238        let a = BigInt::parse_bytes(b"123456789123456789123456789123456789123456789", 10).unwrap();
239        let b = BigInt::parse_bytes(b"12345678", 10).unwrap();
240        assert_eq!(
241            modinv(&a, &b).unwrap(),
242            BigInt::parse_bytes(b"641883", 10).unwrap()
243        );
244    }
245
246    #[test]
247    fn test_sqrtmod() {
248        let a = BigInt::parse_bytes(
249            b"6536923810004159332831702809452452174451353762940761092345538667656658715568",
250            10,
251        )
252        .unwrap();
253        let q = BigInt::parse_bytes(
254            b"7237005577332262213973186563042994240857116359379907606001950938285454250989",
255            10,
256        )
257        .unwrap();
258
259        assert_eq!(
260            (modsqrt(&a, &q).unwrap()).to_string(),
261            "5464794816676661649783249706827271879994893912039750480019443499440603127256"
262        );
263        assert_eq!(
264            (modsqrt_v2(&a, &q).unwrap()).to_string(),
265            "5464794816676661649783249706827271879994893912039750480019443499440603127256"
266        );
267    }
268}