求二次剩余(Quadratic-Residue)
参考文档《ANSI X9.62 ECDSA 1998》
p ≡ 3 mod 4
在素数域 Fp 上，当 p ≡ 3 mod 4时，可以利用 费马小定理 快速计算模平方根，从而高效求解椭圆曲线方程 y2=x3+ax+b mod p 中的 y。
可通过以下公式计算：
y = a<sup>(p+1)/4</sup> mod p (这里的a可以理解为x<sup>3</sup>+ax+b)
推导：
由费马小定理，ap-1 ≡ 1 mod p (这里p为素数，逆元唯一，逆元是a(p-1)/2本身，即a(p-1)/2 mop = 1)
因为a是二次剩余（即存在y使得y<sup>2</sup> ≡ a mod p），所以
(a<sup>(p+1)/4</sup>)<sup>2</sup> = a<sup>(p+1)/2</sup> = a * a<sup>(p-1)/2</sup> = a * (a/p) = a * 1 = a mod p
（其中 (a/p)是勒让德符号，当 a 是二次剩余时等于 1）
//secp256r1 y^2 = x^3 + ax + b
//p
FFFFFFFF00000001000000000000000000000000FFFFFFFFFFFFFFFFFFFFFFFF
//a
FFFFFFFF00000001000000000000000000000000FFFFFFFFFFFFFFFFFFFFFFFC
//b
5AC635D8AA3A93E7B3EBBD55769886BC651D06B0CC53B0F63BCE3C3E27D2604B
//Gx
6B17D1F2E12C4247F8BCE6E563A440F277037D812DEB33A0F4A13945D898C296//actually no use
//Gy
4FE342E2FE1A7F9B8EE7EB4A7C0F9E162BCE33576B315ECECBB6406837BF51F5//actually no use
//n
FFFFFFFF00000000FFFFFFFFFFFFFFFFBCE6FAADA7179E84F3B9CAC2FC632551
//h
1
//给定x，求y(23EB7389BC00546F6081EF82F2C3A6C543DCB7E46883F6AA440C1A4F80E21367)
31AFBD653F927CB14B8B869487FB9E045E93316B08213956B867278012364731
//求y
//求y-step1：求a=x^3^+ax+b
a=x^3+ax+b
=31AFBD653F927CB14B8B869487FB9E045E93316B08213956B867278012364731^3+FFFFFFFF00000001000000000000000000000000FFFFFFFFFFFFFFFFFFFFFFFC*31AFBD653F927CB14B8B869487FB9E045E93316B08213956B867278012364731+5AC635D8AA3A93E7B3EBBD55769886BC651D06B0CC53B0F63BCE3C3E27D2604B mod FFFFFFFF00000001000000000000000000000000FFFFFFFFFFFFFFFFFFFFFFFF
=8AB7B34BB1184807A9AB1161A9ED89DC5C864FE4FA16CAB9541D77B76074EBA4+6AF0C7CF414889ED1D5D6C42680D25F2E4466BBFE79C53FBD6CA897FC95D2A6C+5AC635D8AA3A93E7B3EBBD55769886BC651D06B0CC53B0F63BCE3C3E27D2604B mod FFFFFFFF00000001000000000000000000000000FFFFFFFFFFFFFFFFFFFFFFFF
=F5A87B1AF260D1F4C7087DA411FAAFCF40CCBBA4E1B31EB52AE8013729D21610++5AC635D8AA3A93E7B3EBBD55769886BC651D06B0CC53B0F63BCE3C3E27D2604B mod FFFFFFFF00000001000000000000000000000000FFFFFFFFFFFFFFFFFFFFFFFF
=506EB0F49C9B65DB7AF43AF98893368BA5E9C254AE06CFAB66B63D7551A4765C
//求y-step2：求y
a^((p+1)/4) mod p
=506EB0F49C9B65DB7AF43AF98893368BA5E9C254AE06CFAB66B63D7551A4765C^3FFFFFFFC0000000400000000000000000000000400000000000000000000000 mod FFFFFFFF00000001000000000000000000000000FFFFFFFFFFFFFFFFFFFFFFFF
=23EB7389BC00546F6081EF82F2C3A6C543DCB7E46883F6AA440C1A4F80E21367
//y
23EB7389BC00546F6081EF82F2C3A6C543DCB7E46883F6AA440C1A4F80E21367
p ≡ 5 mod 8
在素数域 Fp 上，当 p ≡ 5 mod 8时，即p=8u+5，可以通过以下步骤求二次剩余y^2 = a mod p (求y）
设g=a，
1.Compute Y = (2g)u mod p
2.Compute i = 2gY2 mod p
3.Compute y = gY(i-1) mod p
4.Compute z = y2 mod p
5.If z==g，output success，otherwise output "no square root exist"
//sm9p256v1 y^2 = x^3 + ax + b
//p
B640000002A3A6F1D603AB4FF58EC74521F2934B1A7AEEDBE56F9B27E351457D
//a
0000000000000000000000000000000000000000000000000000000000000000
//b
0000000000000000000000000000000000000000000000000000000000000005
//Gx
93DE051D62BF718FF5ED0704487D01D6E1E4086909DC3280E8C4E4817C66DDDD//actually no use
//Gy
21FE8DDA4F21E607631065125C395BBC1C1C00CBFA6024350C464CD70A3EA616//actually no use
//n
B640000002A3A6F1D603AB4FF58EC74449F2934B18EA8BEEE56EE19CD69ECF25
//h
1
//给定x，求y(81396BEC33E7CE207DCC9B4DC765A51F139A71CB251B27CA9D2500B3E7C79C49)
943BD683FF6B7E59EDD0825AE79DF7178070F10DAE759770E7ED3538620A509C
//求u，p=8u+5
16C80000005474DE3AC07569FEB1D8E8A43E5269634F5DDB7CADF364FC6A28AF
//求g=a=x^3^+ax+b
=943BD683FF6B7E59EDD0825AE79DF7178070F10DAE759770E7ED3538620A509C^3+0*943BD683FF6B7E59EDD0825AE79DF7178070F10DAE759770E7ED3538620A509C+5 mod B640000002A3A6F1D603AB4FF58EC74521F2934B1A7AEEDBE56F9B27E351457D
=54ABA8C8BE6C4E8ECFD89BB444CC97C93748321F9BCA8A14BFAF00C70EC962CD+5 mod B640000002A3A6F1D603AB4FF58EC74521F2934B1A7AEEDBE56F9B27E351457D
=54ABA8C8BE6C4E8ECFD89BB444CC97C93748321F9BCA8A14BFAF00C70EC962D2
//Y = (2g)^u^ mod p
=(2*54ABA8C8BE6C4E8ECFD89BB444CC97C93748321F9BCA8A14BFAF00C70EC962D2)^16C80000005474DE3AC07569FEB1D8E8A43E5269634F5DDB7CADF364FC6A28AF mod B640000002A3A6F1D603AB4FF58EC74521F2934B1A7AEEDBE56F9B27E351457D
=A95751917CD89D1D9FB1376889992F926E90643F379514297F5E018E1D92C5A4^16C80000005474DE3AC07569FEB1D8E8A43E5269634F5DDB7CADF364FC6A28AF mod B640000002A3A6F1D603AB4FF58EC74521F2934B1A7AEEDBE56F9B27E351457D
=28365169A25D69166EA7BA44E597434F194CF02C623C9E45C92AB60415D91CE3
//i = 2gY^2^ mod p
=(2*54ABA8C8BE6C4E8ECFD89BB444CC97C93748321F9BCA8A14BFAF00C70EC962D2)*(28365169A25D69166EA7BA44E597434F194CF02C623C9E45C92AB60415D91CE3^2) mod p
=A95751917CD89D1D9FB1376889992F926E90643F379514297F5E018E1D92C5A4*15AF7A05C1F5A246DC2C5DFC74BE58C3C15727D5F21EB1622C24FB63B889B5FD mod p
=6C648DE5DC0A3F2CF55ACC93EE0BAF159F9D411806DC5177F5B21FD3DA24D011
//y = gY(i-1) mod p
=54ABA8C8BE6C4E8ECFD89BB444CC97C93748321F9BCA8A14BFAF00C70EC962D2*28365169A25D69166EA7BA44E597434F194CF02C623C9E45C92AB60415D91CE3*(6C648DE5DC0A3F2CF55ACC93EE0BAF159F9D411806DC5177F5B21FD3DA24D011-1) mod p
=7006D6BA5A972CAB589909AAF351FAB145CE7260A10CDBF9DA047D56197D848C*6C648DE5DC0A3F2CF55ACC93EE0BAF159F9D411806DC5177F5B21FD3DA24D010 mod p
=35069413CEBBD8D1583710022E2922260E58217FF55FC711484A9A73FB89A934
//y
35069413CEBBD8D1583710022E2922260E58217FF55FC711484A9A73FB89A934
//p-y
81396BEC33E7CE207DCC9B4DC765A51F139A71CB251B27CA9D2500B3E7C79C49
p ≡ 1 mod 4
在素数域 Fp 上，当 p ≡ 1 mod 4时，即p=4u+1，可以通过以下步骤求二次剩余y^2 = a mod p (求y）
1.Set Q=g=a
2.Generate random P with 0<=P<p
3.compute the Lucas sequence element
U=U2u+1 mod p, V=V2u+1 mod p
4..If V2 ≡ 4Q (mod p) then output y=V/2 mod p and stop
5.If U ≢ ±1 (mod p) then output the message "no square roots exist" and stop
Generating Lucas Sequences
1.Let P and Q be nonzero integers.The Lucas sequences Uk and Vk for P, Q are defined by:
U0 = 0, U1 = 1, and Uk = PUk-1 - QUk-2 for k>=2
V0 = 2, V1 = P, and Vk = PVk-1 - QVk-2 for k>=2
This recursion is adequate for computing Ukand Vk for small values of k. The following algorithm can be used to
efficiently compute Uk and Vk modulo an odd prime p for large values of k.
Input: An odd prime p, integers P and Q, and a positive integer k
Output:Uk mod p and Vk mod p
1.
二、椭圆曲线点计算
y2 = x3 + ax + b ————(1)
椭圆曲线点加
//计算点(x1,y2)+(x2,y2) = (x3, y3)
x3 = s2 - x1 - x2 mod p
y3 = s(x1 - x3) - y1 mod p
//斜率s
s = (y2 - y1)/(x2 - x1) mod p 当P≠Q （不同点相加）
s = (3x12 + a)/(2y1) mod p 当P=Q（相同点相加）
推导
1.P、Q为不同点的斜率
k = (y2 - y1)/(x2 - x1) ————(2)
2.P、Q为相同点的斜率
对曲线表达式两边关于x求导，得
2y*(dy/dx) = 3x2 + a，得
y' = dy/dx = (3x2 + a)/(2y) ————(3)
y'就是曲线上任意一点的切线的斜率
通过点(x1,y1)的直线L(x)可以表示为：
y = k(x - x1) + y1 ————(4)
求(x3,y3)坐标-横坐标x3的计算
将(4)带入(1)，消去y，得
(k(x-x1)+y1)2 = x3 + ax + b
k2(x-x1)2 + 2k(x-x1)y1 + y12 = x3 + ax + b
k2x2 - 2k2xx1 + k2x12 + 2kxy - 2kx1y1 + y12 = x3 + ax + b
k2x2 - 2k2xx1 + k2x12 + 2kxy - 2kx1y1 + y12 - x3 - ax - b = 0 //将等式右边移到左边
-k2x2 + 2*k2xx1 - k2x12 - 2kxy + 2kx1y1 - y12+ x3 + ax + b = 0 //取反
x3 - k2x2 + (a + 2k2x1 - 2ky1)x + b - (k2x12 - 2kx1y1 + y12) = 0 //合并x
————(5)
根据韦达定理(Vieta's Formulas)，三次方程ax3 + bx + cx + d = 0的三个根x1，x2，x3有以下关系
x1 + x2 + x3 = -b/a
x1x2 + x1x3 + x2x3 = -c/a
x1x2x3 = -d/a
对应(5)，x1 + x2 + x3 = -b/a = k2/1
因此，x3 = k2 - x1 - x2
求(x3,y3)坐标-横坐标y3的计算
方法：通过两点斜率直接计算得到，通过所求R(x3, y3)的镜像点R‘(x3, -y3)
k = (y1 - (-y3))/(x1 - x3)，得
y3 = k(x1 - x3) - y1
因此，当P!=Q，计算R
$$
\begin{align}
k &= \frac{y_2 - y_1}{x_2 - x_1} \
x_3 &= k^2 - x_1 - x_2 \
y_3 &= k(x_1 - x_3) - y_1
\end{align}
$$
当进行点倍(P+P)运算时，求2P，只是斜率不同，后续过程是相同的
$$
\begin{align}
k &= \frac{3x_1^2 + a}{2y_1} \
x_3 &= k^2 - 2x_1 \
y_3 &= k(x_1 - x_3) - y_1
\end{align}
$$
//secp256r1 示例demo
//p
FFFFFFFF00000001000000000000000000000000FFFFFFFFFFFFFFFFFFFFFFFF
//a
FFFFFFFF00000001000000000000000000000000FFFFFFFFFFFFFFFFFFFFFFFC
//b
5AC635D8AA3A93E7B3EBBD55769886BC651D06B0CC53B0F63BCE3C3E27D2604B
//P=[1]G
6B17D1F2E12C4247F8BCE6E563A440F277037D812DEB33A0F4A13945D898C296
4FE342E2FE1A7F9B8EE7EB4A7C0F9E162BCE33576B315ECECBB6406837BF51F5
//Q=[2]G
7CF27B188D034F7E8A52380304B51AC3C08969E277F21B35A60B48FC47669978
07775510DB8ED040293D9AC69F7430DBBA7DADE63CE982299E04B79D227873D1
//P+Q
5ECBE4D1A6330A44C8F7EF951D4BF165E6C6B721EFADA985FB41661BC6E7FD6C
8734640C4998FF7E374B06CE1A64A2ECD82AB036384FB83D9A79B127A27D5032
//计算P+P
//k = (3x^2+a)/2y mod p
=(3*6B17D1F2E12C4247F8BCE6E563A440F277037D812DEB33A0F4A13945D898C296*6B17D1F2E12C4247F8BCE6E563A440F277037D812DEB33A0F4A13945D898C296+FFFFFFFF00000001000000000000000000000000FFFFFFFFFFFFFFFFFFFFFFFC)/(2*4FE342E2FE1A7F9B8EE7EB4A7C0F9E162BCE33576B315ECECBB6406837BF51F5) mod p
=(3*98F6B84D29BEF2B281819A5E0E3690D833B699495D694DD1002AE56C426B3F8C+FFFFFFFF00000001000000000000000000000000FFFFFFFFFFFFFFFFFFFFFFFC)/(9FC685C5FC34FF371DCFD694F81F3C2C579C66AED662BD9D976C80D06F7EA3EA) mod p
=(CAE428E87D3CD8168484CF1A2AA3B2889B23CBDB183BE9730080B044C741BEA5+FFFFFFFF00000001000000000000000000000000FFFFFFFFFFFFFFFFFFFFFFFC)*7D13D1ED160030C54147C66B2E0D48CF7FE33DFB45A6DD82DDED53BAE22EA01A mod p
=CAE428E87D3CD8168484CF1A2AA3B2889B23CBDB183BE9730080B044C741BEA2*7D13D1ED160030C54147C66B2E0D48CF7FE33DFB45A6DD82DDED53BAE22EA01A mod p
=A249BCB3191802D3CB56FC75A8731F848E7B2CAF080F3477A4E4D904A671AF9D
//x3 = k^2 - 2x1 mod p
=(A249BCB3191802D3CB56FC75A8731F848E7B2CAF080F3477A4E4D904A671AF9D*A249BCB3191802D3CB56FC75A8731F848E7B2CAF080F3477A4E4D904A671AF9D)-(2*6B17D1F2E12C4247F8BCE6E563A440F277037D812DEB33A0F4A13945D898C296) mod p
=53221EFF4F5BD40D7BCC05CDCBFD9CA8AE9064E3D3C882778F4DBB87F8981EA5-D62FA3E5C258848FF179CDCAC74881E4EE06FB025BD66741E942728BB131852C mod p
=7CF27B188D034F7E8A52380304B51AC3C08969E277F21B35A60B48FC47669978
//y3 = k(x1-x3)-y1
=A249BCB3191802D3CB56FC75A8731F848E7B2CAF080F3477A4E4D904A671AF9D*(6B17D1F2E12C4247F8BCE6E563A440F277037D812DEB33A0F4A13945D898C296-7CF27B188D034F7E8A52380304B51AC3C08969E277F21B35A60B48FC47669978)-4FE342E2FE1A7F9B8EE7EB4A7C0F9E162BCE33576B315ECECBB6406837BF51F5 mod p
=A249BCB3191802D3CB56FC75A8731F848E7B2CAF080F3477A4E4D904A671AF9D*EE2556D95428F2CA6E6AAEE25EEF262EB67A139FB5F9186B4E95F0499132291D-4FE342E2FE1A7F9B8EE7EB4A7C0F9E162BCE33576B315ECECBB6406837BF51F5
=575A97F3D9A94FDBB82586111B83CEF1E64BE13DA81AE0F869BAF8055A37C5C6-4FE342E2FE1A7F9B8EE7EB4A7C0F9E162BCE33576B315ECECBB6406837BF51F5 mod p
=7775510DB8ED040293D9AC69F7430DBBA7DADE63CE982299E04B79D227873D1
//-y3 mod p
= 0-7775510DB8ED040293D9AC69F7430DBBA7DADE63CE982299E04B79D227873D1 mod p
=F888AAEE24712FC0D6C26539608BCF244582521AC3167DD661FB4862DD878C2E
//计算P+Q mod p
//k = (y2-y1)/(x2-x1) mod p
=(07775510DB8ED040293D9AC69F7430DBBA7DADE63CE982299E04B79D227873D1-4FE342E2FE1A7F9B8EE7EB4A7C0F9E162BCE33576B315ECECBB6406837BF51F5)*(7CF27B188D034F7E8A52380304B51AC3C08969E277F21B35A60B48FC47669978-6B17D1F2E12C4247F8BCE6E563A440F277037D812DEB33A0F4A13945D898C296)^(-1) mod p
=B794122CDD7450A59A55AF7C236492C58EAF7A8FD1B8235AD24E7734EAB921DB*(11DAA925ABD70D369195511DA110D9D14985EC614A06E794B16A0FB66ECDD6E2)^(-1) mod p
=B794122CDD7450A59A55AF7C236492C58EAF7A8FD1B8235AD24E7734EAB921DB*5CD14F7B0DB3AE51C4B8A837B2F5EFD55EA52A15F4802742D3248C306AE52D7 mod p
=28CDB0FE6681D140C59FB73F85F16459CA4CEE083742BADA762DE4CA226A8086
//x3 = k^2 - x1 -x2
=28CDB0FE6681D140C59FB73F85F16459CA4CEE083742BADA762DE4CA226A8086*28CDB0FE6681D140C59FB73F85F16459CA4CEE083742BADA762DE4CA226A8086-6B17D1F2E12C4247F8BCE6E563A440F277037D812DEB33A0F4A13945D898C296-7CF27B188D034F7E8A52380304B51AC3C08969E277F21B35A60B48FC47669978 mod p
=46D631DE14629C0A4C070E7D85A54D1C1E539E84958AF85C95EDE85DE6E7597B-6B17D1F2E12C4247F8BCE6E563A440F277037D812DEB33A0F4A13945D898C296-7CF27B188D034F7E8A52380304B51AC3C08969E277F21B35A60B48FC47669978 mod p
=DBBE5FEA333659C3534A279822010C29A7502104679FC4BBA14CAF180E4E96E4-7CF27B188D034F7E8A52380304B51AC3C08969E277F21B35A60B48FC47669978 mod p
=5ECBE4D1A6330A44C8F7EF951D4BF165E6C6B721EFADA985FB41661BC6E7FD6C
//y3 = k(x1-x3)-y1
=28CDB0FE6681D140C59FB73F85F16459CA4CEE083742BADA762DE4CA226A8086*(6B17D1F2E12C4247F8BCE6E563A440F277037D812DEB33A0F4A13945D898C296-5ECBE4D1A6330A44C8F7EF951D4BF165E6C6B721EFADA985FB41661BC6E7FD6C)- 4FE342E2FE1A7F9B8EE7EB4A7C0F9E162BCE33576B315ECECBB6406837BF51F5 mod p
=28CDB0FE6681D140C59FB73F85F16459CA4CEE083742BADA762DE4CA226A8086*C4BED213AF938032FC4F75046584F8C903CC65F3E3D8A1AF95FD32A11B0C52A-4FE342E2FE1A7F9B8EE7EB4A7C0F9E162BCE33576B315ECECBB6406837BF51F5 mod p
=D717A6EF47B37F19C632F2189674410303F8E38DA381170C662FF18FDA3CA227-4FE342E2FE1A7F9B8EE7EB4A7C0F9E162BCE33576B315ECECBB6406837BF51F5 mod p
=8734640C4998FF7E374B06CE1A64A2ECD82AB036384FB83D9A79B127A27D5032
//-y3 mod p
=0-8734640C4998FF7E374B06CE1A64A2ECD82AB036384FB83D9A79B127A27D5032 mod p
=78CB9BF2B6670082C8B4F931E59B5D1327D54FCAC7B047C265864ED85D82AFCD
Affine 仿射坐标
椭圆曲线的仿射方程：
$$
y^2 = x^3 + a x + b
$$
在有限域中，方程：
$$
y^2 = x^3 + a x + b \mod p
$$
在防射坐标下，存在需要使用拓展的欧几里得算法求解乘法逆元的问题
标准射影坐标
在防射坐标中，我们有椭圆曲线方程y2 = x3 + ax + b，对应的点坐标为(x,y)，仿射坐标是三维，点表示为(x,y,z)
在标准射影坐标中，点(x,y,z)对应仿射坐标中的点(x/z,y/z)，同理标准射影坐标(x,y,1)对应仿射坐标(x,y)
将x=x/z，y=y/z带入到仿射方程，得标准射影坐标椭圆曲线方程：
$$
\begin{align}
\left(\frac{y}{z}\right)^2 &= \left(\frac{x}{z}\right)^3 + a\left(\frac{x}{z}\right) + b \notag \
\frac{y2}{z2} &= \frac{x3}{z3} + a \frac{x}{z} + b \notag \
y^2 &= x^3 + axz^2+ b z^3 \quad \text{(两边同乘 $z^3$)}
\end{align}
$$
由于点(x,y,z)对应对应仿射坐标的点(x/z, y/z)，这使得z不能为0，如果z为0，在射影坐标中的点(0,y,0)找不到一个对应的仿射坐标值，但实际上椭圆曲线的点还是以射影坐标为准，也就是说有这样一个点(0,y,0)，在仿射坐标中被称为”无穷远点“，对应映射坐标中z等于0的点。
令P和Q为不相等的两点，P(x1,y1,z1)，Q(x2,y2,z2)，另R(x3,y3,z3)，且R=P+Q，标准射影坐标下的几何加法运算如下：
$$
\begin{align}
λ_1 &= x_1z_2 \
λ_2 &= x_2z_1 \
λ_3 &= λ_1 -λ_2 \
λ_4 &= y_1z_2 \
λ_5 &= y_2z_1 \
λ_6 &= λ_4 -λ_5 \
λ_7 &= λ_1 + λ_2 \
λ_8 &= z_1z_2 \
λ_9 &= z_3^2 \
λ_{10} &= λ_3λ_9 \
λ_{11} &= λ_8λ_6^2 - λ_7λ_9 \
x_3 &= λ_3λ_{11} \
y_3 &= λ_6(λ_9λ_1 - λ_{11}) - λ_4λ_{10} \
z_3 &= λ_{10}λ_8
\end{align}
$$
最后得到的x3，y3，z3即为得到的R点标准射影坐标
对于标量乘法，计算P+P，令P(x1,y1,z1)，R(x3,y3,z3)，标准射影坐标运算过程如下：
$$
\begin{align}
λ_1 &= 3x_12+az_12 \
λ_2 &= 2y_1z_1 \
λ_3 &= y_1^2 \
λ_4 &= λ_3x_1z_1 \
λ_5 &= λ_2^2\
λ_6 &= λ_1^2 - 8λ_4 \
x_3 &= λ_2λ_6 \
y_3 &=λ_1(4λ_4-λ_6)-2λ_5λ_3 \
z_3 &= λ_2λ_5
\end{align}
$$
最后得到的x3，y3，z3即为得到的R点标准射影坐标
Jacobian加重射影坐标
与标准射影坐标方法类似，在Jacobian加重射影坐标中，一个射影点的坐标表示为(x,y,z)，这里规定：
在Jacobian加重射影坐标，点(x,y,z)对应仿射坐标中的点(x/(z2),y/(z3))
同理，在仿射坐标中的点(x,y)对应jacobian加重射影坐标中的点(x,y,1)
将x=x/(z2), y=y/(z3)，代入椭圆曲线方程，得到Jacobian加重射影坐标的椭圆曲线方程
$$
\begin{align}
\left(\frac{y}{z3}\right)2 &= \left(\frac{x}{z2}\right)3 + a\left(\frac{x}{z^2}\right) + b \notag \
\frac{y2}{z6} &= \frac{x3}{z6} + a \frac{x}{z^2} + b \notag \
y^2 &= x^3 + axz^4+ b z^6 \quad \text{(两边同乘 $z^6$)}
\end{align}
$$
同理，当z为0，有射影坐标中的点(x,y,0)，该点同样在仿射坐标中被称为"无穷远点"，即在射影坐标中z=0的点
蒙哥马利模乘
蒙哥马利约减：REDC(x,y)=x * y * R-1
求x * y mod N
1）REDC(x, R2) = xR
2）REDC(y, R2) = yR
3）REDC(xR, yR) = xyR
4）REDC(xyR, 1) = xy
计算约简：TR-1 mod N ≡ T/R mod N
变形：TR-1 mod N ≡ T/R mod N ≡ (T+mN)/R mod N
同时，要使T+mN ≡ 0 mod R，求m
获取：
令N0 ≡ -N-1 mod R
=> -N0 ≡ N-1 mod R
=> -N0N ≡ 1 mod R ----①
令R‘ ≡ R-1 mod N
=>R'R ≡ 1 mod N -----②
由①得：
-N0N + R'R ≡ 1 mod R
T+mN ≡ 0 mod R
=> T(-N0N + R'R)+mN ≡ 0 mod R
=> -TN0N + TR'R + mN ≡ 0 mod R
=> mN ≡ TN0N - TR'R mod R
=> mN ≡TN0N mod R
=> m ≡ TN0 mod R
简单乘法：32560/1000
=((((32 * 0/10) +32 * 6)/10)+325)/10
计算ABR-1 mod N = AB/R mod N
//A和B为超大数，也可以直接(AB+ABN0N)/R mod N，本身大数计算，没有计算优势
从低位分组i=0,1,...,ceil(n_bitlen/mgmr_bitlen)，开始，原分组：res = res + A * Bi
改成蒙哥马利形式：令T = res + A * Bi => TR-1 = (res + A * Bi)/R
①计算m=TN0 mod R = (res + A * Bi) N0 mod R //取低位
②计算TR-1 = (T + mN)/R = [(res+AB1) + (res+ABi)NN0]/R //右移
蒙哥马利模约减demo:
求a mod n，假设n有效bit为n_bits，a有效bit为a_bits，a_bits > n_bits
a mod n
= a + k * n mod n //寻找一个k，使得 a+kn的低n_bits为全0
= a1 * 2n_bits + a0 +kn mod n //将a拆成a1和a0，a0占用n_bits，a1是高位
= > a1 * 2n_bits + a0 +kn ≡ 0 mod 2n_bits
= > a1 * 2n_bits + a0 +kn = x * 2n_bits
= > a0 +kn = (x -a1 ) 2n_bits ≡ 0 mod 2n_bits //x、k未知
= > k = -a0 * n-1 mod 2n_bits
令t = a + k * n (≡ 0 mod 2n_bits) ，R = 2n_bits
t = a + k * n = (a + k * n) * R-1 = a mod n //只需要一次加法、一次乘法
t = t * R-1 = a mod n
=> t = a * R-1 mod n //蒙哥马利约减
多倍点运算实现
1）二进制展开法：
计算Q=[k]P，将k=knkn-1...k0，kn=1 //最高位是有效位
1）Q=P
2）for j = n-1 to 0：
2.1）Q=[2]Q
b.2）if kj=1 then Q=Q+P
c）return Q
霍纳法则：多项式求值运算中的应用
大数运算
加法
减法
乘法
void uint32_Mul(uint32_t* a, uint32_t a_words, uint32_t* b, uint32_t b_words, uint32_t* c, uint32_t* c_words)
{
uint64_t tmp;
uint32_t high;
uint32_t low;
uint32_t carry;
uint32_t i, j;
uint32_clear(c, a_words + b_words);
for (i = 0; i < a_words; i++)
{
carry = 0;
for (j = 0; j < b_words; j++)
{
//(2^32-1)^2 < (2^32)^2 = 2^64
tmp = (uint64_t)a[i] * (uint64_t)b[j];
high = tmp >> 32;
low = (tmp << 32) >> 32;
//此时需要添加上一进位值，以及低32位
tmp = (uint64_t)c[i + j] + (uint64_t)low + (uint64_t)carry;
//计算下一进位值
carry = high + (tmp >> 32);
c[i + j] = (tmp << 32) >> 32;
}
//写入进位，有可能为0
c[i + j] = carry;
}
*c_words = a_words + b_words;
while (0 == c[(*c_words) - 1])
{
(*c_words)--;
}
}
模
大数模小数，利用模乘和模加实现
模加
参考《Guide to Elliptic Curve Cryptography》第二章
模减
模乘
模逆
一种大数除法和取模算法原理和实现
在嵌入式开发中，除法和取模(/和%)执行起来比较慢，一般应尽量避免使用。执行除法或者模运算时，当除数是常数和用同一个除
数重复除法，执行效率会比较高。本文将阐述如何用乘法运算代替除法运算，以及如何减少除法次数。
用乘法代替除法
假设d是U8类型，n是U8类型。约定:
$$
n/d:表示计算机中的除法,比如4/5=0; \
\frac{n}{d}:表示数学意义中的除法,比如4/5=0.8
$$
取A=0xFF，估计n * (A/d)>>8与n/d之间的误差，并计算n%d。
$$
\begin{align}
(1)记 s &= A/d = 0xFF /d \
&=> 2^8 -1 = s * d + t \quad(0 \leq t < d) \
&=> s = \frac{2^8 - 1 - t}{d} = \frac{2^8-1}{d} - \frac{t}{d} \quad (0 \leq \frac{t}{d} < 1) \
&=> s = \frac{2^8 - 1 - t}{d} = \frac{2^8}{d} - \frac{1+t}{d} =\frac{2^8}{d} - e_1 \quad (0 < e_1 \leq 1) \
\end{align}
$$
数学上，一个数除以2^n，与该数右移n位，之间会产生误差e2，其中e2范围[0,1)
$$
\begin{align}
(2)记 t&=(ns)>>8 \
&=> ns\frac{1}{2^8} - e_2 \quad(0\leq e_2 <1) \
&=> n(\frac{2^8}{d} - e_1)\frac{1}{2^8} - e_2 \quad(将(1)中s代入) \
&=> \frac{n}{d} - \frac{ne_1}{2^8} -e_2 \quad(0 < e_1 \leq1,\quad 0\leq e_2 <1) \quad Ⅰ\ \
&\because 0 < e_1 \leq1,\quad 0\leq e_2 <1 \
&\therefore 0< \frac{n*e_1}{2^8} + e_2<2 \quad(0 < \frac{n}{2^8} <1) \quad Ⅱ\ \
&=>\frac{n}{d}-2<t<\frac{n}{d}-0 \quad(结合ⅠⅡ) \quadⅢ\ \
&\because \frac{n}{d} = n/d + e_3 \quad(0\leq e_3 <1) \quad Ⅳ\
&\therefore n/d + e_3-2<t
\end{align}
$$
因此，d是U8类型，预计算0xFF/d后，n * (0xFF/d)>>8可以用来代替除法n/d，从而避免除法。d是其它数据类型的情况可以类似得到。例如，d为U32类型，预计算0xFFFFFFFFFFFF/d，除法n / d可以用n * (0xFFFFFFFFF/d)>>32来代替。
不过，q=n * (0xFF/d)>>8与n/d有可能存在偏差(当d是U8类型)，该偏差可以通过计算余数来校正。计算r=n - q * d，若r≥d，说明q=n/d-1，因此q需要加1才是准确的n/d，r则要减去d才是准确的n%d。
减少除法运算的次数
考虑10进制下的模运算 2551%23
思路1：自然的计算方式，先计算25%23=2，然后25%23=2，最后21%23=21
2551=2 * 103+5 * 102+5 * 101+1 * 100
=(2 * 102+5 * 101+5)10 + 1
=[(2 * 101+5) * 10+5 ]10+1
思路2：预计算100/23=4，100%23=100-423=8
计算25%23=2，然后2 * 8+51=67，67%23=21 （2551=25100+51，2 * 8对应(25%23) * (100%23)）
对于大数运算,思路2比思路1减少了除法(模运算)的次数。
考虑计算中，100比23多了一位，而在8位机或者32位机中，从空间和运算性能上更多地是考虑一个U8类型的数模上另一个U8类型的数，或者一个U32类型的数模上另一个U32类型的数，因此在思路2的基础上，得到思路3。
思路3：预计算s=99/23=4，99%23=99-23*4=7
计算25%23=2，然后2 * (7 + 1)+51=67，67%23=21
大数取模运算的实现
上述3个思路都有取模运算，根据1，取模运算可以用乘法代替。
思路4：如下计算2551%23
（1）预计算s=99/23=4，t=99%23=99-234=7
（2）对于25，为了计算25%23，先计算25/23
2.1 计算q=25 * s=254/100=1，r=25-1 * 23=2<23，所以25%23=r=2
2.2 计算2 * (t+1)=2(7+1)+51=67
2.3 计算q=67 * 4/100=2，r=67-223=21<23，所以2551%23=r=21
注：在2.3步骤中，存在溢出U8的可能，需要额外处理，如"2*(7+1)+91=16+91=107",回绕数据便为99+1+7，如果溢出(res为回绕后的值)可以处理为
res = res + (7 + 1) = 7+8=15，也符合结果。
证明：
res=a*r+b （a/r/b均为32位数，res最大值=(232-1) * (232-1) + (232-1) = 264 - 232）
=t1||t0 + b (计算机运算过程隐式对t1||t0截断为t0)
如下计算2585%23（r=9）
（1）预计算s=99/23=4，t=99%23=99-234=7, t0=t+1=8
（2）对于25，为了计算25%23，先计算25/23
2.1 计算q=25 * s=254/100=1，r=25-1 * 23=2<23，所以25%23=r=2
2.2 计算2 * (t+1)=2(7+1)+85=101，因为101 = (99+1)+1 = (99%23+1)+1%23=(t+1)+1=(7+1)+1=9
2.3 计算q=9 * 4/100=0，r=9-023=9<23，所以r=9
如下计算2189%23（r=4）
（1）预计算s=99/23=4，t=99%23=99-234=7, t0=t+1=8
（2）对于21，为了计算21%23，先计算21/23
2.1 计算q=21 * s=214/100=0，r=21-0 * 23=21<23，所以21%23=r=21
2.2 计算21 * (t+1)=21(7+1)+89=168+89=257
2.2.1 数学运算：253 = (2 * 99+2 * 1)+57 =2 * t0 + 57=73
2.2.2 计算q=73 * 4/100=2，r=73-223=27<23，r=27-23=4, 所以r=4
2.2.1 计算机运算，21 * 8+89=168+89=68（168截断）+89=157=57（截断）
如下计算2185%23（r=0）
（1）预计算s=99/23=4，t=99%23=99-234=7, t0=t+1=8
（2）对于21，为了计算21%23，先计算21/23
2.1 计算q=21 * s=214/100=0，r=21-0 * 23=21<23，所以21%23=r=21
2.2 计算21 * (t+1)=21(7+1)+89=168+85=253，
2.2.1 数学运算253 = (2 * 99+2 * 1)+53 = (2 * 99%23+2 * 1)+53%23=2 * (t+1)+53=28+53=69
2.2.2 计算q=69 * 4/100=2，r=69-2*23=23<23，r = r-23=0, 所以r=0
在RSA大数除小素数表中，会用到上述算法。
32位机器：
rsa代码demo
/* function: get a%prime, prime number count is PTL(internal API)
* parameters:
* a -------------------------- input, uint32_t big odd number a
* aWordLen ------------------- input, word length of a
* r -------------------------- input, r value of PrimeTable
* s -------------------------- input, s value of PrimeTable
* high_result ---------------- output, high result of a%prime, prime number count is PTL
* PTL ------------------------ input, number of little prime
* return:
* PKE_SUCCESS(success); other(error)
* caution:
*/
FLAG_STATIC uint32_t bigint_div_table_high(const uint32_t *const a, const uint32_t aWordLen, const uint16_t *const r,
const double_uint32_t *const s, uint16_t *const high_result, const uint32_t PTL)
{
int32_t j;
uint32_t i;
uint32_t carry, a1, a2, ret = PKE_SUCCESS;
#ifdef CONFIG_PKE_SUPPORT_TASK_DROP
uint32_t is_suspended;
#endif
/*
The following implementation is a fast way to find the remainder of a large integer to a small prime number(a%p):
1.assume a is u32 array, p,t,carry is u32, a1,a2 is u64, a=q*p+r (0<=r<p), a=am||...||a1||a0, a0 is the most significant word
2.pre value s=0xFFFFFFFF/p, r=0xFFFFFFFF%p+1
3.for(i=aWordLen-1;i!=0;i--)
{
a1=carry*r+a[i];
a1=(a1>>32)*r+((a1<<32)>>32);
a2=(a1*s)>>32;
carry=a1-a2*p;
if (carry>=d)
{
carry=carry-p;
}
}
4.return carry
here s[i] is 0xFFFFFFFF/PrimeTable[i], stored in constant PrimeTable_s[i]
here r[i] is 0xFFFFFFFF%PrimeTable[i]+1 , stored in constant PrimeTable_r[i]
*/
for(i=0U; i<PTL; i++)
{
carry = 0U;
for (j = ((int32_t)aWordLen) - 1; j >= 0; j--)
{
/* PRQA S 3384 ++
here carry represents the remainder, carry < PrimeTable[i], and carry and PrimeTable is U16 unsign int,
0 <= carry*(uint32_t)r[i] < 2^32-2^17+1, so carry*(uint32_t)r[i] will not overflow in U32
we assume u32 t = carry*(uint32_t)r[i], so 0<=t+a[j]<2^33-2^17+1<2^32+2^32
if overflow in u32, t+a[j]=0xffffffff+1+x, (t+a[j])%PrimeTable[i] = 0xffffffff%PrimeTable[i]+1+x = r[i]+x
So if there is overflow, we could correct it by add r[i], still ensure a1% PrimeTable[i] equal to (carry*r[i]+a[j])%PrimeTable[i]
*/
a1=(carry*(uint32_t)r[i])+(uint32_t)a[j]; /* big number is represented in little-enidan */
if(a1<a[j])
{
a1+=(uint32_t)r[i];
}
/* PRQA S 3384 -- #
/* PRQA S 3384 ++ #
the '/' symbol represents computer division.
we assume a[i]=PrimeTable[i]*quotient+r (0<=r<PrimeTable[i]), quotient=a[i]/PrimeTable[i]
The following operation calculates the quotient by using (a[i]*s[i])>>32 to replace a[i]/PrimeTable[i],
that use multiplication and right shift instead of computer division.
Mathematically, there will be an error between dividing a number by 2^n and right-shifting the number by n bits.
Here, a2 in the code may not be the actual quotient.there are two possible cases for quotient:
1.a2 = a[i]/PrimeTable[i]
2.a2 +1 = a[i]/PrimeTable[i]
so the multiplication and right shift here are intentional by the algorithm.
*/
a2=(uint32_t)((((uint64_t)a1)*s[i])>>SHIFT_32_BIT);
/* PRQA S 3384 -- #*/
/* PRQA S 3383 ++ #
From the above (a[i]*s[i])>>32, we can know (a2*PrimeTable[i])<=a1, so a1-(a2*PrimeTable[i]) will not overflow
The following operation get carry=a[i]%PrimeTable[i], becase of quotient(here a2 in code), there are two cases of carry:
1.carry = a[i]%PrimeTable[i]
2.carry = a[i]%PrimeTable[i] + PrimeTable[i]
*/
carry=a1-(a2*(uint32_t)PrimeTable[i]);
/* PRQA S 3383 -- #*/
if(carry>=(uint32_t)PrimeTable[i])
{
//correct the remainder of a[i]%PrimeTable[i]
carry -= (uint32_t)PrimeTable[i];
}
}
high_result[i]=(uint16_t)carry;
}
return ret;
}