Skip to main content

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

本文旨在介绍本系列文章第一篇《综述》中涉及的快速模约减算法的详细推导和实现过程

背景

在SM2数字签名算法中,有限域运算是椭圆曲线点运算和数字签名运算的基础,其运算性能将直接影响整体性能表现。而在有限域运算中,模约减运算的整体使用频率最高,在模乘、模平方、模逆等运算中均有运用。以模乘法为例,其定义为:对于素数域参数pp和乘数aa,模乘运算的结果等于 aamodpa * a \mod p。也就是说,对于乘法结果大于pp中间值,模乘法需要将其约减到pp以内,此过程因此被称为模约减
由于pp是一个非常大的素数,且上述运算中间值最大可达p21p^2-1,采用除法来计算模约减的效率将极其低下(在具体实践中,有限域除法一般借助乘法逆元实现,其耗时大约是模乘运算的200-300倍)。目前,铜锁项目中SM2算法使用了通用蒙哥马利约分算法实现模约减过程,暂时没有针对曲线特化的快速模约减算法实现。相应的,部分常用的NIST曲线nistp224nistp256nistp521基于广域梅森素数实现了曲线特化的快速模约减算法,并运用于各自的64-bit平台优化中。

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

预备知识

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

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

pnistp256=22562224+2192+2961p_{nistp256} = 2^{256} − 2^{224} + 2^{192} + 2^{96} − 1

算法原理

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

psm2p256=22562224296+2641p_{sm2p256} = 2^{256} − 2^{224} - 2^{96} + 2^{64} - 1

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

(v[0]20+v[1]232+v[2]264+...+v[15]2480)modp(v[0]*2^0 + v[1]*2^{32} + v[2]*2^{64} + ... + v[15]*2^{480}) \mod p, 其中 0v[i]<2320 \leq v[i]< 2^{32}

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

2^2882^2562^2242^1922^1602^1282^0962^0642^0322^000
1-1000-110-1

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

2^2882^2562^2242^1922^1602^1282^0962^0642^0322^000
0a00000000

注意到,根据模运算的运算规则,a2256modp(a2256xp)modpa * 2^{256} \mod p \equiv (a*2^{256} - x*p) \mod p,也就是说,在原数的基础上增减模数 pp 的倍数并不影响模运算的结果。那么,为了消去高位参数aa,不妨减去一个apa *p,此时有:

2^2882^2562^2242^1922^1602^1282^0962^0642^0322^000
0a00000000
a-a000-aa0-a

上下相减,高位参数v[8]v[8]被成功消除,低位约化的结果为:a2224+a296a264+aa * 2^{224} + a * 2^{96} - a * 2^{64} + a。依次对v[8]v[15]v[8] - v[15]做低位约化即可消除所有高位参数。将所有的低位约化结果合并同类项,可得到如下的sm2快速模约减算法。该算法由白国强等人于2014年首次提出,实际运用中也有多种不同的参数组合形式: image.png 图1:基于参数 psm2p256p_{sm2p256} 的快速模约减算法 最后对约化后的低位参数v[0]v[7]v[0] - v[7]再做一次约减,旨在将模约减结果规约到[0,p)[0,p)范围内。此处最多需要减去13p13p即可得到最终结果,无需使用成本昂贵的模除法。从图1中不难看出,基于广义梅森素数的SM2快速模约减算法仅需要少量的有限域加法和减法即可实现,性能表现明显优于使用模逆元的模约减算法。

算法实现

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

前文提到,在高位参数约化后,还需要对低位参数做约化以获得最终结果,也就是减去参数pp的某个倍数np(0<=p<=13)np \quad (0<=p<=13)。为了进一步提高性能,在代码实现中我们构造一个预计算表以存储113p1-13p,以节约计算npnp的开销。又因为约减结果在[0,p)[0,p)范围内,且p<2256p < 2^{256},因此只需要保留预计算参数的低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*/
};

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

/* 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]v[i]满足0v[i]<2320 \leq v[i]< 2^{32},因此可以使用64位无符号整形变量accumulator作为中间值存储累加结果,从而在v[0]v[7]v[0] - v[7]各参数上直接累加高位参数约化结果。以v[0]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 s14s_1 ~ s_{14},每次还需要记录进位carrycarry:

        /*
* 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,bia_i, b_i,分别使用铜锁提供的默认方法和本文实现的快速模约减方法计算aibimodpa_i * b_i \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快速模约减简单模约减蒙哥马利模乘约减
模乘法耗时(ns/iter)35.53190.4550.69

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