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.youkuaiyun.com/mutourend/article/details/96426020
        let mut x1 = ProjectivePoint {
            U: affine_u, //X
            W: 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~0
            let 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^2
    let t5 = t1.square();   // (U_P - W_P)^2 = U_P^2 - 2 U_P W_P + W_P^2

    let t6 = &t4 - &t5;     // 4 U_P W_P

    let 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_Q
    let 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_Q

    let 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)^2
    let t12 = t10.square(); // 4 (W_P U_Q - U_P W_Q)^2

    let t13 = &APLUS2_OVER_FOUR * &t6; // (A + 2) U_P U_Q

    let t14 = &t4 * &t5;    // ((U_P + W_P)(U_P - W_P))^2 = (U_P^2 - W_P^2)^2
    let t15 = &t13 + &t5;   // (U_P - W_P)^2 + (A + 2) U_P W_P

    let 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)^2
    let t18 = t11;               // W_D * 4 (U_P U_Q - W_P W_Q)^2

    P.U = t14;  // U_{P'} = (U_P + W_P)^2 (U_P - W_P)^2
    P.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)^2
    Q.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》