大数模乘算法选择

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

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

分治乘法

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

safedead的博文

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

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

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

蒙哥马利算法

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

提示

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

本文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

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

$0x08 * 0x01020304 + 0x07 * 0x01020304 * 2^8 + 0x06 * 0x01020304 * 2^{16} + 0x05 * 0x01020304 * 2^{24}$

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

      FFFFFFFF
x           FF
=   FEFFFFFF01

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

$$\frac{FEFFFFFF01}{2^8}$$

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

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

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

$$\frac{a * b}{2^{32}}\pmod {p} \equiv a * b * 2^{-32}\pmod {p}$$

这里有两个问题：

1. 每轮乘法的结果并不一定是$2^8$的倍数，除法就有困难
2. 最终的结果是$a * b * 2^{-32}\pmod {p}$，并不是$a * b \pmod {p}$

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

每轮结果处理

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

$$1 \pmod {10} \equiv 1+10 \pmod {10} \equiv 1+q*p \pmod {p}$$

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

进一步推导$a * b * 2^{-32}\pmod {p}$这个计算过程中，无论每一轮我们加上多少个q * p，结果都是一样的。

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

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

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

中途总结

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

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

每轮结果用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;
}

计算当前轮的乘法结果，加上上一轮的乘法结果，再加上q * p，再除以$2^8$。循环4轮，最后再判断是否大于p，大于再减掉一个p就是最终$a * b * 2^{-32}\pmod {p}$的结果

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

q的计算

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

$q = (a[0]*b[i] + d_0) * t \pmod {2^8}$

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

理论中t的值为： $t = -p^{-1}\pmod {2^8}$

伪代码为：

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

由此可知，$t < 2^8$

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

蒙哥马利域

蒙哥马利乘法定义即是：$a * b * 2^{-32}\pmod {p}$

上述可知，最终结果为$a * b * 2^{-32}\pmod {p}$，并不是我们想要的$a * b \pmod {p}$

那是不是有这样一个变换，比如将输入参数变换一下，输出参数逆变换一下，就可以得到我们想要的$a * b \pmod {p}$呢？

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

$$a′ = a * 2^{32}\pmod {p}$$

$$b′ = b * 2^{32}\pmod {p}$$

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

$$a * 2^{32} * b * 2^{32} * 2^{-32}\pmod {p} \equiv a * b * 2^{32} \pmod {p}$$

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

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

$$(n * 2^{32}) * (1) * 2^{-32}\pmod {p} \equiv n \pmod {p}$$

那么如何将大数n转换到蒙哥马利域呢，也比较简单，乘以$2^{64} \pmod {p}$即可：

$$(n) * (2^{64}) * 2^{-32}\pmod {p} \equiv n * 2^{32} \pmod {p}$$

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

蒙哥马利公式

针对256-bit的大数，r的取值应该是：$2^{256}$

乘法公式：

$a * b * r^{-1}\pmod {p} \equiv a * b * 2^{-256}\pmod {p}$

转换到蒙哥马利域公式：

$(a) * (r^2) * r^{-1}\pmod {p} \equiv a * r \pmod {p} \equiv a * 2^{256}\pmod {p}$

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

$(a * r) * (1) * r^{-1}\pmod {p} \equiv a \pmod {p}$

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

蒙哥马利域转换算法

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

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

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

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

从上述可知，r的取值是$2^{256}$，$r^2$应该是$2^{512}$，那么将a转到蒙哥马利域时$r^2\pmod p$的值应该是什么呢？

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

$$GM\_BN\_MONT\_PRR = 2^{512} \pmod p$$

及

$$GM\_BN\_MONT\_NRR = 2^{512} \pmod n$$

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

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

上一页大数模减

下一页大数蒙哥马利模乘