【零依赖实现国密算法C语言版】国密SM2之蒙哥马利算法

# 大数模乘算法选择

大数模乘算法目前我知道两种，一种是分治乘法后再约减，一种是边乘边约减（蒙哥马利算法）

对于算法的选择，经过比较后，发现还是蒙哥马利算法效率更高一些，当然指的是汇编优化情况下。

# 分治乘法

分治乘法，大家可以参考safedead大牛的博文：

safedead的博文

用x64汇编语言编写384位无符号整数乘法(上) (opens new window)
国密SM2素域椭圆曲线快速约减算法x64编程研究(上) (opens new window)
使用广义梅森素数的有限域快速约减求模算法推导实例 (opens new window)

大家看理论部分就可以，当然使用汇编实现SM2算法的，也可以参考一下大牛的汇编代码

gmssl-v3-dev模乘，对参数p的模乘就是计算完乘法结果后，用这里提到的快速约减算法进行约减。

# 蒙哥马利算法

概念比较复杂，我自己也没完全理解，我就不仔细讲原理了，因为我也讲不清，这里推荐两篇博文给大家

TIP

蒙哥马利算法详解 (opens new window)
如何高效进行模乘、模幂运算？——蒙哥马利算法（Montgomery Algorithm）从入门到精通 (opens new window)

本文a * b mod p以及a * b mod n使用的都是蒙哥马利算法。相关理论需要大家自己去啃。

# 蒙哥马利模乘

我们知道，两个256-bit的大数相乘，结果为512-bit的大数，举个1字节能表示最大值相乘的例子

0xFF * 0xFF = 0xFE01

它的结果是不会超过2字节的。

对于一个512-bit的大数，要计算mod p或mod n是不好计算的，所以就要想，有没有办法边计算乘法结果，边取模，蒙哥马利算法就是这个思路。

我们先拿4字节数来举例，后面再扩展到256bit大数

   0x01020304
x  0x05060708

1
2

按小学数学，结果应该是：

所以需要经过4轮计算，将每轮结果相加就是乘法的结果了，再假设：

      FFFFFFFF
x           FF
=   FEFFFFFF01

1
2
3

这一轮中，都是数值中的最大值，乘法结果为FEFFFFFF01，结果不会超过5字节，如何把它变成4字节呢？

除以即可，这里都是最大值FF来举例的目的，就是让大家清楚，会不会溢出，应该用多大的内存才存的下乘法结果。

也就是说，只要每一轮计算结果都除以，最终的结果也会是4字节，计算4字节mod n或mod p就很简单了（这里暂时假设p和n也就4字节的素数）。

那么这样a * b每轮都除以计算乘法结果，最后再对p取模，结果应该是:

这里有两个问题：

每轮乘法的结果并不一定是的倍数，除法就有困难
最终的结果是，并不是

接下来解决这两个问题即可

# 每轮结果处理

先不拿素数来举例，假设p = 10，是不是可以推导：

也就是说1 % 10 == 11 %10 == (1+q*p) % p，这几个数对10取模的结果都是一样的，因此是不是可以推导，在上述每轮的结果中，我们加上的q * p，只要是p的倍数，那么最终结果对p取模都是一样的。

进一步推导这个计算过程中，无论每一轮我们加上多少个q * p，结果都是一样的。

再进一步推导，只要我们想办法让每一轮加上q * p正好是的倍数，不就可以解决上述提到的第1个问题。

也就是我们计算的重心转移到了计算q这个值去了。

计算每轮结果除以已经很简单了，比如20是10的倍数，它除以10只要把低位抛弃掉即可，同理计算每一轮除以只要把低位抛弃掉即可。

# 中途总结

下面来总结一下，写个伪代码，假设有4字节数a * b，对p取模，我们用uint8_t来表示这个数（忽略中间结果溢出问题先）：

uint8_t a[4]
uint8_t b[4]
uint8_t d[5]

1
2
3

每轮结果用d来存放，回顾一下d用5个字节就可以了吧

d = 0;
for(i = 0; i < 4; i++) {
	d = ((b[i] * a) + d + q * p) * 2^-8;
}
if(d > p) {
	d = d - p;
}

1
2
3
4
5
6
7

计算当前轮的乘法结果，加上上一轮的乘法结果，再加上q * p，再除以。循环4轮，最后再判断是否大于p，大于再减掉一个p就是最终的结果

这里大家可能需要结合理论细品一下

# q的计算

从如何高效进行模乘、模幂运算？——蒙哥马利算法（Montgomery Algorithm）从入门到精通 (opens new window)这篇文章可知：

其中，为d的低2^8位，其实就是d[0]，t如何计算呢

理论中t的值为：

伪代码为：

t = 1;
for(i = 1; i < 8; i++) {
	t = t * t * p mod 2^8;
}
t = 2^8 - t;

1
2
3
4
5

由此可知，

接下来我们解决第2个问题，即乘法结果不是我们想要的结果问题

# 蒙哥马利域

蒙哥马利乘法定义即是：

上述可知，最终结果为，并不是我们想要的

那是不是有这样一个变换，比如将输入参数变换一下，输出参数逆变换一下，就可以得到我们想要的呢？

有的，这个输入变换就是，将每个参数都先乘以，再对p取模，即：

将参数带入蒙哥马利乘法公式，结果为：

设大数n，定义为大数n的蒙哥马利域，从上面的推导我们可知，蒙哥马利域的大数使用蒙哥马利乘法，结果还是属于蒙哥马利域。

蒙哥马利域的大数，如何转换为普通的大数n呢？有了蒙哥马利乘法，其实很简单，乘以1即可:

那么如何将大数n转换到蒙哥马利域呢，也比较简单，乘以即可：

接下来，我们回归理论知识，并将它扩展到256-bit大数

# 蒙哥马利公式

针对256-bit的大数，r的取值应该是：

乘法公式：

转换到蒙哥马利域公式：

蒙哥马利域转换到普通大数公式：

256-bit大数乘法，如果我们还是用4轮来实现，那么每一轮应该要除以

# 蒙哥马利域转换算法

假设我们已经实现了蒙哥马利乘法

// r = a * b * 2^{-256} mod m
void gm_bn_mont_mul(gm_bn_t r, const gm_bn_t a, const gm_bn_t b, const gm_bn_t m);

1
2

接下来实现相应的转换函数就简单了，先上代码：

// 将a转到蒙哥马利域 mod m
void gm_bn_to_mont(gm_bn_t r, const gm_bn_t a, const gm_bn_t m) {
	// 乘以r的平方即可
    gm_bn_mont_mul(r, a, (m[0] == GM_BN_N[0] ? GM_BN_MONT_NRR : GM_BN_MONT_PRR), m);
}

// 将a由蒙哥马利域转为普通大数 mod m
void gm_bn_from_mont(gm_bn_t r, const gm_bn_t a, const gm_bn_t m) {
	// 乘以1即可
    gm_bn_mont_mul(r, a, GM_BN_ONE, m);
}

1
2
3
4
5
6
7
8
9
10
11

从上述可知，r的取值是，应该是，那么将a转到蒙哥马利域时的值应该是什么呢？

这个本文未给出计算方式，本文是借助openssl计算的，即计算：

及

下一章节我们将来实现蒙哥马利模乘

未经本人同意，禁止转载！

← 大数模减大数蒙哥马利模乘 →