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

use super::UInt;
use crate::Limb;

impl<const LIMBS: usize> UInt<LIMBS> {
    /// Computes 1/`self` mod 2^k as specified in Algorithm 4 from
    /// A Secure Algorithm for Inversion Modulo 2k by
    /// Sadiel de la Fe and Carles Ferrer. See
    /// <https://www.mdpi.com/2410-387X/2/3/23>.
    ///
    /// Conditions: `self` < 2^k and `self` must be odd
    pub const fn inv_mod2k(&self, k: usize) -> Self {
        let mut x = Self::ZERO;
        let mut b = Self::ONE;
        let mut i = 0;

        while i < k {
            let mut x_i = Self::ZERO;
            let j = b.limbs[0].0 & 1;
            x_i.limbs[0] = Limb(j);
            x = x.bitor(&x_i.shl_vartime(i));

            let t = b.wrapping_sub(self);
            b = Self::ct_select(b, t, j.wrapping_neg()).shr_vartime(1);
            i += 1;
        }
        x
    }
}

#[cfg(test)]
mod tests {
    use crate::U256;

    #[test]
    fn inv_mod2k() {
        let v = U256::from_be_slice(&[
            0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff,
            0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xfe,
            0xff, 0xff, 0xfc, 0x2f,
        ]);
        let e = U256::from_be_slice(&[
            0x36, 0x42, 0xe6, 0xfa, 0xea, 0xac, 0x7c, 0x66, 0x63, 0xb9, 0x3d, 0x3d, 0x6a, 0x0d,
            0x48, 0x9e, 0x43, 0x4d, 0xdc, 0x01, 0x23, 0xdb, 0x5f, 0xa6, 0x27, 0xc7, 0xf6, 0xe2,
            0x2d, 0xda, 0xca, 0xcf,
        ]);
        let a = v.inv_mod2k(256);
        assert_eq!(e, a);

        let v = U256::from_be_slice(&[
            0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff, 0xff,
            0xff, 0xfe, 0xba, 0xae, 0xdc, 0xe6, 0xaf, 0x48, 0xa0, 0x3b, 0xbf, 0xd2, 0x5e, 0x8c,
            0xd0, 0x36, 0x41, 0x41,
        ]);
        let e = U256::from_be_slice(&[
            0x26, 0x17, 0x76, 0xf2, 0x9b, 0x6b, 0x10, 0x6c, 0x76, 0x80, 0xcf, 0x3e, 0xd8, 0x30,
            0x54, 0xa1, 0xaf, 0x5a, 0xe5, 0x37, 0xcb, 0x46, 0x13, 0xdb, 0xb4, 0xf2, 0x00, 0x99,
            0xaa, 0x77, 0x4e, 0xc1,
        ]);
        let a = v.inv_mod2k(256);
        assert_eq!(e, a);
    }
}