大数减法

接加法章节，已经约定减法中，a大于b，定义如下：

// r = a - b
static void gm_i_bn_sub(gm_bn_t r, const gm_bn_t a, const gm_bn_t b);

先看第一轮：

r[0] = a[0] - b[0];

虽然大数a大于b，但不能保证a[0]是大于b[0]的，所以我们可以先向高位借1，这样就没问题，那么从高位借1，应该是多大的数呢？

应该是0xFFFFFFFF + 1吧，所以调整一下第一轮：

r[0] = 0x100000000ULL + a[0] - b[0];

同样的，后续每一轮都向高位借1，第二轮应该是：

r[1] = 0x100000000ULL + a[0] - b[0] - 1;

这里为什么要减1呢，因为第一轮借了1，所以要减掉1。有遗漏吗？

上面提到a[0]不一定大于b[0]，所以提前向高位借了1，但如果a[0] >= b[0]呢，那向高位借的1就用不上了，就需要加回去：

r[1] = 0x100000000ULL + a[0] - b[0] - 1 + (r[0] >> 32);

简化一下：

r[1] = 0xFFFFFFFF + a[0] - b[0] + (r[0] >> 32); // 会溢出吗？
// 前一轮清除进位
r[0] &= 0x0FFFFFFFFULL;

再看最后一轮，最后一轮没有更高位了，没地方借1了，而且大数a >= b，所以也不需要借位

r[7] = a[7] - b[7] - 1 + (r[6] >> 32);

因为前一轮借了1位，所以要减1，并且要加上前一轮可能产生的进位

转换为循环，最终实现为：

static void gm_i_bn_sub(gm_bn_t r, const gm_bn_t a, const gm_bn_t b) {
    int i;
    r[0] = 0x100000000ULL + a[0] - b[0];
    for (i = 1; i < 7; i++) {
        r[i] = 0x0FFFFFFFFULL + a[i] - b[i] + (r[i - 1] >> 32);
        r[i - 1] &= 0x0FFFFFFFFULL;
    }
    r[i] = a[i] - b[i] + (r[i - 1] >> 32) - 1;
    r[i - 1] &= 0x0FFFFFFFFULL;
}

合并循环后：

static void gm_i_bn_sub(gm_bn_t r, const gm_bn_t a, const gm_bn_t b) {
    int i;
    uint64_t t = 1;
    for (i = 0; i < 8; i++) {
        t = 0x0FFFFFFFFULL + a[i] - b[i] + t;
        r[i] = t & 0x0FFFFFFFFULL;
        t = t >> 32;
    }
    r[i] += (t - 1) << 32;
}

实测合并循环的版本效率不如未合并的版本，执行2亿次，相差1秒左右。

大数模减

实现了大数减法后，再来实现模减，定义及实现如下：

// r = (a - b) mod m
void gm_bn_sub(gm_bn_t r, const gm_bn_t a, const gm_bn_t b, const gm_bn_t m) {
    if (gm_bn_cmp(a, b) >= 0) {
        gm_i_bn_sub(r, a, b);
    } else {
        gm_bn_t t;
        gm_i_bn_sub(t, m, b);
        gm_i_bn_add(r, t, a);
    }
}

因为我们实现的大数减法a必须大于等于b，所以要先进行判断

如果a < b，则可以先用m减去b，再与a相加，推导：

r = a - b + m = m - b + a;

因为a < b，所以a - b之后是一个负数，对m取模的话，加上m即可，m加一个负数，结果不可能超过m，所以不用考虑溢出进位问题。

单元测试

我们已经实现了大数模加及模减，为了验证算法的正确性及相应效率，需要写一个单元测试用的代码来进行测试，这里不详细介绍了，我先贴代码：

# include "gm.h"
#include <stdio.h>

#define TEST_BN_ALG(alg_name, a, b, r) \
    do { \
        if(strncmp(argv[1], alg_name, strlen(alg_name)) == 0) { \
            test_bn(a, b, r, alg_name); \
        } \
    }while(0)

/**
 * 大数单元测试
 * @param a 大数a
 * @param b 大数b
 * @param res 预期结果
 * @param ralg 测试的指令add、sub、mul等
 */
void test_bn(const char * a, const char * b, const char * res, const char * ralg) {
    int i, j;
    gm_bn_t bna, bnb, bnr;
    char bnr_hex[65] = {0};
    const uint64_t *m;
    const char *alg = ralg + 4;

    if (gm_bn_from_hex(bna, a) < 0) {
        printf("convert a to bn failed.\n");
    }

    if (gm_bn_from_hex(bnb, b) < 0) {
        printf("convert b to bn failed.\n");
    }

    if (strncmp(ralg, "gmp", 3) == 0) {
        m = GM_BN_P;
    } else {
        m = GM_BN_N;
    }

    if(strcmp(alg, "add") == 0) { // 2亿
        gm_bn_copy(bnr, bna);
        for (i = 0; i < 200000000; i++) {
            gm_bn_add(bnr, bnr, bnb, m);
        }
    }else if(strcmp(alg, "sub") == 0) { // 2亿
        gm_bn_copy(bnr, bna);
        for (i = 0; i < 200000000; i++) {
            gm_bn_sub(bnr, bnr, bnb, m);
        }
    }

    gm_bn_to_hex(bnr, bnr_hex);

    printf("r = %s\n", bnr_hex);
    printf("test result: %s\n", (strcmp(res, bnr_hex) == 0 ? "ok" : "fail"));
}

int main(int argc, char ** argv) {
    size_t cmdLen = strlen(argv[1]);
    if (cmdLen < 2) return -1;

    /**  bn mod p ops **/
    TEST_BN_ALG("gmp_add",
                "32C4AE2C1F1981195F9904466A39C9948FE30BBFF2660BE1715A4589334C74C7",
                "9DDD52AF95B748A553D1B1E106627F901CD453F067A0D50202C672130C90F607",
                "6564F159CE84BEC13D3CA303596C95F1E0FC286C6B6EC0B35CD378B33A84A721");

    TEST_BN_ALG("gmp_sub",
                "32C4AE2C1F1981195F9904466A39C9948FE30BBFF2660BE1715A4589334C74C7",
                "9DDD52AF95B748A553D1B1E106627F901CD453F067A0D50202C672130C90F607",
                "00246AFE6FAE437181F565897B06FD373EC9EF13795D570F85E1125F2C14426D");

    /**  bn mod n ops **/
    TEST_BN_ALG("gmn_add",
                "32C4AE2C1F1981195F9904466A39C9948FE30BBFF2660BE1715A4589334C74C7",
                "9DDD52AF95B748A553D1B1E106627F901CD453F067A0D50202C672130C90F607",
                "6564F159CE84BEC13D3CA3035D805623B71F26583B35FD3A6F9D86075603B079");

    TEST_BN_ALG("gmn_sub",
                "32C4AE2C1F1981195F9904466A39C9948FE30BBFF2660BE1715A4589334C74C7",
                "9DDD52AF95B748A553D1B1E106627F901CD453F067A0D50202C672130C90F607",
                "00246AFE6FAE437181F5658976F33D0568A6F127A9961A887317050B10953915");
}

约定

gmp_开头的指令，代表大数对国密参数p取模
gmn_开头的指令，代表大数对国密参数n取模

如果大家写的算法，跟我这里计算预期的结果不一致，就要去看看算法哪里有问题了。

其中gm_bn_copy其实就是内存复制：

#define gm_bn_copy(r, a) memcpy((r), (a), sizeof(gm_bn_t))

下面来看一下各算法对应的效率：

saintdeMacBook-Pro:bn saint$ time ./a.out gmp_add
r = 6564F159CE84BEC13D3CA303596C95F1E0FC286C6B6EC0B35CD378B33A84A721
test result: ok

real	0m12.510s
user	0m12.119s
sys	0m0.109s
saintdeMacBook-Pro:bn saint$ time ./a.out gmp_sub
r = 00246AFE6FAE437181F565897B06FD373EC9EF13795D570F85E1125F2C14426D
test result: ok

real	0m12.369s
user	0m12.175s
sys	0m0.079s
saintdeMacBook-Pro:bn saint$ time ./a.out gmn_add
r = 6564F159CE84BEC13D3CA3035D805623B71F26583B35FD3A6F9D86075603B079
test result: ok

real	0m12.279s
user	0m12.070s
sys	0m0.080s
saintdeMacBook-Pro:bn saint$ 
saintdeMacBook-Pro:bn saint$ time ./a.out gmn_sub
r = 00246AFE6FAE437181F5658976F33D0568A6F127A9961A887317050B10953915
test result: ok

real	0m12.051s
user	0m11.804s
sys	0m0.087s

我的电脑是15年买的MacBook Pro，那种厚厚的老款Pro，cpuinfo如下：

saintdeMacBook-Pro:bn saint$ sysctl machdep.cpu.brand_string
machdep.cpu.brand_string: Intel(R) Core(TM) i5-3210M CPU @ 2.50GHz

Linux环境未测试过，理论是没问题的，Windows也未测试过，可能不需要改动，也可能需要小改，这个应该不难。

后续单元测试

后续单元测试，就只写关键部分，不贴完整代码了，太长，跨大章节会把对应的代码贴出来。

途中感想

有没有感觉，其实大数相关操作也没想像中复杂，想想openssl的bn目录，里面那么多内容，感觉很复杂，还依赖了一堆的内存操作OPENSSL_malloc、OPENSSL_free，还有asn1相关内容，头大。