铜锁 SM2 算法性能优化实践(二)|快速模约减算法实现

原创
08/02 16:00
阅读数 114

本文是对 《铜锁 SM2 算法性能优化实践(一)|综述》 的深入扩展,详细介绍快速模约减算法实现
图片可能无法显示,欢迎阅读原文链接:https://www.yuque.com/tsdoc/ts/pacggl7s29yzsm19

背景

在 SM2 数字签名算法中,有限域运算是椭圆曲线点运算和数字签名运算的基础,其运算性能将直接影响整体性能表现。而在有限域运算中,模约减运算的整体使用频率最高,在模乘、模平方、模逆等运算中均有运用。以模乘法为例,其定义为:对于素数域参数 p 和乘数 a,模乘运算的结果等于 a*a mod p。也就是说,对于乘法结果大于 p 的中间值,模乘法需要将其约减到 p 以内,此过程因此被称为模约减

由于 p 是一个非常大的素数,且上述运算中间值最大可达 p2-1,采用除法来计算模约减的效率将极其低下 在具体实践中,有限域除法一般借助乘法逆元实现,其耗时大约是模乘运算的 200-300 倍) 。目前,铜锁项目中 SM2 算法使用了通用蒙哥马利约分算法实现模约减过程,暂时没有针对曲线特化的快速模约减算法实现。相应的,部分常用的 NIST 曲线 nistp224nistp256 和 nistp521 基于广域梅森素数实现了曲线特化的快速模约减算法,并运用于各自的 64-bit 平台优化中。

考虑到 SM2 曲线同样采用广域梅森素数作为素数域参数,且实现 64-bit 平台优化也需要使用模约减算法,因此为 SM2 曲线实现曲线参数特化的快速模约减算法是非常有必要的。同时,为了保证对部分 32 位平台的兼容性和后续扩展,在功能开发的过程中保留了原有的条件编译模式,同时支持 32 位和 64 位平台使用该快速模约减算法。测试结果表明,使用了快速模约减的 SM2 模乘运算与简单模乘和蒙哥马利模乘相比,时间开销分别降低 81.34% 和 29.91%

预备知识

在数论中,梅森素数 (Mersenne Prime) 是形如 2n-1 的质数,其中 n 为整数。梅森素数的数量非常稀少,直到今天 (2023 年 8 月) 仅有 51 个梅森素数被发现。由于梅森素数的构成与 2 的指数幂相关,在以二进制为数据表示方式的现代计算机中,它不仅便于计算和存储,也能够加速有限域运算,因此在密码学领域广受青睐。NIST 曲线家族中,nistp521 曲线的素数域参数 p 就选择了梅森素数 2521-1。

然而,梅森素数的数量实在是过于稀少,以至于在椭圆曲线素数域参数 p 的合理取值范围内仅有一个梅森素数 2521-1,其他的梅森素数要么过大,要么太小,不满足使用要求。因此,许多椭圆曲线的素数域参数 p 选用了梅森素数的拓展版本——广义梅森素数 (General Mersenne Prime) 。广义梅森素数是形如 f(2n) 的素数,其中 f(x) 是具有较小整数系数的低次多项式。例如,nistp256 曲线采用如下所示的广义梅森素数作为素数域参数 p:

算法原理

在国标 GM/T 0003.5-2012 中,给出了 SM2 曲线的推荐素数域参数 p。p 是一个广义梅森素数,表示如下:

下面以大整数 v(v<p2<2512) 和参数 psm2p256 为例,推导基于广义梅森素数的快速模约减算法。该算法利用参数 psm2p256 各项系数之间最小差为 232 的特点,将表示中间值的大整数按照 32 位分割,并转化为如下表示形式:

然后将 v[8]-v[15] 的高位参数约化到 v[0]-v[7]  的低位参数上。在约化之前,我们先根据将参数 psm2p256 转化为按照 32 位分割的表示形式:

接下来处理所有高位参数,此处以 v[8] 为例,设 v[8]=a,那么 v[8] 可以表示为:

注意到,根据模运算的运算规则:

也就是说,在原数的基础上增减模数 p 的倍数并不影响模运算的结果。那么,为了消去高位参数 a,不妨减去一个 a*p,此时有:

上下相减,高位参数 v[8] 被成功消除,低位约化的结果为:

依次对 v[8]-v[15] 做低位约化即可消除所有高位参数。将所有的低位约化结果合并同类项,可得到如下的 SM2 快速模约减算法。该算法由白国强[1]等人于 2014 年首次提出,实际运用中也有多种不同的参数组合形式:

基于参数 psm2p256 的快速模约减算法

最后对约化后的低位参数 v[0]-v[7] 再做一次约减,旨在将模约减结果规约到 [0,p)范围内。 此处最多需要减去 13p 即可得到最终结果,无需使用成本昂贵的模除法。从图1中不难看出,基于广义梅森素数的 SM2 快速模约减算法仅需要少量的有限域加法和减法即可实现,性能表现明显优于使用模逆元的模约减算法。

算法实现

本节将介绍所实现算法的核心部分,完整算法实现请参见 crypto/bn/bn_sm2.c

前文提到,在高位参数约化后,还需要对低位参数做约化以获得最终结果,也就是减去参数 p 的某个倍数 np(0≤p≤13)。为了进一步提高性能,在代码实现中我们构造一个预计算表以存储 1-13p,以节约计算 np 的开销。又因为约减结果在 【0,p) 范围内,且 p2256 ,因此只需要保留预计算参数的低 256 位:

/*
 * /* Pre-computed tables are "carry-less" values of modulus*(i+1),
 * all values are in little-endian format.
 */
static const BN_ULONG _sm2_p_256[][BN_SM2_256_TOP] = {
    {0xFFFFFFFFFFFFFFFFull, 0xFFFFFFFF00000000ull,
     0xFFFFFFFFFFFFFFFFull, 0xFFFFFFFEFFFFFFFFull},/*p*/
    {0xFFFFFFFFFFFFFFFEull, 0xFFFFFFFE00000001ull,
     0xFFFFFFFFFFFFFFFFull, 0xFFFFFFFDFFFFFFFFull},/*2p*/
    ...
    {0xFFFFFFFFFFFFFFF4ull, 0xFFFFFFF40000000Bull,
     0xFFFFFFFFFFFFFFFFull, 0xFFFFFFF3FFFFFFFFull},/*12p*/
    {0xFFFFFFFFFFFFFFF3ull, 0xFFFFFFF30000000Cull,
     0xFFFFFFFFFFFFFFFFull, 0xFFFFFFF2FFFFFFFFull}/*13p*/
};

同时,为了保证被约减数 v 满足 0≤v<p2 ,在模约减前还需要对输入进行验证。基于同样的原因,我们还需要预计算 p2 的值:

/* pre-compute the value of p^2 check if the input satisfies input < p^2. */
static const BN_ULONG _sm2_p_256_sqr[] = {
    0x0000000000000001ULL, 0x00000001FFFFFFFEULL,
    0xFFFFFFFE00000001ULL, 0x0000000200000000ULL,
    0xFFFFFFFDFFFFFFFEULL, 0xFFFFFFFE00000003ULL,
    0xFFFFFFFFFFFFFFFFULL, 0xFFFFFFFE00000000ULL
};

在 64 位平台和部分支持 64 位无符号整型的 32 位平台上,由于 v[i] 满足 0≤v[i]<232,因此可以使用 64 位无符号整形变量 accumulator 作为中间值存储累加结果,从而在 v[0]-v[7] 参数上直接累加高位参数约化结果。以 v[0] 为例,其相应代码如下:

    	SM2_INT64 acc;         /* accumulator */

        acc = rp[0];
        acc += bp[8 - 8];
        acc += bp[9 - 8];
        acc += bp[10 - 8];
        acc += bp[11 - 8];
        acc += bp[12 - 8];
        acc += bp[13 - 8];
        acc += bp[14 - 8];
        acc += bp[15 - 8];
        acc += bp[13 - 8];
        acc += bp[14 - 8];
        acc += bp[15 - 8];	/*累加高位参数约化结果*/
        rp[0] = (unsigned int)acc; /*获得结果的低32位*/
        acc >>= 32; /*进位部分右移32位,作为v[1]输入的一部分*/

对于部分不支持 64 位无符号整型的 32 位平台,32 位整形变量由于数据溢出无法作为中间值保留进位结果,上述代码将不能正确运行。此时需要借助现有工具函数和铜锁底层大数运算函数的支持,按照前文展示的算法逐一累加 s1~s14,每次还需要记录进位 carry

        /*
         * r_d += s2 + s6 + s7 + s8 + s9
         * s2 = (c15, c14, c13, c12, c11, 0, c9, c8)
         */
        sm2_set_256(t_d, buf.bn, 15, 14, 13, 12, 11, 0, 9, 8);
        carry += (int)bn_add_words(r_d, r_d, t_d, BN_SM2_256_TOP);
        /*
         * s6 = (c11, c11, c10, c15, c14, 0, c13, c12)
         */
        sm2_set_256(t_d, buf.bn, 11, 11, 10, 15, 14, 0, 13, 12);
        carry += (int)bn_add_words(r_d, r_d, t_d, BN_SM2_256_TOP);
        /*
         * s7 = (c10, c15, c14, c13, c12, 0, c11, c10)
         */
        sm2_set_256(t_d, buf.bn, 10, 15, 14, 13, 12, 0, 11, 10);
        carry += (int)bn_add_words(r_d, r_d, t_d, BN_SM2_256_TOP);
        /*
         * s8 = (c9, 0, 0, c9, c8, 0, c10, c9)
         */
        sm2_set_256(t_d, buf.bn, 9, 0, 0, 9, 8, 0, 10, 9);
        carry += (int)bn_add_words(r_d, r_d, t_d, BN_SM2_256_TOP);
        /*
         * s9 = (c8, 0, 0, 0, c15, 0, c12, c11)
         */
        sm2_set_256(t_d, buf.bn, 8, 0, 0, 0, 15, 0, 12, 11);
        carry += (int)bn_add_words(r_d, r_d, t_d, BN_SM2_256_TOP);

最后处理进位,得到最终结果,此处变量 mask 的作用是减少代码中不必要的分支,尽可能提高其抗侧信道攻击强度:

    mask =
        0 - (PTR_SIZE_INT) (*u.f) (c_d, r_d, _sm2_p_256[0], BN_SM2_256_TOP);
    mask &= 0 - (PTR_SIZE_INT) carry;
    res = c_d;
    res = (BN_ULONG *)(((PTR_SIZE_INT) res & ~mask) |
                       ((PTR_SIZE_INT) r_d & mask));
    sm2_cp_bn(r_d, res, BN_SM2_256_TOP);
    r->top = BN_SM2_256_TOP;
    bn_correct_top(r);

测试

正确性测试

为了验证快速模约减算法的正确性,我们在测试环节设计了相应的单元测试。该测试的核心思路是,随机生成一批有限域大整数 ai,bi,分别使用铜锁提供的默认方法和本文实现的快速模约减方法计算 ai*bi mod p,然后比较二者计算结果是否相等:

    /*
     * Randomly generate two numbers in the prime field and multiply them,then compare the results of fast modular reduction and conventional algorithms.
     */
    for (i = 0; i < 100; i++){
        if (!TEST_true(BN_priv_rand_range_ex(a, p, 0, ctx))
            || !TEST_true(BN_priv_rand_range_ex(b, p, 0, ctx))
            || !TEST_true(BN_mul(ab_fast, a, b, ctx))
            || !TEST_true(BN_sm2_mod_256(ab_fast, ab_fast, p, ctx))
            || !TEST_true(BN_mod_mul(ab, a, b, p, ctx))) {
            goto done;
        }

        if (!TEST_int_eq(BN_ucmp(ab, ab_fast), 0)){
            goto done;
        }
    }

可以使用如下命令调用测试:

make test TESTS="test_mod_sm2"

测试结果:

➜  Tongsuo git:(sm2p256) ✗ make test TESTS="test_mod_sm2"
/Library/Developer/CommandLineTools/usr/bin/make depend && /Library/Developer/CommandLineTools/usr/bin/make _tests
10-test_mod_sm2.t .. ok   
All tests successful.
Files=1, Tests=1,  1 wallclock secs ( 0.00 usr  0.01 sys +  0.07 cusr  0.02 csys =  0.10 CPU)
Result: PASS

性能测试

为了充分比较快速模约减算法的性能,我们选择比较 SM2 有限域运算中较常用的模乘法运算作为测试单元,测试使用了 SM2 快速模约减函数 BN_sm2_mod_256、简单模约减函数 BN_nnmod 的模乘函数 ossl_ec_GFp_field_mul 和通用蒙哥马利模乘函数 BN_mod_mul_montgomery 的模乘法运算耗时。测试在 Apple M2 pro 平台上进行,测试数据如下表所示:

测试结果表明,使用 SM2 快速模约减函数的模乘法耗时最少,相较于简单模约减能够减少 81.34% 的时间开销,即使对使用了汇编优化的蒙哥马利模乘运算,在计算时间上仍有 29.91% 的优势。

下篇 SM2 算法性能优化实践系列文章我们将介绍“快速模逆元算法实现”,敬请期待!

文中链接

[1]白国强等人提出的快速模约减算法:https://ieeexplore.ieee.org/document/7011249

展开阅读全文
加载中
点击引领话题📣 发布并加入讨论🔥
打赏
0 评论
0 收藏
0
分享
返回顶部
顶部