secp256k1算法详解三（点操作关键理论及源码分析）
1 基础概念1.1 Short Weierstrass Curve
椭圆曲线（Elliptic Curve，EC）是密码学中非常重要的代数结构，在几何上，椭圆曲线是由三次方程描述的一些曲线。不同的公式给出了不同曲线形式，对于维尔斯特拉斯形式的椭圆曲线，其定义在域K上的椭圆曲线E由下述方程定义：
例如椭圆曲线y2 + xy + y = x3 - x + 1图形如下所示：
详细内容可以参考：张焕国著《椭圆曲线密码学导论》。由于一般形式的维尔斯特拉斯方程系数太多不利于分析，一般多基于短维尔斯特拉斯曲线（当域K的特征≠2，3时，一般形式维尔斯特拉斯可以简化为短维尔斯特拉斯）进行研究，它是有以下方程描述的椭圆曲线：
y2 = x3 + ax + b
参数a，b满足条件：4a3+27b2≠0，限定条件是为了保证曲线不包含奇点（singularities某些特定方面不再表现良好的点，例如缺乏可微性或可分析性）。在椭圆曲线密码学中x，y以及参数a，b都属于有限域Fp，如对于secp256k1曲线，参数a=0，b=7，即公式进一步简化为y2 = x3 + 7。
1.2 点加及倍点
在椭圆曲线上定义了两个点的加法运算，两个点P和点Q之和点为R，则R是过P，Q点直线与椭圆曲线交点关于x轴的对称点，这里P和Q的横坐标不同时，当P和Q重合时（横纵坐标都相同），则R是过P点切线与椭圆曲线交点关于x轴的对称点，如下图所示：
为了理论的完整性，还引入无穷远点的概念，这时对应P和Q横坐标相同纵坐标不同，则它们和点即为无穷远点（此时P和Q互为对方逆元），这里不再详细介绍相关内容，请自行查阅有关资料。
1.3 仿射坐标(Affine Coordinates)
维尔斯特拉斯形式曲线上点的坐标就是椭圆曲线点的仿射坐标(x, y)，接下来基于仿射坐标分析下点加及倍点运算。对于椭圆曲线上两个点P(x1, y1)，Q(x2, y2)，令R(x3, y3) = P + Q，当x1 ≠ x2时，过P，Q点直线坐标为y=λx + m，代入到椭圆曲线方程得
x3 - λ2x2 + (a - 2λm)x + b - m2 = 0
设直线与椭圆曲线交点为R'(x', y')，则可知x1，x2，x'是以上一元三次方程的根，由韦达定理可知：-x1 - x2 - x' = -λ2，则有：
其中λ是直线斜率。因为R是R'关于x轴对称点，所以R坐标为(λ2 - x1 - x2, (2x1 + x2 - λ2)λ -y1)，带入λ有以下计算公式：
当P和Q重合时（且y ≠ 0），对椭圆曲线公式两边关于x求导得斜率：
进而得倍点公式：
1.4 射影坐标(Projective Coordinates)
在数学和密码学领域，射影坐标（Projective Coordinates） 是一种通过引入额外维度来避免除法运算的坐标表示方法，在传统的仿射坐标（Affine Coordinates） 中，平面上的点用二维坐标(x, y)表示，对应笛卡尔坐标系中的位置。但仿射坐标在处理涉及除法的运算（如椭圆曲线的点加法、点加倍）时，会引入计算成本较高的逆元操作（有限域中的除法等价于乘以逆元）。射影坐标通过引入一个额外的非零参数Z，将仿射坐标(x, y)扩展为三维坐标(X, Y, Z)，其中Z ≠ 0，且满足x = X / Z，y = Y / Z（或其他等价变换，具体形式取决于坐标系定义）。这一变换的核心是用乘法和加法替代除法：通过将仿射坐标中的除法转化为射影坐标中对Z分量的乘法操作，避免了计算逆元的高成本步骤，尤其在有限域等代数结构中能显著提升运算效率。
1.4.1 标准射影坐标（Standard Projective Coordinates，简称P3)
定义：对于仿射点(x, y)，射影坐标表示为(X, Y, Z)，满足：x = X/Z，y = Y/Z，其中Z ≠ 0。
等价性：若λ ≠ 0，则(X, Y, Z)与(λX, λY, λZ)表示同一个仿射点（即坐标具有齐次性）。
1.4.2 雅可比射影坐标（Jacobian Projective Coordinates，简称J3)
定义：为进一步优化椭圆曲线运算（尤其是点加倍）而设计，满足：x = X/Z2，y = Y/Z3，其中Z ≠ 0。
优势：在椭圆曲线的点加倍运算中，能减少乘法和平方操作的次数，是 ECC 实现中最常用的坐标形式之一。
2 公式推导
由于secp256k1射影坐标用的是雅可比坐标，所以只介绍其公式推导。
2.1 倍点公式推导
雅可比坐标用(X, Y, Z)表示仿射坐标(x, y) = (X/Z2, Y/Z3)，需要将仿射坐标得倍点公式转换为雅可比坐标。
2.1.1 将斜率用雅可比表示
仿射斜率：
2.1.2 计算X3, Y3, Z3
将以上以及x1雅可比表示代入公式（3）中有：
进一步进行通分得：
同理由公式（4）可得：
观察以上两式，只要令：
则有x3 = X3/Z32，y3 = Y3/Z33，完全满足推导出来的等式，所以以上即雅可比倍点公式。
2.1.3 优化计算步骤
计算中间变量：
则X3，Y3，Z3由以下计算公式给出：
其实以上D=4X1Y12，该计算方式只有在(X1+B)2比单独计算 X1⋅B更快才有价值，否则可直接使用无D版本公式：
对于secp256k1来说因为a=0，所以hyperelliptic.org给出了两种简化后的计算对比，M：乘法（Multiplication），指有限域上的元素乘法操作；S：平方（Squaring），指有限域上的元素平方操作（通常比一般乘法成本低，故单独统计）；add：加法（Addition），指有限域上的元素加法或减法操作（成本较低，有时可忽略，但此处单独列出）；3*2、1*3、1*8：表示特定常数乘法的次数（如乘以常数 2、3、8 的操作，这类操作可通过优化实现，成本通常低于一般乘法）。因此，第一个公式表示该算法共需：2次乘法 + 5次平方 + 6次加法 + 3次乘以2的操作 + 1次乘以3的操作 + 1次乘以8的操作。
2.2 点加公式推导
2.2.1 x3雅可比坐标表示
雅可比坐标用(X, Y, Z)表示仿射坐标(x, y) = (X/Z2, Y/Z3)，将雅可比表示的坐标直接代入公式（1）可知：
2.2.2 y3雅可比坐标表示
求y3时，由于：y3 = λ(x1 - x3) - y1 且 y3 = λ(x2 - x3) - y2，等式两侧分别相加可知：
等式两边除于2即得y3。
2.2.3 两点和公式
观察以上得出x3和y3表达式，只要令：
则有：
为了使得计算流程更为清晰，可引入如下中间变量：
则求解公式可简化为：
对于secp256k1来说因为a=0，hyperelliptic.org给出详细的计算复杂度，如下所示：
在源码中使用的正是add-1998-cmo方法，很容易验算它和add-1986-cc方法是一致的。
2.3 统一倍点及点和公式
2.3.1 侧信道攻击
维基百科是这么定义的：“在密码学中，侧信道攻击（Side-channel attack）是一种攻击方式，它基于从密码系统的物理实现中获取的信息，而非暴力破解法或是算法中的理论性缺陷，例如利用时间信息、功率消耗、电磁泄露或甚是声音可以提供额外的信息，来对系统的破解”。简单一点理解就是利用一些“旁门左道”的手法来获取我们需要的机密信息；比如经常在电视剧中看到的一个场景：一个窃贼将听诊器压在保险柜的前面板上，通过内部的机械声来打开保险柜。小偷会慢慢地转动转盘，听着内部机械结构所泄露出的咔哒声或阻力声，来分析保险箱齿轮的内部运作，并从而得知其密码的组合。除了拨号盘上的数字和保险柜“是”或“否”的打开状态以外，这个保险柜并不会给用户任何反馈，但保险箱的物理机械所产生的那些微小的触动和声音线索，这也是一个典型的侧信道攻击。
SPA-Like Attacks 通常是指类似简单功率分析（Simple Power Analysis，SPA）的攻击方式。SPA 是一种侧信道攻击方法，它通过分析密码设备在运行过程中的功率消耗情况，来获取设备操作和密钥等信息。其攻击原理是：密码设备在执行不同操作时，其功率消耗会有所不同。例如，乘法运算通常比加法运算消耗更多的功率，对不同数据的处理也会导致不同的功率消耗模式。攻击者通过测量密码设备的功率消耗，观察其功率波形，就可以推断出设备正在执行的操作，进而分析出密钥等敏感信息。在椭圆曲线密码体制中，标量乘算法常采用 “倍点 - 点加” 算法。由于倍点和点加运算的时间和消耗能量不同，攻击者可以通过分析功率波形，识别出这两种运算，从而获取关于标量（通常与密钥相关）的信息，对椭圆曲线密码系统构成威胁。标量乘法算法如下，如果从功率追踪显示倍点运算后又进行了一次加法运算，则当前位ki应该等于1，否则等于0。
2.3.2 统一公式
分析倍点及点和公式，它们最大的不同就在于关联直线（“chord-and-tangent”弦与切线）的斜率λ表示不同，如果能用相同的公式表达λ，则公式（1）、（2）和公式（3）、（4）就能统一起来。
首先给出一个基于一般维尔斯特拉斯方程的命题：E是基于域K的椭圆曲线，其方程为：y2 + a1xy + a3y = x3 + a2x2 + a4x + a6，同样P = (x1, y1)，Q = (x2, y2)是椭圆曲线上的点，且y(P) ≠ y(-Q)（保证两点不关于曲线对称，避免点加退化为倍点或无穷远点），令P + Q = (x3, y3)，则有：
下面给出以上结论的证明过程，对于Q(x2, y2)假设其负元-Q坐标为(x', y')，则有x'=x2，即它们是直线x=x2和椭圆曲线的两个交点，将直线方程代入椭圆方程，并将其看作是关于y的一元二次方程，由韦达定理可知y2+y'=-a1x2 - a3，即y' = -y2 - a1x2 - a3，所以条件y(P) ≠ y(-Q)等价于y1 ≠ -y2 - a1x2 - a3，首先看下P ≠ Q时λ推导过程：
假如P = Q，则可以将x2用x1取代，y2用y1取代，则上述λ公式可以转换为：
这正对于倍点计算时求切线时的求导公式，即已将倍点和点和求斜率的公式统一了起来。进而，对于椭圆曲线y2 = x3 + ax + b来说，统一公式可以简化为：
对于secp256k1来说，统一公式可进一步将a忽略。
3 源码分析
本节给出典型函数的源码，并有选择的对函数进行分析。
3.1 secp256k1_gej_double_nonzero
该函数用于执行非零点的倍点计算，第13~14行相当于计算无D版本公式中Z3=2*Y1*Z1；第15~16行求得无D版本公式中E，第17行求得T2=E2，第18行求得B，第19~21行求得8B2，第22~25行求得-2M，第26行求得X3；最后第27~32行给出了Y3。
1 static SECP256K1_INLINE void secp256k1_gej_double_nonzero(secp256k1_gej *r, const secp256k1_gej *a) {
2 /* Operations: 3 mul, 4 sqr, 0 normalize, 12 mul_int/add/negate.
3 *
4 * Note that there is an implementation described at
5 * https://hyperelliptic.org/EFD/g1p/auto-shortw-jacobian-0.html#doubling-dbl-2009-l
6 * which trades a multiply for a square, but in practice this is actually slower,
7 * mainly because it requires more normalizations.
8 */
9 secp256k1_fe t1,t2,t3,t4;
10
11 r->infinity = 0;
12
13 secp256k1_fe_mul(&r->z, &a->z, &a->y);
14 secp256k1_fe_mul_int(&r->z, 2); /* Z' = 2*Y*Z (2) */
15 secp256k1_fe_sqr(&t1, &a->x);
16 secp256k1_fe_mul_int(&t1, 3); /* T1 = 3*X^2 (3) */
17 secp256k1_fe_sqr(&t2, &t1); /* T2 = 9*X^4 (1) */
18 secp256k1_fe_sqr(&t3, &a->y);
19 secp256k1_fe_mul_int(&t3, 2); /* T3 = 2*Y^2 (2) */
20 secp256k1_fe_sqr(&t4, &t3);
21 secp256k1_fe_mul_int(&t4, 2); /* T4 = 8*Y^4 (2) */
22 secp256k1_fe_mul(&t3, &t3, &a->x); /* T3 = 2*X*Y^2 (1) */
23 r->x = t3;
24 secp256k1_fe_mul_int(&r->x, 4); /* X' = 8*X*Y^2 (4) */
25 secp256k1_fe_negate(&r->x, &r->x, 4); /* X' = -8*X*Y^2 (5) */
26 secp256k1_fe_add(&r->x, &t2); /* X' = 9*X^4 - 8*X*Y^2 (6) */
27 secp256k1_fe_negate(&t2, &t2, 1); /* T2 = -9*X^4 (2) */
28 secp256k1_fe_mul_int(&t3, 6); /* T3 = 12*X*Y^2 (6) */
29 secp256k1_fe_add(&t3, &t2); /* T3 = 12*X*Y^2 - 9*X^4 (8) */
30 secp256k1_fe_mul(&r->y, &t1, &t3); /* Y' = 36*X^3*Y^2 - 27*X^6 (1) */
31 secp256k1_fe_negate(&t2, &t4, 2); /* T2 = -8*Y^4 (3) */
32 secp256k1_fe_add(&r->y, &t2); /* Y' = 36*X^3*Y^2 - 27*X^6 - 8*Y^4 (4) */
33 }之后的secp256k1_gej_double_var函数正式调用了该函数实现倍点操作，只不过在secp256k1_gej_double_var函数中对无穷远点做了特殊处理。
3.2 secp256k1_gej_add_var
2.2节已经详细分析了公式推导相关过程，所以不再对该函数进行详细分析，请根据注释自行理解。
1 /*
2 * lambda1 = x1*z2^2
3 * lambda2 = x2*z1^2
4 * lambda3 = lambda2 - lambda1
5 * lambda4 = y1*z2^3
6 * lambda5 = y2*z1^3
7 * lambda6 = lambda5 - lambda4
8 * lambda7 = lambda1 + lambda2
9 * lambda8 = lambda4 + lambda5
10 * x3 = lambda6^2 - (2*lambda1*lambda3^2 + lambda3^3)
11 * lambda9 = lambda7*lambda3^2 - 2*x3
12 * y3 = (lambda5 - lambda4)*(lambda1*lambda3^2 - x3) - lambda4*lambda3^3
13 * z3 = z1*z2*lambda3
14 */
15 static void secp256k1_gej_add_var(secp256k1_gej *r, const secp256k1_gej *a, const secp256k1_gej *b, secp256k1_fe *rzr) {
16 /* Operations: 12 mul, 4 sqr, 2 normalize, 12 mul_int/add/negate */
17 secp256k1_fe z22, z12, u1, u2, s1, s2, h, i, i2, h2, h3, t;
18
19 if (a->infinity) {
20 *r = *b;
21 return;
22 }
23
24 if (b->infinity) {
25 if (rzr != NULL) {
26 secp256k1_fe_set_int(rzr, 1);
27 }
28 *r = *a;
29 return;
30 }
31
32 r->infinity = 0;
33 secp256k1_fe_sqr(&z22, &b->z);
34 secp256k1_fe_sqr(&z12, &a->z);
35 secp256k1_fe_mul(&u1, &a->x, &z22); // u1=lambda1
36 secp256k1_fe_mul(&u2, &b->x, &z12); // u2=lambda2
37
38 secp256k1_fe_mul(&s1, &a->y, &z22); secp256k1_fe_mul(&s1, &s1, &b->z); // s1=lambda4
39 secp256k1_fe_mul(&s2, &b->y, &z12); secp256k1_fe_mul(&s2, &s2, &a->z); // s2=lambda5
40 secp256k1_fe_negate(&h, &u1, 1); secp256k1_fe_add(&h, &u2); // lambda3
41 secp256k1_fe_negate(&i, &s1, 1); secp256k1_fe_add(&i, &s2); // i=lambda6
42 if (secp256k1_fe_normalizes_to_zero_var(&h)) {
43 if (secp256k1_fe_normalizes_to_zero_var(&i)) {
44 secp256k1_gej_double_var(r, a, rzr);
45 } else {
46 if (rzr != NULL) {
47 secp256k1_fe_set_int(rzr, 0);
48 }
49 r->infinity = 1;
50 }
51 return;
52 }
53 secp256k1_fe_sqr(&i2, &i); // i2=lambda6^2
54 secp256k1_fe_sqr(&h2, &h); // h2=lambda3^2
55 secp256k1_fe_mul(&h3, &h, &h2); // h3=lambda3^3
56 secp256k1_fe_mul(&h, &h, &b->z); // Z2*lambda3
57 if (rzr != NULL) {
58 *rzr = h;
59 }
60 secp256k1_fe_mul(&r->z, &a->z, &h); // Z1*Z2*lambda3
61 secp256k1_fe_mul(&t, &u1, &h2); // t=lambda1*lambda3^2
62 r->x = t; secp256k1_fe_mul_int(&r->x, 2); secp256k1_fe_add(&r->x, &h3); secp256k1_fe_negate(&r->x, &r->x, 3); secp256k1_fe_add(&r->x, &i2);
63 secp256k1_fe_negate(&r->y, &r->x, 5); secp256k1_fe_add(&r->y, &t); secp256k1_fe_mul(&r->y, &r->y, &i);
64 secp256k1_fe_mul(&h3, &h3, &s1); secp256k1_fe_negate(&h3, &h3, 1);
65 secp256k1_fe_add(&r->y, &h3);
66 }3.3 secp256k1_gej_add_ge
该函数用于计算射影坐标和仿射坐标相加，由2.3节可知统一公式如下：
代入x1=X1/Z12，y1=Y1/Z13，x2=X2/Z22，y2=Y2/Z23，可得：
函数源码如下：
1 static void secp256k1_gej_add_ge(secp256k1_gej *r, const secp256k1_gej *a, const secp256k1_ge *b) {
2 /* Operations: 7 mul, 5 sqr, 4 normalize, 21 mul_int/add/negate/cmov */
3 static const secp256k1_fe fe_1 = SECP256K1_FE_CONST(0, 0, 0, 0, 0, 0, 0, 1);
4 secp256k1_fe zz, u1, u2, s1, s2, t, tt, m, n, q, rr;
5 secp256k1_fe m_alt, rr_alt;
6 int infinity, degenerate;
7 if (b->infinity) {
8 *r = *a;
9 return;
10 }
11
12 /** In:
13 * Eric Brier and Marc Joye, Weierstrass Elliptic Curves and Side-Channel Attacks.
14 * In D. Naccache and P. Paillier, Eds., Public Key Cryptography, vol. 2274 of Lecture Notes in Computer Science, pages 335-345. Springer-Verlag, 2002.
15 *we find as solution for a unified addition/doubling formula:
16 * lambda = ((x1 + x2)^2 - x1 * x2 + a) / (y1 + y2), with a = 0 for secp256k1's curve equation.
17 * x3 = lambda^2 - (x1 + x2)
18 * 2*y3 = lambda * (x1 + x2 - 2 * x3) - (y1 + y2).
19 *
20 *Substituting x_i = Xi / Zi^2 and yi = Yi / Zi^3, for i=1,2,3, gives:
21 * U1 = X1*Z2^2, U2 = X2*Z1^2
22 * S1 = Y1*Z2^3, S2 = Y2*Z1^3
23 * Z = Z1*Z2
24 * T = U1+U2
25 * M = S1+S2
26 * Q = T*M^2
27 * R = T^2-U1*U2
28 * X3 = 4*(R^2-Q)
29 * Y3 = 4*(R*(3*Q-2*R^2)-M^4)
30 * Z3 = 2*M*Z
31 *(Note that the paper uses xi = Xi / Zi and yi = Yi / Zi instead.)
32 *
33 *This formula has the benefit of being the same for both addition
34 *of distinct points and doubling. However, it breaks down in the
35 *case that either point is infinity, or that y1 = -y2. We handle
36 *these cases in the following ways:
37 *
38 * - If b is infinity we simply bail by means of a VERIFY_CHECK.
39 *
40 * - If a is infinity, we detect this, and at the end of the
41 * computation replace the result (which will be meaningless,
42 * but we compute to be constant-time) with b.x : b.y : 1.
43 *
44 * - If a = -b, we have y1 = -y2, which is a degenerate case.
45 * But here the answer is infinity, so we simply set the
46 * infinity flag of the result, overriding the computed values
47 * without even needing to cmov.
48 *
49 * - If y1 = -y2 but x1 != x2, which does occur thanks to certain
50 * properties of our curve (specifically, 1 has nontrivial cube
51 * roots in our field, and the curve equation has no x coefficient)
52 * then the answer is not infinity but also not given by the above
53 * equation. In this case, we cmov in place an alternate expression
54 * for lambda. Specifically (y1 - y2)/(x1 - x2). Where both these
55 * expressions for lambda are defined, they are equal, and can be
56 * obtained from each other by multiplication by (y1 + y2)/(y1 + y2)
57 * then substitution of x^3 + 7 for y^2 (using the curve equation).
58 * For all pairs of nonzero points (a, b) at least one is defined,
59 * so this covers everything.
60 */
61
62 secp256k1_fe_sqr(&zz, &a->z); /* z = Z1^2 */
63 u1 = a->x; secp256k1_fe_normalize_weak(&u1); /* u1 = U1 = X1*Z2^2 (1) */
64 secp256k1_fe_mul(&u2, &b->x, &zz); /* u2 = U2 = X2*Z1^2 (1) */
65 s1 = a->y; secp256k1_fe_normalize_weak(&s1); /* s1 = S1 = Y1*Z2^3 (1) */
66 secp256k1_fe_mul(&s2, &b->y, &zz); /* s2 = Y2*Z1^2 (1) */
67 secp256k1_fe_mul(&s2, &s2, &a->z); /* s2 = S2 = Y2*Z1^3 (1) */
68 t = u1; secp256k1_fe_add(&t, &u2); /* t = T = U1+U2 (2) */
69 m = s1; secp256k1_fe_add(&m, &s2); /* m = M = S1+S2 (2) */
70 secp256k1_fe_sqr(&rr, &t); /* rr = T^2 (1) */
71 secp256k1_fe_negate(&m_alt, &u2, 1); /* Malt = -X2*Z1^2 */
72 secp256k1_fe_mul(&tt, &u1, &m_alt); /* tt = -U1*U2 (2) */
73 secp256k1_fe_add(&rr, &tt); /* rr = R = T^2-U1*U2 (3) */
74 /** If lambda = R/M = 0/0 we have a problem (except in the "trivial"
75 *case that Z = z1z2 = 0, and this is special-cased later on). */
76 degenerate = secp256k1_fe_normalizes_to_zero(&m) &
77 secp256k1_fe_normalizes_to_zero(&rr);
78 /* This only occurs when y1 == -y2 and x1^3 == x2^3, but x1 != x2.
79 * This means either x1 == beta*x2 or beta*x1 == x2, where beta is
80 * a nontrivial cube root of one. In either case, an alternate
81 * non-indeterminate expression for lambda is (y1 - y2)/(x1 - x2),
82 * so we set R/M equal to this. */
83 rr_alt = s1;
84 secp256k1_fe_mul_int(&rr_alt, 2); /* rr = Y1*Z2^3 - Y2*Z1^3 (2) */
85 secp256k1_fe_add(&m_alt, &u1); /* Malt = X1*Z2^2 - X2*Z1^2 */
86
87 secp256k1_fe_cmov(&rr_alt, &rr, !degenerate);
88 secp256k1_fe_cmov(&m_alt, &m, !degenerate);
89 /* Now Ralt / Malt = lambda and is guaranteed not to be 0/0.
90 * From here on out Ralt and Malt represent the numerator
91 * and denominator of lambda; R and M represent the explicit
92 * expressions x1^2 + x2^2 + x1x2 and y1 + y2. */
93 secp256k1_fe_sqr(&n, &m_alt); /* n = Malt^2 (1) */
94 secp256k1_fe_mul(&q, &n, &t); /* q = Q = T*Malt^2 (1) */
95 /* These two lines use the observation that either M == Malt or M == 0,
96 * so M^3 * Malt is either Malt^4 (which is computed by squaring), or
97 * zero (which is "computed" by cmov). So the cost is one squaring
98 * versus two multiplications. */
99 secp256k1_fe_sqr(&n, &n);
100 secp256k1_fe_cmov(&n, &m, degenerate); /* n = M^3 * Malt (2) */
101 secp256k1_fe_sqr(&t, &rr_alt); /* t = Ralt^2 (1) */
102 secp256k1_fe_mul(&r->z, &a->z, &m_alt); /* r->z = Malt*Z (1) */
103 infinity = secp256k1_fe_normalizes_to_zero(&r->z) * (1 - a->infinity);
104 secp256k1_fe_mul_int(&r->z, 2); /* r->z = Z3 = 2*Malt*Z (2) */
105 secp256k1_fe_negate(&q, &q, 1); /* q = -Q (2) */
106 secp256k1_fe_add(&t, &q); /* t = Ralt^2-Q (3) */
107 secp256k1_fe_normalize_weak(&t);
108 r->x = t; /* r->x = Ralt^2-Q (1) */
109 secp256k1_fe_mul_int(&t, 2); /* t = 2*x3 (2) */
110 secp256k1_fe_add(&t, &q); /* t = 2*x3 - Q: (4) */
111 secp256k1_fe_mul(&t, &t, &rr_alt); /* t = Ralt*(2*x3 - Q) (1) */
112 secp256k1_fe_add(&t, &n); /* t = Ralt*(2*x3 - Q) + M^3*Malt (3) */
113 secp256k1_fe_negate(&r->y, &t, 3); /* r->y = Ralt*(Q - 2x3) - M^3*Malt (4) */
114 secp256k1_fe_normalize_weak(&r->y);
115 secp256k1_fe_mul_int(&r->x, 4); /* r->x = X3 = 4*(Ralt^2-Q) */
116 secp256k1_fe_mul_int(&r->y, 4); /* r->y = Y3 = 4*Ralt*(Q - 2x3) - 4*M^3*Malt (4) */
117
118 /** In case a->infinity == 1, replace r with (b->x, b->y, 1). */
119 secp256k1_fe_cmov(&r->x, &b->x, a->infinity);
120 secp256k1_fe_cmov(&r->y, &b->y, a->infinity);
121 secp256k1_fe_cmov(&r->z, &fe_1, a->infinity);
122 r->infinity = infinity;
123 }注释中第22~30行正是由之前的公式而来，源码中参数b是一个仿射坐标的曲线点(x2, y2)，其对应雅可比射影坐标可表示为(x2, y2, 1)。源码第33~59行注释说了如何对4中特殊情况进行处理的，尤其是第4中情况degenerate case，它特指如下情形：
下图展示了一个典型的secp256k1椭圆曲线 y2=x3+7中的退化情况示例，P和Q满足以上1和3两个条件，却在实数域中不可能满足第2个条件，但是有限域Fp中该条件却有可能成立，所以代码中也把这种情况考虑进去了，其他特殊情况比较显而易见，不再详细说明。
参考：
weierstrass elliptic curves and side-channel attacks
https://blog.csdn.net/qq_50680426/article/details/120940244
来源：程序园用户自行投稿发布，如果侵权，请联系站长删除
免责声明：如果侵犯了您的权益，请联系站长，我们会及时删除侵权内容，谢谢合作！
页:
[1]