大数乘法

上一章节已简单的介绍了蒙哥马利相关算法，接上一章可知，我们将256-bit的乘法拆成了四轮，即每轮是$2^{64}$乘以一个256-bit的大数。

接下来我们先定义一个乘法函数：

// 乘法结果 r = k * b
static void gm_i_bn_mul(uint64_t * r, const uint64_t * k, const gm_bn_t b);

其中k是一个$2^{64}$的数，因为我们是用数组来表示，数组每个元素只存$2^{32}$的数，所以这里的k是一个含有2个元素的数组，乘法的结果理论uint64_t r[10]是可以存的下的，下面我们看一下乘法算法：

   b7b6b5b4b3b2b1b0
x              k1k0

结果（数组下标从0开始）：k0b0, k0b1+k1b0, .. ,k0b7+k1b6, k1b7

将这个乘法拆成两轮，先来看第一轮：

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;

大家可以把循环展开看一下，这里用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;

第二轮，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;
    }
}

t的计算

在计算q值过程中，还有一个t的计算，先回顾q的公式：

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

注意这里，我们对$2^{64}$取模，因为我们是将256-bit拆成四轮乘法，所以每轮是除以$2^{64}$

伪代码为：

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

代码实现：

uint64_t t = 1;
for(i = 1; i < 64; i++) {
    t = t * t;
    t = t * (m[0] | m[1] << 32);
}
t = 0xFFFFFFFFFFFFFFFF - t + 1;

这里的m即国密推荐曲线中的p或n。

最终计算结果分别是：

$t = -p^{-1}\pmod {2^{64}} = 0x01$

$t = -n^{-1}\pmod {2^{64}} = 0x327F9E8872350975ULL$

蒙哥马利模乘

先回顾一下伪代码：

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

再看看代码实现：

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

拆解来看，先是计算b[i] * a，即：

// 结果存放在temp_mul
gm_i_bn_mul(temp_mul, &b[i * 2], a);

然后加上上一轮的结果

// 结果存放在temp_r
gm_i_bn_add_x(temp_r, temp_r, temp_mul, 11);

再计算q值，回顾公式：

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

$(a[0]*b[i] + d_0)$其实就是(temp_r[0] | temp_r[1] << 32)。

然后是计算q * m，再将结果加到temp_r，循环体内最后是除以$2^{64}$，因为加上q*m后，结果是$2^{64}$倍数，所以直接把低64位抛弃掉即可。

函数最后判断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);
}

单元测试

补充我们的单元测试案例，加上模乘

main函数增加：

TEST_BN_ALG("gmp_mul",
                "32C4AE2C1F1981195F9904466A39C9948FE30BBFF2660BE1715A4589334C74C7",
                "9DDD52AF95B748A553D1B1E106627F901CD453F067A0D50202C672130C90F607",
                "64DD9339D3DFA3D15B581B1DD13E3D9202982F62473372E76B5D591A38F193CD");

TEST_BN_ALG("gmn_mul",
                "32C4AE2C1F1981195F9904466A39C9948FE30BBFF2660BE1715A4589334C74C7",
                "9DDD52AF95B748A553D1B1E106627F901CD453F067A0D50202C672130C90F607",
                "6FE789E58A88991A4600A167FAF7F4F49058EF92ED5DDD701F1971356A674484");

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

注意，这里进行模乘前，需要将参数转化为蒙哥马利域，最后结果需要转化回来

算法效率：

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

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

上一页蒙哥马利算法

下一页大数模平方