根据资料【1】中定义,
The Montgomery inverse of an integer a ∈ [1, p−1] is defined by Kaliski [3] as the integer
x
=
a
−
1
2
n
(
m
o
d
p
)
x = a^{−1}2^n (mod\ p)
x=a−12n(mod p). Similarly,
we will use the notation
x
:
=
M
o
n
I
n
v
(
a
)
=
a
−
1
2
n
(
m
o
d
p
)
x := MonInv(a) = a^{−1}2^n (mod\ p)
x:=MonInv(a)=a−12n(mod p)
1. Montgomery inverse
取R=2n,若需求a的倒数x,使得ax = 1 mod p,则有
1)计算 b = aR2。
2)计算montgomery reduction of b:
c = bR-1 mod p = (aR2)R-1 mod p = aR mod p
3)计算montgomery inversion of c:
x:= MonInv© = c-1R mod p = (aR)-1R mod p = a-1 mod p.
根据步骤3),则可获得a的倒数x,使得ax = 1 mod p。
2. 有限域内倒数
有限域内的乘法具有以下特征:
x(p-2) * x = x(p-1) = 1 (mod p)
由此可推测出,求有限域的x值的倒数可转换为求x(p-2)的值。求p-2次幂的运算可基于理论基础为 https://briansmith.org/ecc-inversion-addition-chains-01#curve25519_scalar_inversion 。
3. 实际代码实现
https://github.com/dalek-cryptography/curve25519-dalek/blob/a659b92305cd3286be270cba39e1ae7c8cf59b73/src/scalar.rs 的具体实现如下 :
其中步骤1)和2)对应to_montgomery()函数实现;
步骤3)对应montgomery_invert()和from_montgomery()函数实现。
实际to_montgomery()中针对输入值a,求解的是:
c=(aR2)R-1 mod p = aR mod p.
实际montgomery_invert()中针对输入值c=aR mod p,求解的是:
y = c(p-2)R-(p-3) mod p = (cR-1)(p-2)R mod p = (cR-1)-1R mod p = R2/c mod p.
实际from_montgomery()输入的为y,求解的是:
x = yR-1 mod p = (cR-1)-1 mod p = a-1 mod p.
由此,x即为a的mod p 倒数。
ax=1 mod p.
/// Inverts an UnpackedScalar not in Montgomery form.
pub fn invert(&self) -> UnpackedScalar {
self.to_montgomery().montgomery_invert().from_montgomery()
}
/// Puts a Scalar52 in to Montgomery form, i.e. computes `a*R (mod l)`
#[inline(never)]
pub fn to_montgomery(&self) -> Scalar52 {
Scalar52::montgomery_mul(self, &constants::RR)
}
/// Takes a Scalar52 out of Montgomery form, i.e. computes `a/R (mod l)`
#[inline(never)]
pub fn from_montgomery(&self) -> Scalar52 {
let mut limbs = [0u128; 9];
for i in 0..5 {
limbs[i] = self[i] as u128;
}
Scalar52::montgomery_reduce(&limbs)
}
// 针对步骤3):
/// Inverts an UnpackedScalar in Montgomery form.
pub fn montgomery_invert(&self) -> UnpackedScalar {
// Uses the addition chain from
// https://briansmith.org/ecc-inversion-addition-chains-01#curve25519_scalar_inversion
let _1 = self;
let _10 = _1.montgomery_square();
let _100 = _10.montgomery_square();
let _11 = UnpackedScalar::montgomery_mul(&_10, &_1);
let _101 = UnpackedScalar::montgomery_mul(&_10, &_11);
let _111 = UnpackedScalar::montgomery_mul(&_10, &_101);
let _1001 = UnpackedScalar::montgomery_mul(&_10, &_111);
let _1011 = UnpackedScalar::montgomery_mul(&_10, &_1001);
let _1111 = UnpackedScalar::montgomery_mul(&_100, &_1011);
// _10000
let mut y = UnpackedScalar::montgomery_mul(&_1111, &_1);
#[inline]
fn square_multiply(y: &mut UnpackedScalar, squarings: usize, x: &UnpackedScalar) {
for _ in 0..squarings {
*y = y.montgomery_square();
}
*y = UnpackedScalar::montgomery_mul(y, x);
}
square_multiply(&mut y, 123 + 3, &_101);
square_multiply(&mut y, 2 + 2, &_11);
square_multiply(&mut y, 1 + 4, &_1111);
square_multiply(&mut y, 1 + 4, &_1111);
square_multiply(&mut y, 4, &_1001);
square_multiply(&mut y, 2, &_11);
square_multiply(&mut y, 1 + 4, &_1111);
square_multiply(&mut y, 1 + 3, &_101);
square_multiply(&mut y, 3 + 3, &_101);
square_multiply(&mut y, 3, &_111);
square_multiply(&mut y, 1 + 4, &_1111);
square_multiply(&mut y, 2 + 3, &_111);
square_multiply(&mut y, 2 + 2, &_11);
square_multiply(&mut y, 1 + 4, &_1011);
square_multiply(&mut y, 2 + 4, &_1011);
square_multiply(&mut y, 6 + 4, &_1001);
square_multiply(&mut y, 2 + 2, &_11);
square_multiply(&mut y, 3 + 2, &_11);
square_multiply(&mut y, 3 + 2, &_11);
square_multiply(&mut y, 1 + 4, &_1001);
square_multiply(&mut y, 1 + 3, &_111);
square_multiply(&mut y, 2 + 4, &_1111);
square_multiply(&mut y, 1 + 4, &_1011);
square_multiply(&mut y, 3, &_101);
square_multiply(&mut y, 2 + 4, &_1111);
square_multiply(&mut y, 3, &_101);
square_multiply(&mut y, 1 + 2, &_11);
y
}
4. batch_invert算法
// 本算法要求inputs中所有元素值均为非零。其中返回值为Scalar类型为inputs所有成员乘积的倒数。
// 输入inputs值为mut,返回时其中存储的为每个成员相应的倒数。
#[cfg(feature = "alloc")]
pub fn batch_invert(inputs: &mut [Scalar]) -> Scalar {
// This code is essentially identical to the FieldElement
// implementation, and is documented there. Unfortunately,
// it's not easy to write it generically, since here we want
// to use `UnpackedScalar`s internally, and `Scalar`s
// externally, but there's no corresponding distinction for
// field elements.
use clear_on_drop::ClearOnDrop;
use clear_on_drop::clear::ZeroSafe;
// Mark UnpackedScalars as zeroable.
unsafe impl ZeroSafe for UnpackedScalar {}
let n = inputs.len();
let one: UnpackedScalar = Scalar::one().unpack().to_montgomery(); //为R mod p
// Wrap the scratch storage in a ClearOnDrop to wipe it when
// we pass out of scope.
let scratch_vec = vec![one; n];
let mut scratch = ClearOnDrop::new(scratch_vec);
// Keep an accumulator of all of the previous products
let mut acc = Scalar::one().unpack().to_montgomery(); //为R mod p
// Pass through the input vector, recording the previous
// products in the scratch space //假设inputs=[x_0,x_1,...x_(n-1)]
for (input, scratch) in inputs.iter_mut().zip(scratch.iter_mut()) { //从0~n-1
*scratch = acc; //即为x_0=R mod p, scratch_i=x_(i-1)*x_(i-2)*...*x_1*x_0*R mod p
// Avoid unnecessary Montgomery multiplication in second pass by
// keeping inputs in Montgomery form
let tmp = input.unpack().to_montgomery(); //即为tmp=x_i*R mod p
*input = tmp.pack(); // x_i = x_i*R mod p
acc = UnpackedScalar::montgomery_mul(&acc, &tmp); // 即为acc= R*x_0*x_1....*x_(n-1)*R/R mod p = x_0*x_1....*x_(n-1)*R mod p
}
// acc is nonzero iff all inputs are nonzero
debug_assert!(acc.pack() != Scalar::zero()); //本算法要求inputs中所有元素值均为非零。
// Compute the inverse of all products
acc = acc.montgomery_invert().from_montgomery(); //即为acc=R^2/acc/R mod p = R/acc mod p = 1/(x_0*x_1....*x_(n-1)) mod p
// We need to return the product of all inverses later
let ret = acc.pack(); // 即为ret = 1/(x_0*x_1....*x_(n-1)) mod p
// Pass through the vector backwards to compute the inverses
// in place
for (input, scratch) in inputs.iter_mut().rev().zip(scratch.into_iter().rev()) { //注意此处,从n-1~0,以下先以i=n-1为例,其他同理。
let tmp = UnpackedScalar::montgomery_mul(&acc, &input.unpack()); //即为tmp = acc*x_(n-1)*R/R mod p = 1/(x_0*x_1....*x_(n-2)) mod p
*input = UnpackedScalar::montgomery_mul(&acc, &scratch).pack(); // 即为x_(n-1) = acc*x_(n-2)*x_(n-3)*...*x_1*x_0*R/R mod p = 1/x_(n-1)
acc = tmp; //即为acc = 1/(x_0*x_1....*x_(n-2)) mod p
}
ret
}
参考资料:
[1] http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.75.8377&rep=rep1&type=pdf
[2] https://briansmith.org/ecc-inversion-addition-chains-01#curve25519_scalar_inversion
[3] https://briansmith.org/ecc-inversion-addition-chains-01#curve25519_scalar_inversion
[4] https://github.com/dalek-cryptography/curve25519-dalek/blob/a659b92305cd3286be270cba39e1ae7c8cf59b73/src/scalar.rs