# 大数乘法
上一章节已简单的介绍了蒙哥马利相关算法,接上一章可知,我们将256-bit的乘法拆成了四轮,即每轮是
接下来我们先定义一个乘法函数:
// 乘法结果 r = k * b
static void gm_i_bn_mul(uint64_t * r, const uint64_t * k, const gm_bn_t b);
2
其中k是一个
b7b6b5b4b3b2b1b0
x k1k0
结果(数组下标从0开始):k0b0, k0b1+k1b0, .. ,k0b7+k1b6, k1b7
2
3
4
将这个乘法拆成两轮,先来看第一轮:
uint64_t t = 0;
for(j = 0; j < 8; j++) {
t = k[0] * b[j] + t;
r[j] = t & 0x0FFFFFFFFULL;
t >>= 32;
}
r[8] = t;
2
3
4
5
6
7
大家可以把循环展开看一下,这里用t存过程计算结果,r每个元素只存乘法的低32位,高32位加到下一轮的结果中,循环结束后,最后的高32位加到r[8]中。
下面看第二轮:
uint64_t t = 0;
for(j = 0; j < 8; j++) {
t = r[1 + j] + k[1] * b[j] + t;
r[1 + j] = t & 0x0FFFFFFFFULL;
t >>= 32;
}
r[1 + 8] = t;
2
3
4
5
6
7
第二轮,j=0时,乘法结果应该存在r[1]中吧,另外还得加上第一轮r[1]中的数,这里的t会溢出吗?
结果是不会溢出的,因为r[1+j]是一个32位的数,t也是一个32位的数,乘法的结果加上两32位的数,最大值是每个字节刚好都为FF,大家可以用计算器算一下。
合并两轮循环,最终算法为:
static void gm_i_bn_mul(uint64_t * r, const uint64_t * k, const gm_bn_t b) {
int i, j;
// k0b0, k0b1+k1b0, .. ,k0b7+k1b6, k1b7
uint64_t t;
for(i = 0; i < 2; i++) {
t = 0;
for (j = 0; j < 8; j++) {
t = r[i + j] + k[i] * b[j] + t;
r[i + j] = t & 0x0FFFFFFFFULL;
t >>= 32;
}
r[i + 8] = t;
}
}
2
3
4
5
6
7
8
9
10
11
12
13
14
15
# t的计算
在计算q值过程中,还有一个t的计算,先回顾q的公式:
注意这里,我们对
伪代码为:
t = 1;
for(i = 1; i < 64; i++) {
t = t * t * p mod 2^64;
}
t = 2^64 - t;
2
3
4
5
代码实现:
uint64_t t = 1;
for(i = 1; i < 64; i++) {
t = t * t;
t = t * (m[0] | m[1] << 32);
}
t = 0xFFFFFFFFFFFFFFFF - t + 1;
2
3
4
5
6
这里的m即国密推荐曲线中的p或n。
最终计算结果分别是:
# 蒙哥马利模乘
先回顾一下伪代码:
d = 0;
for(i = 0; i < 4; i++) {
d = ((b[i] * a) + d + q * p) * 2^-64;
}
if(d > p) {
d = d - p;
}
2
3
4
5
6
7
再看看代码实现:
// 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) {
uint64_t temp_r[11] = {0};
uint64_t temp_mul[11] = {0};
uint64_t temp_q[2] = {0};
int i, j;
// t的值
uint64_t mod_t = (m[0] == GM_BN_N[0] ? 0x327F9E8872350975ULL : 0x01);
for(i = 0; i < 4; i++) {
// bi*a
memset(temp_mul, 0, sizeof(uint64_t) * 11);
gm_i_bn_mul(temp_mul, &b[i * 2], a);
gm_i_bn_add_x(temp_r, temp_r, temp_mul, 11);
// cal q
uint64_t q = (temp_r[0] | temp_r[1] << 32) * mod_t;
temp_q[0] = q & 0x0FFFFFFFFULL;
temp_q[1] = q >> 32;
// q * m
memset(temp_mul, 0, sizeof(uint64_t) * 11);
gm_i_bn_mul(temp_mul, temp_q, m);
gm_i_bn_add_x(temp_r, temp_r, temp_mul, 11);
// reduce temp_r
memmove(temp_r, temp_r + 2, sizeof(uint64_t) * 9);
temp_r[9] = 0;
temp_r[10] = 0;
}
// temp_r[8]为进位,将进位存放在temp_r[7]的高32位
temp_r[7] += temp_r[8] << 32;
if(gm_bn_cmp(temp_r, m) >= 0) {
gm_i_bn_sub(temp_r, temp_r, m);
}
gm_bn_copy(r, temp_r);
}
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
拆解来看,先是计算b[i] * a
,即:
// 结果存放在temp_mul
gm_i_bn_mul(temp_mul, &b[i * 2], a);
2
然后加上上一轮的结果
// 结果存放在temp_r
gm_i_bn_add_x(temp_r, temp_r, temp_mul, 11);
2
再计算q值,回顾公式:
然后是计算q * m,再将结果加到temp_r,循环体内最后是除以
函数最后判断temp_r是否大于等于m,条件成立就再减掉一个m即可。
这里的gm_i_bn_add_x是将gm_i_bn_add变化一下,因为我们的大数是8个uint64_t组成,但是这里每轮的乘法结果需要10个uint64_t才能存下,所以原来的加法就不满足了,另外这里用11个uint64_t来存放加法结果,是将最后可能的进位存放在下标为10的元素中。
static void gm_i_bn_add_x(gm_bn_t r, const gm_bn_t a, const gm_bn_t b, int count) {
int i;
r[0] = a[0] + b[0];
for(i = 1; i < count; i++) {
r[i] = a[i] + b[i] + (r[i - 1] >> 32);
r[i - 1] &= 0x0FFFFFFFFULL;
}
}
static void gm_i_bn_add(gm_bn_t r, const gm_bn_t a, const gm_bn_t b) {
gm_i_bn_add_x(r, a, b, 8);
}
2
3
4
5
6
7
8
9
10
11
12
# 单元测试
补充我们的单元测试案例,加上模乘
main函数增加:
TEST_BN_ALG("gmp_mul",
"32C4AE2C1F1981195F9904466A39C9948FE30BBFF2660BE1715A4589334C74C7",
"9DDD52AF95B748A553D1B1E106627F901CD453F067A0D50202C672130C90F607",
"64DD9339D3DFA3D15B581B1DD13E3D9202982F62473372E76B5D591A38F193CD");
TEST_BN_ALG("gmn_mul",
"32C4AE2C1F1981195F9904466A39C9948FE30BBFF2660BE1715A4589334C74C7",
"9DDD52AF95B748A553D1B1E106627F901CD453F067A0D50202C672130C90F607",
"6FE789E58A88991A4600A167FAF7F4F49058EF92ED5DDD701F1971356A674484");
2
3
4
5
6
7
8
9
test_bn函数增加:
else if(strcmp(alg, "mul") == 0){ // 1千万
gm_bn_to_mont(bnr, bna, m);
gm_bn_to_mont(bnb, bnb, m);
for (i = 0; i < 10000000; i++) {
gm_bn_mont_mul(bnr, bnr, bnb, m);
}
gm_bn_from_mont(bnr, bnr, m);
}
2
3
4
5
6
7
8
注意,这里进行模乘前,需要将参数转化为蒙哥马利域,最后结果需要转化回来
算法效率:
saintdeMacBook-Pro:bn saint$ time ./a.out gmp_mul
r = 64DD9339D3DFA3D15B581B1DD13E3D9202982F62473372E76B5D591A38F193CD
test result: ok
real 0m12.381s
user 0m11.873s
sys 0m0.123s
saintdeMacBook-Pro:bn saint$ time ./a.out gmn_mul
r = 6FE789E58A88991A4600A167FAF7F4F49058EF92ED5DDD701F1971356A674484
test result: ok
real 0m12.324s
user 0m11.944s
sys 0m0.124s
2
3
4
5
6
7
8
9
10
11
12
13
14