求解 x x x值，满足： x 2 ≡ a ( m o d p ) x^2 \equiv a\quad (mod \quad p) x2≡a(modp) 其中 p p p为odd素数。 需要做两层判断： 1）在有限域 p p p内，是否存在平方值为 a a a的解； 2）若存在解，则相应的解 x x x值如何计算？
对于1），可通过Legendre_symbol来判断 a a a是否为quadratic residue modulo p. 对于2），可通过1917年的Pocklington算法来计算相应的 x x 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。 
若已确定 a a a为quadratic residue，在有限域 p p p内存在平方根解 x x x，可利用Pocklington算法来求解。 Pocklington对不同的有限域值 p p p分了以下三种情况进行处理： 1）若 p = 4 m + 3 , m ∈ N p=4m+3,m\in \mathbb{N} p=4m+3,m∈N，则解为 x ≡ ± a m + 1 x\equiv \pm a^{m+1} x≡±am+1. 2）若 p = 8 m + 5 , m ∈ N p=8m+5,m\in \mathbb{N} p=8m+5,m∈N，分情况为：
- 当 a 2 m + 1 ≡ 1 a^{2m+1}\equiv 1 a2m+1≡1时，则解为 x ≡ ± a m + 1 x\equiv \pm a^{m+1} x≡±am+1.
- 当 a 2 m + 1 ≡ − 1 a^{2m+1}\equiv -1 a2m+1≡−1，且2为quadratic non-residue所以有 4 2 m + 1 ≡ − 1 4^{2m+1}\equiv -1 42m+1≡−1，所以有 ( 4 a ) 2 m + 1 ≡ 1 (4a)^{2m+1}\equiv 1 (4a)2m+1≡1。 ⇒ ( ± ( 4 a ) m + 1 ) 2 ≡ 4 a \Rightarrow (\pm(4a)^{m+1})^2\equiv 4a ⇒(±(4a)m+1)2≡4a，假设 y ≡ ( 4 a ) m + 1 y\equiv(4a)^{m+1} y≡(4a)m+1，则 y 为 y为 y为 y 2 ≡ 4 a y^2\equiv 4a y2≡4a的解。 ⇒ x ≡ ± y / 2 或 者 当 y 为 奇 数 时 ， x ≡ ± ( p + y ) / 2 \Rightarrow x\equiv \pm y/2或者当y为奇数时，x\equiv \pm (p+y)/2 ⇒x≡±y/2或者当y为奇数时，x≡±(p+y)/2.
3）若 p = 8 m + 1 p=8m+1 p=8m+1，假设 D ≡ − a D\equiv -a D≡−a，转换为求解 x 2 + D ≡ 0 x^2+D\equiv 0 x2+D≡0。寻找相应的 t 1 和 u 1 t_1和u_1 t1和u1值， N = t 1 2 − D u 1 2 N=t_1^2-Du_1^2 N=t12−Du12，使得 N N N值为quadratic non-residue，即相应的 L e g e n d r e S y m b o l ( N , p ) LegendreSymbol(N, p) LegendreSymbol(N,p)值应为-1。 进一步假设： t n = ( t 1 + u 1 D ) n + ( t 1 − u 1 D ) n 2 , u n = ( t 1 + u 1 D ) n − ( t 1 − u 1 D ) n 2 D 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} tn=2(t1+u1D )n+(t1−u1D )n,un=2D (t1+u1D )n−(t1−u1D )n 从而使得以下等式均成立： 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 t_{m+n}=t_mt_n+Du_mu_n, u_{m+n}=t_mu_n+t_nu_m, t_n^2-Du_n^2=N^n tm+n=tmtn+Dumun,um+n=tmun+tnum,tn2−Dun2=Nn 注意，当 p = 8 m + 1 p=8m+1 p=8m+1成立时，其实 p = 4 m ′ + 1 p=4m'+1 p=4m′+1亦成立，从而有a为quadratic residue，且-a亦即D也为quadratic residue， ⇒ 存 在 x ′ 使 得 x ′ 2 ≡ D 成 立 \Rightarrow 存在x'使得x'^2\equiv D成立 ⇒存在x′使得x′2≡D成立。具体magma脚本验证如下：
clear;
LegendreSymbol(-4, 41); //1，为quadratic residue
LegendreSymbol(4, 41); //1，为quadratic residue
注意有限域内有 t p − 1 ≡ 1 t^{p-1}\equiv 1 tp−1≡1，所以有： t p ≡ t 1 p ≡ t 1 , u p ≡ u 1 p D ( p − 1 ) / 2 ≡ u 1 t_p\equiv t_1^p\equiv t_1, u_p\equiv u_1^pD^{(p-1)/2}\equiv u_1 tp≡t1p≡t1,up≡u1pD(p−1)/2≡u1 分别取 m = p − 1 , n = 1 m=p-1,n=1 m=p−1,n=1，套用到上面的等式从而有： t p ≡ t 1 ≡ t p − 1 t 1 + D u p − 1 u 1 , u p ≡ u 1 ≡ t p − 1 u 1 + t 1 u p − 1 t_p\equiv t_1\equiv t_{p-1}t_1+Du_{p-1}u_1,u_p\equiv u_1\equiv t_{p-1}u_1+t_1u_{p-1} tp≡t1≡tp−1t1+Dup−1u1,up≡u1≡tp−1u1+t1up−1 以上等式成立的解为 t p − 1 ≡ 1 , u p − 1 ≡ 0 t_{p-1}\equiv 1, u_{p-1}\equiv 0 tp−1≡1,up−1≡0。 因为 p p p为odd素数，所以有 p − 1 = 2 r p-1=2r p−1=2r。 分别取 m = r , n = r m=r,n=r m=r,n=r，则有： 0 ≡ u p − 1 ≡ u 2 r ≡ t r u r + t r u r ≡ 2 t r u r 0\equiv u_{p-1}\equiv u_{2r}\equiv t_ru_r+t_ru_r\equiv 2t_ru_r 0≡up−1≡u2r≡trur+trur≡2trur ⇒ t r 或 者 u r 可 整 除 p . \Rightarrow t_r或者u_r可整除p. ⇒tr或者ur可整除p. 假设 u r u_r ur可整除 p p p，即 u r ≡ 0 u_r\equiv 0 ur≡0。若r为偶数，则取 r = 2 s , m = s , n = s r=2s, m=s,n=s r=2s,m=s,n=s，从而有 0 ≡ 2 t s u s 0\equiv 2t_su_s 0≡2tsus。但是并不是每一个 u i u_i ui都可整除 p p p，如 u 1 u_1 u1就不可以。若r为奇数，因 t r 2 − D u r 2 = N r t_r^2-Du_r^2=N^r tr2−Dur2=Nr等式成立而同时N要求为quadratic non-residue，所以 t r 2 t_r^2 tr2也应为quadratic non-residue，从而自相矛盾。所以，这个循环将在特定的 l l l值处停止，使得 t l ≡ 0 t_l\equiv 0 tl≡0，从而有： − D u l 2 ≡ N l -Du_l^2\equiv N^l −Dul2≡Nl 因为 − D -D −D亦为quadratic residue而 N N N为quadratic non-residue，所以 l l l必须为偶数。假设 l = 2 k , m = k , n = k l=2k,m=k,n=k l=2k,m=k,n=k，从而有 0 ≡ t l ≡ t k 2 + D u k 2 0\equiv t_l\equiv t_k^2+Du_k^2 0≡tl≡tk2+Duk2. ∵ x 2 ≡ − D \because x^2\equiv-D ∵x2≡−D ∴ t k 2 + D u k 2 ≡ t k 2 − x 2 u k 2 ≡ 0 \therefore t_k^2+Du_k^2\equiv t_k^2-x^2u_k^2\equiv 0 ∴tk2+Duk2≡tk2−x2uk2≡0 ∴ t k 2 ≡ x 2 u k 2 \therefore t_k^2\equiv x^2u_k^2 ∴tk2≡x2uk2 ∴ u k x ≡ ± t k \therefore u_kx\equiv \pm t_k ∴ukx≡±tk ∴ x ≡ ± t k u k \therefore x\equiv \pm \frac{t_k}{u_k} ∴x≡±uktk
4. curve25519 p=2255-19域sqrt_ratio_i实现对于curve25519，其 p = 2 255 − 19 p=2^{255}-19 p=2255−19，满足以上第二个条件，即 p = 8 m + 5 , m ∈ N p=8m+5,m\in \mathbb{N} p=8m+5,m∈N。
- 当 a 2 m + 1 ≡ 1 a^{2m+1}\equiv 1 a2m+1≡1时，则解为 x ≡ ± a m + 1 x\equiv \pm a^{m+1} x≡±am+1.
- 当 a 2 m + 1 ≡ − 1 a^{2m+1}\equiv -1 a2m+1≡−1，且2为quadratic non-residue所以有 4 2 m + 1 ≡ − 1 4^{2m+1}\equiv -1 42m+1≡−1，所以有 ( 4 a ) 2 m + 1 ≡ 1 (4a)^{2m+1}\equiv 1 (4a)2m+1≡1。 ⇒ ( ± ( 4 a ) m + 1 ) 2 ≡ 4 a \Rightarrow (\pm(4a)^{m+1})^2\equiv 4a ⇒(±(4a)m+1)2≡4a，假设 y ≡ ( 4 a ) m + 1 y\equiv(4a)^{m+1} y≡(4a)m+1，则 y 为 y为 y为 y 2 ≡ 4 a y^2\equiv 4a y2≡4a的解。 ⇒ x ≡ ± y / 2 或 者 当 y 为 奇 数 时 ， x ≡ ± ( p + y ) / 2 \Rightarrow x\equiv \pm y/2或者当y为奇数时，x\equiv \pm (p+y)/2 ⇒x≡±y/2或者当y为奇数时，x≡±(p+y)/2.
其中
m
=
(
p
−
5
)
/
8
,
m
+
1
=
(
p
+
3
)
/
8
m=(p-5)/8, m+1=(p+3)/8
m=(p−5)/8,m+1=(p+3)/8，同时有限域内任意值
a
，
均
有
a
p
−
1
≡
1
a，均有a^{p-1}\equiv 1
a，均有ap−1≡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
