# 倍点公式
若P1 == P2，则:
这里a的值是P-3，所以为了计算方便，可以把λ1的计算公式变化一下：
我们根据公式实现代码即可：
# 算法实现
void gm_point_double(gm_point_t * r, const gm_point_t * p) {
gm_bn_t tmp1, tmp2;
gm_bn_t X3, Y3, Z3;
gm_bn_t R1, R2, R3;
if (gm_is_at_infinity(p)) {
gm_point_copy(r, p);
return;
}
// λ1 = 3X^2 + aZ^4 = 3X^2 - 3Z^4 = 3 * ( (X + Z^2) * (X - Z^2) )
// z^2
gm_bn_sqr(tmp1, p->Z, GM_BN_P);
// X + Z^2
gm_bn_add(tmp2, p->X, tmp1, GM_BN_P);
// X - Z^2
gm_bn_sub(tmp1, p->X, tmp1, GM_BN_P);
// (X + Z^2) * (X - Z^2)
gm_bn_mont_mul(R1, tmp1, tmp2, GM_BN_P);
// λ1 = 3 * (X + Z^2) * (X - Z^2)
gm_bn_add(tmp1, R1, R1, GM_BN_P);
gm_bn_add(R1, tmp1, R1, GM_BN_P);
// λ2 = X4Y^2 = X * (2Y)^2
// Z3 = 2YZ
// 2Y
gm_bn_add(tmp1, p->Y, p->Y, GM_BN_P);
// Z3 = 2YZ
gm_bn_mont_mul(Z3, tmp1, p->Z, GM_BN_P);
// λ2 = X * (2Y)^2
gm_bn_sqr(R2, tmp1, GM_BN_P);
gm_bn_mont_mul(R2, p->X, R2, GM_BN_P);
// λ3 = 8Y^4 = 2 * 4Y^4 = 2 * (2 * Y^2) ^ 2
// Y^2
gm_bn_sqr(tmp1, p->Y, GM_BN_P);
// 2Y^2
gm_bn_add(tmp1, tmp1, tmp1, GM_BN_P);
// (2 * Y^2) ^ 2
gm_bn_sqr(tmp1, tmp1, GM_BN_P);
// λ3 = 2 * (2 * Y^2) ^ 2
gm_bn_add(R3, tmp1, tmp1, GM_BN_P);
// X3 = λ1^2 − 2λ2
// R1^2
gm_bn_sqr(tmp1, R1, GM_BN_P);
// 2R2
gm_bn_add(tmp2, R2, R2, GM_BN_P);
// X3 = λ1^2 − 2λ2
gm_bn_sub(X3, tmp1, tmp2, GM_BN_P);
// Y3 = λ1(λ2 − X3) − λ3
// λ2 − X3
gm_bn_sub(tmp1, R2, X3, GM_BN_P);
// λ1(λ2 − X3)
gm_bn_mont_mul(tmp2, R1, tmp1, GM_BN_P);
// Y3 = λ1(λ2 − X3) − λ3
gm_bn_sub(Y3, tmp2, R3, GM_BN_P);
// output
gm_bn_copy(r->X, X3);
gm_bn_copy(r->Y, Y3);
gm_bn_copy(r->Z, Z3);
}
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
# 单元测试
/**
* 点单元测试
* @param pa 点a
* @param pb 点b
* @param res 预期结果
* @param alg 测试算法dbl、add、mul等
*/
void test_ec(const char * pa, const char * pb, const char * res, const char * alg) {
int i, j;
gm_point_t p1, p2, r;
char e_r_hex[129] = {0};
gm_bn_t k;
gm_point_from_hex(&p1, pa);
gm_point_from_hex(&p2, pb);
gm_point_from_hex(&r, pa);
if(strcmp(alg, "point_dbl") == 0){ // 10万
int i = 0;
for (i = 0; i < 100000; i++) {
gm_point_double(&r, &r);
}
}
gm_point_to_hex(&r, e_r_hex);
printf("r = %s\n", e_r_hex);
printf("test result: %s\n", (strcmp(res, e_r_hex) == 0 ? "ok" : "fail"));
}
#define TEST_EC_ALG(alg_name, a, b, r) \
do { \
if(strncmp(argv[1], alg_name, strlen(alg_name)) == 0) { \
test_ec(a, b, r, alg_name); \
} \
}while(0)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
main函数增加
TEST_EC_ALG("point_dbl",
"EB04AAE0D53FBA1E3611D5B9ED6EFA3EE5BA57C41AA7A09DDC5816AF09057757CE6FA0678392F4716E45F58E7322C76D5997B1FE44C36D8A5A59B146EE162B93",
"00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000",
"976A325416F812A692E6880A72287604C32FE61839EA4732040451E831E1128F6654AC40E0754A30095AC2FC4C49D4147D1E2395190AEFFB555BD55AC736B2BA");
1
2
3
4
2
3
4
算法效率：
saintdeMacBook-Pro:bn saint$ time ./a.out point_dbl
r = 976A325416F812A692E6880A72287604C32FE61839EA4732040451E831E1128F6654AC40E0754A30095AC2FC4C49D4147D1E2395190AEFFB555BD55AC736B2BA
test result: ok
real 0m1.232s
user 0m1.164s
sys 0m0.013s
1
2
3
4
5
6
7
2
3
4
5
6
7
未经本人同意，禁止转载！