铜锁 SM2 算法性能优化实践（二）｜快速模约减算法实现
本文是对 《铜锁 SM2 算法性能优化实践（一）｜综述》 的深入扩展，详细介绍快速模约减算法实现。
原文链接：https://www.yuque.com/tsdoc/ts/pacggl7s29yzsm19
背景
在 SM2 数字签名算法中，有限域运算是椭圆曲线点运算和数字签名运算的基础，其运算性能将直接影响整体性能表现。而在有限域运算中，模约减运算的整体使用频率最高，在模乘、模平方、模逆等运算中均有运用。以模乘法为例，其定义为：对于素数域参数 p 和乘数 a，模乘运算的结果等于 a*a mod p。也就是说，对于乘法结果大于 p 的中间值，模乘法需要将其约减到 p 以内，此过程因此被称为模约减。
由于 p 是一个非常大的素数，且上述运算中间值最大可达 p2-1，采用除法来计算模约减的效率将极其低下 （ 在具体实践中，有限域除法一般借助乘法逆元实现，其耗时大约是模乘运算的 200-300 倍） 。目前，铜锁项目中 SM2 算法使用了通用蒙哥马利约分算法实现模约减过程，暂时没有针对曲线特化的快速模约减算法实现。相应的，部分常用的 NIST 曲线 nistp224、nistp256 和 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 年首次提出，实际运用中也有多种不同的参数组合形式：
图 1：基于参数 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) 范围内，且
，因此只需要保留预计算参数的低 256 位：
同时，为了保证被约减数 v 满足 0≤v<p2 ，在模约减前还需要对输入进行验证。基于同样的原因，我们还需要预计算 p2 的值：
在 64 位平台和部分支持 64 位无符号整型的 32 位平台上，由于 v[i] 满足 0≤v[i]<232，因此可以使用 64 位无符号整形变量 accumulator 作为中间值存储累加结果，从而在 v[0]-v[7] 参数上直接累加高位参数约化结果。以 v[0] 为例，其相应代码如下：
对于部分不支持 64 位无符号整型的 32 位平台，32 位整形变量由于数据溢出无法作为中间值保留进位结果，上述代码将不能正确运行。此时需要借助现有工具函数和铜锁底层大数运算函数的支持，按照前文展示的算法逐一累加 s1~s14，每次还需要记录进位 carry：
最后处理进位，得到最终结果，此处变量 mask 的作用是减少代码中不必要的分支，尽可能提高其抗侧信道攻击强度：
测试
正确性测试
为了验证快速模约减算法的正确性，我们在测试环节设计了相应的单元测试。该测试的核心思路是，随机生成一批有限域大整数 ai，bi，分别使用铜锁提供的默认方法和本文实现的快速模约减方法计算 ai*bi mod p，然后比较二者计算结果是否相等：
可以使用如下命令调用测试：
测试结果：
性能测试
为了充分比较快速模约减算法的性能，我们选择比较 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
版权声明: 本文为 InfoQ 作者【铜锁开源密码库】的原创文章。
原文链接:【http://xie.infoq.cn/article/0b1438f87f5b400474dffb13c】。文章转载请联系作者。
还未添加个人签名 2023-04-24 加入
铜锁/Tongsuo是一个提供现代密码学算法和安全通信协议的开源基础密码库。
评论