curve25519-dalek中field reduce原理分析

对于Curve25519，其Field域内的module F_p = 2²⁵⁵-19。
对于64位系统：

/// A `FieldElement51` represents an element of the field
/// \\( \mathbb Z / (2\^{255} - 19)\\).
///
/// In the 64-bit implementation, a `FieldElement` is represented in
/// radix \\(2\^{51}\\) as five `u64`s; the coefficients are allowed to
/// grow up to \\(2\^{54}\\) between reductions modulo \\(p\\).
///
/// # Note
///
/// The `curve25519_dalek::field` module provides a type alias
/// `curve25519_dalek::field::FieldElement` to either `FieldElement51`
/// or `FieldElement2625`.
///
/// The backend-specific type `FieldElement51` should not be used
/// outside of the `curve25519_dalek::field` module.
#[derive(Copy, Clone)]
pub struct FieldElement51(pub (crate) [u64; 5]);

src/backend/serial/u64/field.rs中的reduce函数，是将[u64;5] low-reduce成h, $h\in [0,2\ast p),p=2^{255}-19$。具体的原理如下：

1. field reduce原理分析

如要求某整数 $umod(2^{255}-19)$，可将u整数用多项式做如下表示：
$u={\sum }_i^{}{u}_i{2}^{51i}{x}^i，其中，u_i\in N$
设置x=1，通过对u多项式求值即可代表域F_p内的值。
如需求一个值：
$a\in [0,2^{320}-1]mod(2^{255}-19)=？$

将$a\in [0,2^{320}-1]$以上图表示，同时根据$a_i$分别取相应的$b_i,c_i$:
$b_i=a_i\&(2<<51-1),c_i=a_i>>51，其中，c_i<=2^{13}$
采用parallel carry-out方式进行，对应有：

以多项式方式表示时，其中的最高项为:
$c_4\ast 2^{255}\ast x^{5}mod(2^{255}-19)$
$\because 2^{255}\ast x^5\equiv 19mod(2^{255}-19)$
$\therefore c_4\ast 2^{255}\ast x^5\equiv c_4\ast 19mod(2^{255}-19)$
所以上图可演化为：

$\because c_i<=2^{13}$
$\because 2^{51}+2^{13}<2^{51}+2^{13}\ast 19<2^{51.0000000001}$
$\because 2^{(51.0000000001\ast 5)}<2\ast (2^{255}-19)=2\ast p,p=2^{255}-19$

sage: (2^51)+(2^13)*19<2^51.0000000001
True
sage: 2^(51.0000000001*5) < 2*(2^255-19)
True

$\therefore h_4=b_4+c_3,h_3=b_3+c_2,h_2=b_2+c_1,h_1=b_1+c_0,h_0=b_0+c_4\ast 19$
$\therefore h={\sum }_{i=0}^4{h}_i\ast 2^{(51.0000000001\ast i)}<2\ast p$

2. field reduce代码实现

/// Given 64-bit input limbs, reduce to enforce the bound 2^(51 + epsilon).
    #[inline(always)]
    fn reduce(mut limbs: [u64; 5]) -> FieldElement51 {
        const LOW_51_BIT_MASK: u64 = (1u64 << 51) - 1;

        // Since the input limbs are bounded by 2^64, the biggest
        // carry-out is bounded by 2^13.
        //
        // The biggest carry-in is c4 * 19, resulting in
        //
        // 2^51 + 19*2^13 < 2^51.0000000001
        //
        // Because we don't need to canonicalize, only to reduce the
        // limb sizes, it's OK to do a "weak reduction", where we
        // compute the carry-outs in parallel.

        let c0 = limbs[0] >> 51;
        let c1 = limbs[1] >> 51;
        let c2 = limbs[2] >> 51;
        let c3 = limbs[3] >> 51;
        let c4 = limbs[4] >> 51;

        limbs[0] &= LOW_51_BIT_MASK;
        limbs[1] &= LOW_51_BIT_MASK;
        limbs[2] &= LOW_51_BIT_MASK;
        limbs[3] &= LOW_51_BIT_MASK;
        limbs[4] &= LOW_51_BIT_MASK;

        limbs[0] += c4 * 19;
        limbs[1] += c0;
        limbs[2] += c1;
        limbs[3] += c2;
        limbs[4] += c3;

        FieldElement51(limbs)
    }

3. field reduce结果$\in [0,2\ast p),p=2^{255}-19$

src/backend/serial/u64/field.rs中的sub,neg等函数都只是调用了reduce函数，将最终的返回结果限定在了$\in [0,2\ast p),p=2^{255}-19$，只有调用to_bytes函数，才会将结果限定在$\in [0,p),p=2^{255}-19$。

		// -1 + (2^255 - 19) = 2^255 - 20 = p-1
	    let a: [u8;32] = [0xec, 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, 0xff, 0xff, 0xff, 0x7f];
    
        let mut b = FieldElement::from_bytes(&a);
        println!("zyd-before negate b:{:?}", b);
        b.conditional_negate(Choice::from(1)); //求倒数，用的是reduce结果，对应结果为：-(p-1) mod p = p+1
        println!("zyd--b:{:?}", b);
        println!("zyd--b.to_bytes():{:?}", b.to_bytes()); //to_bytes()函数会对结果进行再次module，(p+1) mod p = 1.

对应输出为：

zyd-before negate b:FieldElement51([2251799813685228, 2251799813685247, 2251799813685247, 2251799813685247, 2251799813685247])
zyd--b:FieldElement51([2251799813685230, 2251799813685247, 2251799813685247, 2251799813685247, 2251799813685247])
zyd--b.to_bytes():[1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]

4. 一些常量值表示

1）-1对应modulo值为2²⁵⁵-19-1=57896044618658097711785492504343953926634992332820282019728792003956564819948，用FieldElement51表示为：

FieldElement51([2251799813685228, 2251799813685247, 2251799813685247, 2251799813685247, 2251799813685247])

sage: 2251799813685228+2251799813685247*(2^51)+2251799813685247*(2^102)+2251799813685247*(2^153)+2251799813685247*(2^204)==2^255-19-1
True

2）0值对应FieldElement51表示为：

FieldElement51([ 0, 0, 0, 0, 0 ])

3）1值对应FieldElement51表示为：

FieldElement51([ 1, 0, 0, 0, 0 ])

4）16*p值对应FieldElement51表示为：

FieldElement51([ 36028797018963664, 36028797018963952, 36028797018963952, 36028797018963952, 36028797018963952 ])

sage: p=2^255-19
sage: 36028797018963664+36028797018963952*(2^51)+36028797018963952*(2^102)+36028
....: 797018963952*(2^153)+36028797018963952*(2^204)==16*p
True
sage: