有限域内的平方根求解原理解析及curve25519-dalek中的实现

1. 引言

求解$x$值，满足：
$$x^2\equiv a(modp)$$
其中$p$为odd素数。
需要做两层判断：
1）在有限域$p$内，是否存在平方值为$a$的解；
2）若存在解，则相应的解$x$值如何计算？

对于1），可通过Legendre_symbol来判断$a$是否为quadratic residue modulo p.
对于2），可通过1917年的Pocklington算法来计算相应的$x$值。

2. Legendre_symbol

在Montgomery curves 和 Curve25519满足Montgomery curve特征的magma脚本验证中有利用Legendre_symbol来判断判断B(A+2)、B(A-2)、A2-4 三者至少有一个是 module p的 quadratic residue。
可转换为求Legendre symbol值来判断，若该值为-1，则说明不是quadratic residue，若为1，则是quadratic residue。

3. Pocklington算法

若已确定$a$为quadratic residue，在有限域$p$内存在平方根解$x$，可利用Pocklington算法来求解。
Pocklington对不同的有限域值$p$分了以下三种情况进行处理：
1）若$p=4m+3,m\in \mathbb{N}$，则解为$x\equiv \pm a^{m+1}$.
2）若$p=8m+5,m\in \mathbb{N}$，分情况为：

• 当$a^{2m+1}\equiv 1$时，则解为$x\equiv \pm a^{m+1}$.
• 当$a^{2m+1}\equiv -1$，且2为quadratic non-residue所以有$4^{2m+1}\equiv -1$，所以有$(4a{)}^{2m+1}\equiv 1$。$\Rightarrow (\pm (4a{)}^{m+1}{)}^2\equiv 4a$，假设$y\equiv (4a{)}^{m+1}$，则$y为$$y^2\equiv 4a$的解。$\Rightarrow x\equiv \pm y/2或者当y为奇数时，x\equiv \pm (p+y)/2$.

3）若$p=8m+1$，假设$D\equiv -a$，转换为求解$x^2+D\equiv 0$。寻找相应的$t_1和u_1$值，$N=t_1^2-D{u}_1^2$，使得$N$值为quadratic non-residue，即相应的$LegendreSymbol(N,p)$值应为-1。
进一步假设：
$$t_n=\frac{(t_1+u_1\sqrt{D}{)}^n+(t_1-u_1\sqrt{D}{)}^n}{2},u_n=\frac{(t_1+u_1\sqrt{D}{)}^n-(t_1-u_1\sqrt{D}{)}^n}{2\sqrt{D}}$$
从而使得以下等式均成立：
$$t_{m+n}=t_m{t}_n+D{u}_m{u}_n,u_{m+n}=t_m{u}_n+t_n{u}_m,t_n^2-D{u}_n^2=N^n$$
注意，当$p=8m+1$成立时，其实$p=4m^{\prime }+1$亦成立，从而有a为quadratic residue，且-a亦即D也为quadratic residue，$\Rightarrow 存在x^{\prime }使得x^{\prime 2}\equiv D成立$。具体magma脚本验证如下：

clear;
LegendreSymbol(-4, 41);  //1，为quadratic residue
LegendreSymbol(4, 41);  //1，为quadratic residue

注意有限域内有$t^{p-1}\equiv 1$，所以有：
$$t_p\equiv t_1^p\equiv t_1,u_p\equiv u_1^p{D}^{(p-1)/2}\equiv u_1$$
分别取$m=p-1,n=1$，套用到上面的等式从而有：
$$t_p\equiv t_1\equiv t_{p-1}{t}_1+D{u}_{p-1}{u}_1,u_p\equiv u_1\equiv t_{p-1}{u}_1+t_1{u}_{p-1}$$
以上等式成立的解为$t_{p-1}\equiv 1,u_{p-1}\equiv 0$。
因为$p$为odd素数，所以有$p-1=2r$。
分别取$m=r,n=r$，则有：
$$0\equiv u_{p-1}\equiv u_{2r}\equiv t_r{u}_r+t_r{u}_r\equiv 2t_r{u}_r$$
$\Rightarrow t_r或者u_r可整除p.$
假设$u_r$可整除$p$，即$u_r\equiv 0$。若r为偶数，则取$r=2s,m=s,n=s$，从而有$0\equiv 2t_s{u}_s$。但是并不是每一个$u_i$都可整除$p$，如$u_1$就不可以。若r为奇数，因$t_r^2-D{u}_r^2=N^r$等式成立而同时N要求为quadratic non-residue，所以$t_r^2$也应为quadratic non-residue，从而自相矛盾。所以，这个循环将在特定的$l$值处停止，使得$t_l\equiv 0$，从而有：
$-D{u}_l^2\equiv N^l$
因为$-D$亦为quadratic residue而$N$为quadratic non-residue，所以$l$必须为偶数。假设$l=2k,m=k,n=k$，从而有$0\equiv t_l\equiv t_k^2+D{u}_k^2$.
$\because x^2\equiv -D$
$\therefore t_k^2+D{u}_k^2\equiv t_k^2-x^2{u}_k^2\equiv 0$
$\therefore t_k^2\equiv x^2{u}_k^2$
$\therefore u_{k}x\equiv \pm t_k$
$\therefore x\equiv \pm \frac{t_k}{u_k}$

4. curve25519 p=2²⁵⁵-19域sqrt_ratio_i实现

对于curve25519，其$p=2^{255}-19$，满足以上第二个条件，即$p=8m+5,m\in \mathbb{N}$。

• 当$a^{2m+1}\equiv 1$时，则解为$x\equiv \pm a^{m+1}$.
• 当$a^{2m+1}\equiv -1$，且2为quadratic non-residue所以有$4^{2m+1}\equiv -1$，所以有$(4a{)}^{2m+1}\equiv 1$。$\Rightarrow (\pm (4a{)}^{m+1}{)}^2\equiv 4a$，假设$y\equiv (4a{)}^{m+1}$，则$y为$$y^2\equiv 4a$的解。$\Rightarrow x\equiv \pm y/2或者当y为奇数时，x\equiv \pm (p+y)/2$.

其中$m=(p-5)/8,m+1=(p+3)/8$，同时有限域内任意值$a，均有a^{p-1}\equiv 1$。不难理解sqrt_ratio_i函数中的注释，均只取非负数解：

// Using the same trick as in ed25519 decoding, we merge the
        // inversion, the square root, and the square test as follows.
        //
        // To compute sqrt(α), we can compute β = α^((p+3)/8).
        // Then β^2 = ±α, so multiplying β by sqrt(-1) if necessary
        // gives sqrt(α).
        //
        // To compute 1/sqrt(α), we observe that
        //    1/β = α^(p-1 - (p+3)/8) = α^((7p-11)/8)
        //                            = α^3 * (α^7)^((p-5)/8).
        //
        // We can therefore compute sqrt(u/v) = sqrt(u)/sqrt(v)
        // by first computing
        //    r = u^((p+3)/8) v^(p-1-(p+3)/8)
        //      = u u^((p-5)/8) v^3 (v^7)^((p-5)/8)
        //      = (uv^3) (uv^7)^((p-5)/8).
        //
        // If v is nonzero and u/v is square, then r^2 = ±u/v,
        //                                     so vr^2 = ±u.
        // If vr^2 =  u, then sqrt(u/v) = r.
        // If vr^2 = -u, then sqrt(u/v) = r*sqrt(-1).
        //
        // If v is zero, r is also zero.
		// 并通过flipped_sign_sqrt_i 来判断是否有解。替代 Legendre_symbol判断。
    /// - `(Choice(1), +sqrt(u/v))  ` if `v` is nonzero and `u/v` is square;
    /// - `(Choice(1), zero)        ` if `u` is zero;
    /// - `(Choice(0), zero)        ` if `v` is zero and `u` is nonzero;
    /// - `(Choice(0), +sqrt(i*u/v))` if `u/v` is nonsquare (so `i*u/v` is square).
    ///

具体的代码如下：

/// Given `FieldElements` `u` and `v`, compute either `sqrt(u/v)`
    /// or `sqrt(i*u/v)` in constant time.
    ///
    /// This function always returns the nonnegative square root.
    ///
    /// # Return
    ///
    /// - `(Choice(1), +sqrt(u/v))  ` if `v` is nonzero and `u/v` is square;
    /// - `(Choice(1), zero)        ` if `u` is zero;
    /// - `(Choice(0), zero)        ` if `v` is zero and `u` is nonzero;
    /// - `(Choice(0), +sqrt(i*u/v))` if `u/v` is nonsquare (so `i*u/v` is square).
    ///
    pub fn sqrt_ratio_i(u: &FieldElement, v: &FieldElement) -> (Choice, FieldElement) {
        // Using the same trick as in ed25519 decoding, we merge the
        // inversion, the square root, and the square test as follows.
        //
        // To compute sqrt(α), we can compute β = α^((p+3)/8).
        // Then β^2 = ±α, so multiplying β by sqrt(-1) if necessary
        // gives sqrt(α).
        //
        // To compute 1/sqrt(α), we observe that
        //    1/β = α^(p-1 - (p+3)/8) = α^((7p-11)/8)
        //                            = α^3 * (α^7)^((p-5)/8).
        //
        // We can therefore compute sqrt(u/v) = sqrt(u)/sqrt(v)
        // by first computing
        //    r = u^((p+3)/8) v^(p-1-(p+3)/8)
        //      = u u^((p-5)/8) v^3 (v^7)^((p-5)/8)
        //      = (uv^3) (uv^7)^((p-5)/8).
        //
        // If v is nonzero and u/v is square, then r^2 = ±u/v,
        //                                     so vr^2 = ±u.
        // If vr^2 =  u, then sqrt(u/v) = r.
        // If vr^2 = -u, then sqrt(u/v) = r*sqrt(-1).
        //
        // If v is zero, r is also zero.

        let v3 = &v.square()  * v;
        let v7 = &v3.square() * v;
        let mut r = &(u * &v3) * &(u * &v7).pow_p58();
        let check = v * &r.square();

        let i = &constants::SQRT_M1;

        let correct_sign_sqrt   = check.ct_eq(        u);
        let flipped_sign_sqrt   = check.ct_eq(     &(-u));
        let flipped_sign_sqrt_i = check.ct_eq(&(&(-u)*i));

        let r_prime = &constants::SQRT_M1 * &r;
        r.conditional_assign(&r_prime, flipped_sign_sqrt | flipped_sign_sqrt_i);

        // Choose the nonnegative square root.
        let r_is_negative = r.is_negative();
        r.conditional_negate(r_is_negative);

        let was_nonzero_square = correct_sign_sqrt | flipped_sign_sqrt;

        (was_nonzero_square, r)
    }

参考资料：
[1] https://math.stackexchange.com/questions/3280271/compute-sqrtx-pmod-p-working-on-finite-fields
[2] https://en.wikipedia.org/wiki/Pocklington's_algorithm
[3] https://en.wikipedia.org/wiki/Legendre_symbol