curve25519-dalek中的MontgomeryPoint及与Scalar乘积

1. Curve25519定义

The Curve25519 function is Fp-restricted x-coordinate scalar multiplication on E(Fp² ), where p is the prime number p = 2²⁵⁵ − 19 and E is the elliptic curve y² = x³ + 486662x² + x.
对应的base point为(9,…)。
The order of the base point is a prime,
namely 2²⁵² + 27742317777372353535851937790883648493.
FieldElement对应的域空间为[0, p-1)，即为[0, 2²⁵⁵-19 - 1)。

https://github.com/dalek-cryptography/curve25519-dalek/blob/master/src/montgomery.rs 中的相应定义为

pub struct MontgomeryPoint(pub [u8;32]); //采用little-endian方式存储。/// The X25519 basepoint, in `MontgomeryPoint` format.
pub const X25519_BASEPOINT: MontgomeryPoint =MontgomeryPoint([0x09, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00]);/// `BASEPOINT_ORDER` is the order of the Ristretto group and of the Ed25519 basepoint, i.e.,
/// $$
/// \ell = 2^\{252\} + 27742317777372353535851937790883648493.
/// $$
/// sage: hex(2^252+27742317777372353535851937790883648493)
/// '1000000000000000000000000000000014def9dea2f79cd65812631a5cf5d3ed'
pub const BASEPOINT_ORDER: Scalar = Scalar{bytes: [0xed, 0xd3, 0xf5, 0x5c, 0x1a, 0x63, 0x12, 0x58,0xd6, 0x9c, 0xf7, 0xa2, 0xde, 0xf9, 0xde, 0x14,0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00,0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x00, 0x10,],
};

2. MontgomeryPoint*Scalar

/// Multiply this `MontgomeryPoint` by a `Scalar`.
impl<'a, 'b> Mul<&'b Scalar> for &'a MontgomeryPoint {type Output = MontgomeryPoint;/// Given `self` \\( = u\_0(P) \\), and a `Scalar` \\(n\\), return \\( u\_0([n]P) \\).fn mul(self, scalar: &'b Scalar) -> MontgomeryPoint {// Algorithm 8 of Costello-Smith 2017。此处采用的是《Montgomery curves and their arithmetic》论文中的算法8。let affine_u = FieldElement::from_bytes(&self.0); // 将[u8; 32]转换位[u64; 5]数组，且数组内每个元素仅存储51个bit。let mut x0 = ProjectivePoint::identity(); //采用O和P点作为基础点，可保证loop为constant lenght.具体原因参见https://blog.csdn.net/mutourend/article/details/96426020let mut x1 = ProjectivePoint {U: affine_u, //XW: FieldElement::one(), //Z,归一化了};let bits: [i8; 256] = scalar.bits(); //将[u8;32]转换为二进制位按little-endian形式存储[i8;256]，虽然用的i8，其实元素只为0或1。for i in (0..255).rev() { //254~0。scalar的bitlen l=256，循环的范围是从l-2~0let choice: u8 = (bits[i + 1] ^ bits[i]) as u8; //异或debug_assert!(choice == 0 || choice == 1); //保证bits中存储的元素仅为0或1.ProjectivePoint::conditional_swap(&mut x0, &mut x1, choice.into());differential_add_and_double(&mut x0, &mut x1, &affine_u);}ProjectivePoint::conditional_swap(&mut x0, &mut x1, Choice::from(bits[0] as u8));x0.to_affine() //将输入(X:Z)，转换为：$x=X/Z$，以[u8;32]格式存储。}
}

其中FieldElement的定义如下：
FieldElement用数组形式存储[u64;5]，因为Curve25519的p为2²⁵⁵-19，所以FieldElement每个数组元素仅需要用其中的51位即可，分别均匀的存储在5个元素中，即可表示域p内的所有值。

/// A `FieldElement` represents an element of the field
/// \\( \mathbb Z / (2\^{255} - 19)\\).
///
/// The `FieldElement` type is an alias for one of the platform-specific
/// implementations.
#[cfg(feature = "u64_backend")]
pub type FieldElement = backend::serial::u64::field::FieldElement51;pub struct FieldElement51(pub (crate) [u64; 5]);/// Load a `FieldElement51` from the low 255 bits of a 256-bit //256位是因为输入位[u8; 32] = 32*8 = 256 bit。/// input./// 因为p=2^255-19，所以(2^255 - 18) mod p = 1/// # Warning////// This function does not check that the input used the canonical/// representative.  It masks the high bit, but it will happily/// decode 2^255 - 18 to 1.  Applications that require a canonical/// encoding of every field element should decode, re-encode to/// the canonical encoding, and check that the input was/// canonical.///pub fn from_bytes(bytes: &[u8; 32]) -> FieldElement51 {let load8 = |input: &[u8]| -> u64 {(input[0] as u64)| ((input[1] as u64) << 8)| ((input[2] as u64) << 16)| ((input[3] as u64) << 24)| ((input[4] as u64) << 32)| ((input[5] as u64) << 40)| ((input[6] as u64) << 48)| ((input[7] as u64) << 56)};let low_51_bit_mask = (1u64 << 51) - 1;FieldElement51( //将[u8; 32]转换位[u64; 5]数组，且数组内每个元素仅存储51个bit。// load bits [  0, 64), no shift[  load8(&bytes[ 0..])        & low_51_bit_mask// load bits [ 48,112), shift to [ 51,112), (load8(&bytes[ 6..]) >>  3) & low_51_bit_mask// load bits [ 96,160), shift to [102,160), (load8(&bytes[12..]) >>  6) & low_51_bit_mask// load bits [152,216), shift to [153,216), (load8(&bytes[19..]) >>  1) & low_51_bit_mask// load bits [192,256), shift to [204,112), (load8(&bytes[24..]) >> 12) & low_51_bit_mask])}/// Construct zero.pub fn zero() -> FieldElement51 {FieldElement51([ 0, 0, 0, 0, 0 ])}/// Construct one.pub fn one() -> FieldElement51 {FieldElement51([ 1, 0, 0, 0, 0 ])}

ProjectivePoint的定义如下：
其中ProjectivePoint的对应的是F映射（见Montgomery curve的运算（1）——add/double运算 第2节说明）
$$F:P\mapsto \{\begin{array}{ccc}(x_p:1)& if& P=(x_p:y_p:1)\\ (1:0)& if& P=O=(0:1:0)\end{array}$$
即$F((X:Y:Z))=(X:Z)$仅对Z!=0的情况成立。

/// A `ProjectivePoint` holds a point on the projective line
/// \\( \mathbb P(\mathbb F\_p) \\), which we identify with the Kummer
/// line of the Montgomery curve.
#[derive(Copy, Clone, Debug)]
struct ProjectivePoint {pub U: FieldElement,pub W: FieldElement,
}impl Identity for ProjectivePoint {fn identity() -> ProjectivePoint {ProjectivePoint {U: FieldElement::one(),W: FieldElement::zero(),}}
}

let bits: [i8; 256] = scalar.bits();对应为：

	/// Get the bits of the scalar.pub(crate) fn bits(&self) -> [i8; 256] { //将[u8;32]转换为二进制位按little-endian形式存储[i8;256]，虽然用的i8，其实元素只为0或1。let mut bits = [0i8; 256];for i in 0..256 {// As i runs from 0..256, the bottom 3 bits index the bit,// while the upper bits index the byte.bits[i] = ((self.bytes[i>>3] >> (i&7)) & 1u8) as i8; //每个byte内共8个bit元素，元素编号为从0~7，所以需要i>>3和i&7}bits}/// The `Scalar` struct holds an integer \\(s < 2\^{255} \\) which
/// represents an element of \\(\mathbb Z / \ell\\).
#[derive(Copy, Clone)]
pub struct Scalar {/// `bytes` is a little-endian byte encoding of an integer representing a scalar modulo the/// group order.////// # Invariant////// The integer representing this scalar must be bounded above by \\(2\^{255}\\), or/// equivalently the high bit of `bytes[31]` must be zero.////// This ensures that there is room for a carry bit when computing a NAF representation.//// XXX This is pub(crate) so we can write literal constants.  If const fns were stable, we could//     make the Scalar constructors const fns and use those instead.pub(crate) bytes: [u8; 32],
}

differential_add_and_double中，将xADD和xDBL放一个函数内实现，进一步优化了执行效率。

/// Perform the double-and-add step of the Montgomery ladder.
///
/// Given projective points
/// \\( (U\_P : W\_P) = u(P) \\),
/// \\( (U\_Q : W\_Q) = u(Q) \\),
/// and the affine difference
/// \\(      u\_{P-Q} = u(P-Q) \\), set
/// $$
///     (U\_P : W\_P) \gets u([2]P)
/// $$
/// and
/// $$
///     (U\_Q : W\_Q) \gets u(P + Q).
/// $$
fn differential_add_and_double(P: &mut ProjectivePoint,Q: &mut ProjectivePoint,affine_PmQ: &FieldElement,
) {let t0 = &P.U + &P.W;let t1 = &P.U - &P.W;let t2 = &Q.U + &Q.W;let t3 = &Q.U - &Q.W;let t4 = t0.square();   // (U_P + W_P)^2 = U_P^2 + 2 U_P W_P + W_P^2let t5 = t1.square();   // (U_P - W_P)^2 = U_P^2 - 2 U_P W_P + W_P^2let t6 = &t4 - &t5;     // 4 U_P W_Plet t7 = &t0 * &t3;     // (U_P + W_P) (U_Q - W_Q) = U_P U_Q + W_P U_Q - U_P W_Q - W_P W_Qlet t8 = &t1 * &t2;     // (U_P - W_P) (U_Q + W_Q) = U_P U_Q - W_P U_Q + U_P W_Q - W_P W_Qlet t9  = &t7 + &t8;    // 2 (U_P U_Q - W_P W_Q)let t10 = &t7 - &t8;    // 2 (W_P U_Q - U_P W_Q)let t11 =  t9.square(); // 4 (U_P U_Q - W_P W_Q)^2let t12 = t10.square(); // 4 (W_P U_Q - U_P W_Q)^2let t13 = &APLUS2_OVER_FOUR * &t6; // (A + 2) U_P U_Qlet t14 = &t4 * &t5;    // ((U_P + W_P)(U_P - W_P))^2 = (U_P^2 - W_P^2)^2let t15 = &t13 + &t5;   // (U_P - W_P)^2 + (A + 2) U_P W_Plet t16 = &t6 * &t15;   // 4 (U_P W_P) ((U_P - W_P)^2 + (A + 2) U_P W_P)let t17 = affine_PmQ * &t12; // U_D * 4 (W_P U_Q - U_P W_Q)^2let t18 = t11;               // W_D * 4 (U_P U_Q - W_P W_Q)^2P.U = t14;  // U_{P'} = (U_P + W_P)^2 (U_P - W_P)^2P.W = t16;  // W_{P'} = (4 U_P W_P) ((U_P - W_P)^2 + ((A + 2)/4) 4 U_P W_P)Q.U = t18;  // U_{Q'} = W_D * 4 (U_P U_Q - W_P W_Q)^2Q.W = t17;  // W_{Q'} = U_D * 4 (W_P U_Q - U_P W_Q)^2
}

to_affine函数将输入(X:Z)，转换为：$x=X/Z$，以[u8;32]格式存储。

impl ProjectivePoint {/// Dehomogenize this point to affine coordinates.////// # Return////// * \\( u = U / W \\) if \\( W \neq 0 \\);/// * \\( 0 \\) if \\( W \eq 0 \\);pub fn to_affine(&self) -> MontgomeryPoint {let u = &self.U * &self.W.invert(); //此处invert为FieldElement内的invert，乘积为两个FieldElement的乘积。MontgomeryPoint(u.to_bytes())}
}

3. Scalar*MontgomeryPoint

impl<'a, 'b> Mul<&'b MontgomeryPoint> for &'a Scalar {type Output = MontgomeryPoint;fn mul(self, point: &'b MontgomeryPoint) -> MontgomeryPoint {point * self //转换调用2算法，MontgomeryPoint*Scalar}
}

4. MontgomeryPoint *=Scalar

impl<'b> MulAssign<&'b Scalar> for MontgomeryPoint {fn mul_assign(&mut self, scalar: &'b Scalar) {*self = (self as &MontgomeryPoint) * scalar; //转换调用2算法，MontgomeryPoint*Scalar}
}

参考资料：
[1] 论文《Montgomery curves and their arithmetic》