忆信-个人博客

外观

【零依赖实现国密算法C语言版】国密SM2之倍点

约 579 字大约 2 分钟

2023-05-13

倍点公式

若P1 == P2，则:

λ1=3x2+az4,λ2=4xy2,λ3=8y4,x3=λ12−2λ2,y3=λ1(λ2−x3)−λ3,z3=2yzλ1 = 3x^2 + az^4, λ2 = 4xy^2, λ3 = 8y^4, x3 = λ1^2 − 2λ2, y3 = λ1(λ2 − x3) − λ3, z3 = 2yz

这里a的值是P-3，所以为了计算方便，可以把λ1的计算公式变化一下：

λ1=3X2+aZ4=3X2−3Z4=3∗((X+Z2)∗(X−Z2))λ1 = 3X^2 + aZ^4 = 3X^2 - 3Z^4 = 3 * ( (X + Z^2) * (X - Z^2) )

我们根据公式实现代码即可：

算法实现

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);
}

单元测试

/**
 * 点单元测试
 * @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)

main函数增加

TEST_EC_ALG("point_dbl",
            "EB04AAE0D53FBA1E3611D5B9ED6EFA3EE5BA57C41AA7A09DDC5816AF09057757CE6FA0678392F4716E45F58E7322C76D5997B1FE44C36D8A5A59B146EE162B93",
            "00000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000",
            "976A325416F812A692E6880A72287604C32FE61839EA4732040451E831E1128F6654AC40E0754A30095AC2FC4C49D4147D1E2395190AEFFB555BD55AC736B2BA");

算法效率：

saintdeMacBook-Pro:bn saint$ time ./a.out point_dbl
r = 976A325416F812A692E6880A72287604C32FE61839EA4732040451E831E1128F6654AC40E0754A30095AC2FC4C49D4147D1E2395190AEFFB555BD55AC736B2BA
test result: ok

real	0m1.232s
user	0m1.164s
sys	0m0.013s