在Go，C#，Java，Erlang，Python等语言中都内置有大数的计算，笔者曾经写过一篇博文专门对比测试它们的性能：使用斐波那契（Fibonacci）数列来测试各大语言的性能
在C/C++中虽然标准库并没有内置大数的计算，但是有可选择的高性能的第三方库，比较有名的是GMP库和OpenSSL的大数计算。其中GMP库号称是全球最快的大数计算库，没有之一。
本文就简单介绍一下如何使用它，然后再使用它写一个斐波那契（Fibonacci）数列来看看性能，可以对比一下前述的几大语言的测试。
一、安装
在MinGW下使用命令安装：
pacman -S mingw-w64-x86_64-gmp
Ubuntu下使用命令安装：
sudo apt install libgmp-dev
MSVC使用，可以去官方下载源码自行编译。
二、常用函数简介
gmp库的类型分为整数mpz_t，有理数mpq_t，浮点数mpf_t，自然数mp_limb_t。整数类的函数以mpz_开头，有理数以mpq_开头，浮点数以mpf_开头，还有一些非常底层的操作自然数的函数以mpn_开头，它们快速，高效，但是难以使用。
由于GMP库是C语言API，所以需要遵循C语言的使用规范。所有的数据在使用前都需要初始化（分配内存）再使用，使用完后还必须释放资源（释放内存）。
下面以整数为例介绍一些常用的函数，其它的可以类推。
1. 初始化
void mpz_init (mpz_t x)
初始化x为0。
void mpz_inits (mpz_t x, …)
如果有多个变量需要初始化，可以使用此函数一次性初始化完，要求以0结尾，比如mpz_inits (a, b, c, 0)。
2. 释放
void mpz_clear (mpz_t x)
释放x使用的内存
void mpz_clears (mpz_t x, …)
同初始化一样，如果有多个变量需要释放内存，可以使用此函数一次性释放完，同样要求以0结尾，mpz_clears (a, b, c, 0)。
3. 赋值
void mpz_set (mpz_t rop, const mpz_t op)
void mpz_set_ui (mpz_t rop, unsigned long int op)
void mpz_set_si (mpz_t rop, signed long int op)
void mpz_set_d (mpz_t rop, double op)
void mpz_set_q (mpz_t rop, const mpq_t op)
void mpz_set_f (mpz_t rop, const mpf_t op)
这些函数都是将op赋值给rop。其中 mpz_set_d, mpz_set_q 和 mpz_set_f 会截断 op为一个整数.
int mpz_set_str (mpz_t rop, const char *str, int base)
使用字符串str中的数据来赋值给rop，base代表基数，可以是2到62，如果是0，则使用字符串str中的数据前导来决定，0x/0X表示十六进制，0b/0B表示二进制，0表示八进制，其它表示十进制。对于基数不超过36的情况，字母大小写不做区分；大写字母与小写字母的取值相同。对于基数在37至62之间的情况，大写字母对应常规的10…35，小写字母对应36…61。
若整个字符串是指定进制下的有效数字，该函数返回0；否则返回-1。
void mpz_swap (mpz_t rop1, mpz_t rop2)
交换两个数的值。
4. 初始与赋值的组合
void mpz_init_set (mpz_t rop, const mpz_t op)
void mpz_init_set_ui (mpz_t rop, unsigned long int op)
void mpz_init_set_si (mpz_t rop, signed long int op)
void mpz_init_set_d (mpz_t rop, double op)
int mpz_init_set_str (mpz_t rop, const char *str, int base)
在定义了mpz_t整数后，如果想初始化与赋值一起，可以使用上面这些函数。但是需要注意的是不能对已经初始化了的mpz_t整数调用上面这些函数，这样会导致内存泄漏。
5. 转换
unsigned long int mpz_get_ui (const mpz_t op)
signed long int mpz_get_si (const mpz_t op)
double mpz_get_d (const mpz_t op)
char * mpz_get_str (char *str, int base, const mpz_t op)
前面的赋值是将C语言中的标准类型转换为GMP中的mpz_t类型，如果想要转换回来就使用上面这些函数。
double mpz_get_d_2exp (signed long int *exp, const mpz_t op)
将操作数转换为双精度浮点数，必要时进行截断（即向零舍入），并单独返回指数。
返回值的范围为 0.5<=abs(d)<1，指数会存储至 *exp。d * 2^exp 即为经截断处理的操作数数值。若操作数为零，则返回值为 0.0，且 *exp 中存储的值为 0。
这与标准 C 语言的 frexp 函数类似。
6. 算术运算
void mpz_add (mpz_t rop, const mpz_t op1, const mpz_t op2)
void mpz_add_ui (mpz_t rop, const mpz_t op1, unsigned long int op2)
加法运算：
r
o
p
=
o
p
1
+
o
p
2
rop = op1 + op2
rop=op1+op2
void mpz_sub (mpz_t rop, const mpz_t op1, const mpz_t op2)
void mpz_sub_ui (mpz_t rop, const mpz_t op1, unsigned long int op2)
void mpz_ui_sub (mpz_t rop, unsigned long int op1, const mpz_t op2)
减法运算：
r
o
p
=
o
p
1
−
o
p
2
rop = op1 - op2
rop=op1−op2
Function: void mpz_mul (mpz_t rop, const mpz_t op1, const mpz_t op2)
Function: void mpz_mul_si (mpz_t rop, const mpz_t op1, long int op2)
Function: void mpz_mul_ui (mpz_t rop, const mpz_t op1, unsigned long int op2)
乘法运算：
r
o
p
=
o
p
1
∗
o
p
2
rop = op1 * op2
rop=op1∗op2
Function: void mpz_addmul (mpz_t rop, const mpz_t op1, const mpz_t op2)
Function: void mpz_addmul_ui (mpz_t rop, const mpz_t op1, unsigned long int op2)
加乘积运算：
r
o
p
+
=
o
p
1
∗
o
p
2
rop += op1 * op2
rop+=op1∗op2
Function: void mpz_submul (mpz_t rop, const mpz_t op1, const mpz_t op2)
Function: void mpz_submul_ui (mpz_t rop, const mpz_t op1, unsigned long int op2)
减乘积运算：
r
o
p
−
=
o
p
1
∗
o
p
2
rop -= op1 * op2
rop−=op1∗op2
Function: void mpz_mul_2exp (mpz_t rop, const mpz_t op1, mp_bitcnt_t op2)
运算：
r
o
p
=
o
p
1
∗
2
o
p
2
rop = op1 * 2^{op2}
rop=op1∗2op2，该操作也可定义为按op2位进行左移。
Function: void mpz_neg (mpz_t rop, const mpz_t op)
负数运算：
r
o
p
=
−
o
p
rop = -op
rop=−op
Function: void mpz_abs (mpz_t rop, const mpz_t op)
取绝对值运算
void mpz_cdiv_q (mpz_t q, const mpz_t n, const mpz_t d)
void mpz_cdiv_r (mpz_t r, const mpz_t n, const mpz_t d)
void mpz_cdiv_qr (mpz_t q, mpz_t r, const mpz_t n, const mpz_t d)
unsigned long int mpz_cdiv_q_ui (mpz_t q, const mpz_t n, unsigned long int d)
unsigned long int mpz_cdiv_r_ui (mpz_t r, const mpz_t n, unsigned long int d)
Function: unsigned long int mpz_cdiv_qr_ui (mpz_t q, mpz_t r, const mpz_t n, unsigned long int d)
unsigned long int mpz_cdiv_ui (const mpz_t n, unsigned long int d)
void mpz_cdiv_q_2exp (mpz_t q, const mpz_t n, mp_bitcnt_t b)
void mpz_cdiv_r_2exp (mpz_t r, const mpz_t n, mp_bitcnt_t b)
void mpz_fdiv_q (mpz_t q, const mpz_t n, const mpz_t d)
void mpz_fdiv_r (mpz_t r, const mpz_t n, const mpz_t d)
void mpz_fdiv_qr (mpz_t q, mpz_t r, const mpz_t n, const mpz_t d)
unsigned long int mpz_fdiv_q_ui (mpz_t q, const mpz_t n, unsigned long int d)
unsigned long int mpz_fdiv_r_ui (mpz_t r, const mpz_t n, unsigned long int d)
unsigned long int mpz_fdiv_qr_ui (mpz_t q, mpz_t r, const mpz_t n, unsigned long int d)
unsigned long int mpz_fdiv_ui (const mpz_t n, unsigned long int d)
void mpz_fdiv_q_2exp (mpz_t q, const mpz_t n, mp_bitcnt_t b)
void mpz_fdiv_r_2exp (mpz_t r, const mpz_t n, mp_bitcnt_t b)
void mpz_tdiv_q (mpz_t q, const mpz_t n, const mpz_t d)
void mpz_tdiv_r (mpz_t r, const mpz_t n, const mpz_t d)
void mpz_tdiv_qr (mpz_t q, mpz_t r, const mpz_t n, const mpz_t d)
unsigned long int mpz_tdiv_q_ui (mpz_t q, const mpz_t n, unsigned long int d)
unsigned long int mpz_tdiv_r_ui (mpz_t r, const mpz_t n, unsigned long int d)
unsigned long int mpz_tdiv_qr_ui (mpz_t q, mpz_t r, const mpz_t n, unsigned long int d)
unsigned long int mpz_tdiv_ui (const mpz_t n, unsigned long int d)
void mpz_tdiv_q_2exp (mpz_t q, const mpz_t n, mp_bitcnt_t b)
void mpz_tdiv_r_2exp (mpz_t r, const mpz_t n, mp_bitcnt_t b)
除法运算
将 n 除以 d，得到商 q 和余数 r。对于 2exp 函数，
d
=
2
b
d=2^b
d=2b。舍入共有三种模式，分别适用于不同场景：
- cdiv 函数会将 q 向正无穷方向向上取整，且 r 的符号与 d 相反。其中字母 c 是“ceil（向上取整）”的缩写。
- fdiv 会将商 q 向负无穷方向向下取整，余数 r 的符号将与除数 d 一致。其中的 f 代表“floor（向下取整）”。
- tdiv 会将 q 向零方向取整，且 r 的符号将与 n 一致。其中 t 代表“截断（truncate）”。
在所有情况下，q 和 r 都满足 n = q ∗ d + r n = q*d + r n=q∗d+r，且 r 满足 0 ≤ ∣ r ∣ < ∣ d ∣ 0 ≤ |r| < |d| 0≤∣r∣<∣d∣。
7.取模运算
void mpz_mod (mpz_t r, const mpz_t n, const mpz_t d)
unsigned long int mpz_mod_ui (mpz_t r, const mpz_t n, unsigned long int d)
将 r 设为 n 对 d 取模的结果。该运算会忽略除数的符号，所得结果始终为非负值。
8. 指数运算
void mpz_pow_ui (mpz_t rop, const mpz_t base, unsigned long int exp)
void mpz_ui_pow_ui (mpz_t rop, unsigned long int base, unsigned long int exp)
将rop设置为base的exp次幂，即：
r
o
p
=
b
a
s
e
e
x
p
rop=base^{exp}
rop=baseexp。其中0的0次方的结果为1。
void mpz_powm (mpz_t rop, const mpz_t base, const mpz_t exp, const mpz_t mod)
void mpz_powm_ui (mpz_t rop, const mpz_t base, unsigned long int exp, const mpz_t mod)
将rop设置为（base的exp次方）对mod取模的结果。
9. 求根运算
void mpz_sqrt (mpz_t rop, const mpz_t op)
将 rop 设置为操作数平方根的截断整数部分，即：
r
o
p
=
o
p
rop=\sqrt{op}
rop=op
int mpz_root (mpz_t rop, const mpz_t op, unsigned long int n)
void mpz_rootrem (mpz_t root, mpz_t rem, const mpz_t u, unsigned long int n)
void mpz_sqrtrem (mpz_t rop1, mpz_t rop2, const mpz_t op)
int mpz_perfect_power_p (const mpz_t op)
int mpz_perfect_square_p (const mpz_t op)
10. 比较运算
int mpz_cmp (const mpz_t op1, const mpz_t op2)
int mpz_cmp_d (const mpz_t op1, double op2)
int mpz_cmp_si (const mpz_t op1, signed long int op2)
int mpz_cmp_ui (const mpz_t op1, unsigned long int op2)
比较 op1 与 op2。若 op1 > op2 则返回正值，若 op1 = op2 则返回零，若 op1 < op2 则返回负值。
mpz_cmp_ui 和 mpz_cmp_si 均为宏，会对其参数进行多次求值。调用 mpz_cmp_d 时可传入无穷大值，但传入 NaN 时结果未定义。
int mpz_cmpabs (const mpz_t op1, const mpz_t op2)
int mpz_cmpabs_d (const mpz_t op1, double op2)
int mpz_cmpabs_ui (const mpz_t op1, unsigned long int op2)
比较操作数op1与op2的绝对值：若abs(op1) > abs(op2)则返回正值，若abs(op1) = abs(op2)则返回零，若abs(op1) < abs(op2)则返回负值。
mpz_cmpabs_d 可传入无穷大参数调用，但传入 NaN 时结果未定义。
int mpz_sgn (const mpz_t op)
若 op > 0 则返回 +1，若 op = 0 则返回 0，若 op < 0 则返回 -1。
该函数实际是以宏的形式实现的，会对其参数进行多次求值。
11. 逻辑与位运算
void mpz_and (mpz_t rop, const mpz_t op1, const mpz_t op2)
位与运算：rop = op1 & op2
void mpz_ior (mpz_t rop, const mpz_t op1, const mpz_t op2)
位或运算：rop = op1 | op2
void mpz_xor (mpz_t rop, const mpz_t op1, const mpz_t op2)
位异或运算：rop = op1 ^ op2
void mpz_setbit (mpz_t rop, mp_bitcnt_t bit_index)
置位运算，将rop的第bit_index位设置为1
void mpz_clrbit (mpz_t rop, mp_bitcnt_t bit_index)
清位运算，将rop的第bit_index位设置为0
int mpz_tstbit (const mpz_t op, mp_bitcnt_t bit_index)
测试rop的第bit_index位，根据值返回1或者0
12. 格式化输出
A、格式化输出字符
GMP库的格式化输出函数与标准C的格式化输出非常相似。只是添加了自己的格式化输出字符：
GMP 分别为 mpz_t、mpq_t 和 mpf_t 新增了类型“Z”、“Q”和“F”，为 mp_limb_t 新增了类型“M”，为 mp_limb_t 数组新增了类型“N”。“Z”、“Q”、“M”和“N”的表现与整数类似。若有需要，“Q”会打印出“/”及分母。“F”的表现与浮点数类似。
B、格式化输出函数
int gmp_printf (const char *fmt, …)
int gmp_vprintf (const char *fmt, va_list ap)
输出至标准输出流 stdout。返回写入的字符数，若发生错误则返回 -1。
int gmp_fprintf (FILE *fp, const char *fmt, …)
int gmp_vfprintf (FILE *fp, const char *fmt, va_list ap)
输出至流 fp。返回写入的字符数，若出错则返回 -1。
int gmp_sprintf (char *buf, const char *fmt, …)
int gmp_vsprintf (char *buf, const char *fmt, va_list ap)
在缓冲区（buf）中构造一个以空字符结尾的字符串，返回写入的字符数，不包含结尾的空字符。
buf 指向的存储空间与 fmt 字符串之间不允许存在重叠区域。
不建议使用这些函数，因为没有任何防护机制可防止缓冲区 buf 的可用空间被超出。
int gmp_snprintf (char *buf, size_t size, const char *fmt, …)
int gmp_vsnprintf (char *buf, size_t size, const char *fmt, va_list ap)
在 buf 中构造一个以空字符结尾的字符串。写入的字节数不会超过 size。若要获得完整输出，size 的值必须足够容纳该字符串及其结尾的空字符。
返回值为应生成的字符总个数，不包含结尾的空字符。若返回值 retval 大于等于 size，则实际输出会被截断为前 size-1 个字符，并在末尾追加一个空字符。
区域 {buf,size} 与 fmt 字符串之间不允许存在重叠。
请注意，返回值采用的是 ISO C99 snprintf 风格，即便 C 库中的 vsnprintf 是较早的 GLIBC 2.0.x 风格，情况也是如此。
三、计算斐波那契（Fibonacci）数列
C源码，main.c：
#include <gmp.h>
int main(int argc, char* argv[]) {
mpz_t a;
mpz_t b;
mpz_t sum;
mpz_init(sum);
mpz_init_set_si(a, 1);
mpz_init_set_si(b, 1);
for (int i = 0; i < 1000000; ++i)
{
mpz_add(sum, a, b);
mpz_set(b, a);
mpz_set(a, sum);
}
gmp_printf("%Zd\n", sum);
mpz_clears(a, b, sum, 0);
return 0;
}
CMakeLists.txt：
cmake_minimum_required(VERSION 3.20)
project(t)
find_package(PkgConfig)
pkg_check_modules(pkgLibs REQUIRED IMPORTED_TARGET gmp)
aux_source_directory(. SRC)
add_executable(${PROJECT_NAME} ${SRC})
target_include_directories(${PROJECT_NAME} PRIVATE ${pkgLibs_INCLUDE_DIRS})
target_link_libraries(${PROJECT_NAME} PRIVATE PkgConfig::pkgLibs)
或者直接使用命令编译：
gcc -O2 -o t main.c -lgmp
在MinGW下使用：time ./t运行结果(多次运行，以最少的一次为准)：
测试了一下Debug以及Release版本，发现基本没差别。结果有点出乎意料啊，不得不佩服Go的性能！
为了公平，这里的C源码并没有使用GMP计算斐波那契数的专用函数：mpz_fib_ui/mpz_fib2_ui，因为其它语言没有这样的专用函数，使用了就不公平。
虽说Go语言在循环中使用了Set：
b.Set(a)
a.Set(sum)
C语言也使用了set：
mpz_set(b, a);
mpz_set(a, sum);
但Go语言应该对这种情况的内存分配做了极致优化，而mpz_set是内存的深度复制，会严重影响性能，而它有另一个API：mpz_swap，可以有效避免深度复制，做到重复利用。
改一下测试代码：
#include <gmp.h>
int main(int argc, char* argv[]) {
mpz_t a;
mpz_t b;
mpz_t c;
mpz_init(c);
mpz_init_set_si(a, 1);
mpz_init_set_si(b, 1);
for (int i = 0; i < 1000000; ++i) {
mpz_add(c, a, b);
mpz_swap(b, a);
mpz_swap(a, c);
}
gmp_printf("%Zd\n", a); // 最后结果在a变量中
mpz_clears(a, b, c, 0);
return 0;
}
测试结果：
这回终于符合预期了，只用了1.932秒，比Go快三分之一左右。
那如果使用GMP库计算斐波那契数的专用函数进行测试呢：
#include <gmp.h>
int main(int argc, char* argv[]) {
mpz_t a;
mpz_t b;
mpz_init(a);
mpz_init(b);
mpz_fib2_ui(a, b, 1000002); // 由于原来的计算是预先设置了两项：1，1，所以这里需要加2
gmp_printf("%Zd\n", a);
mpz_clears(a, b, 0);
return 0;
}
结果也是出乎意料，这回是出乎意料的快：只需要0.095秒，果然GMP是大数计算的性能之王。
各测试对比表：
| 语言 | Code Runner结果 | Code Runner 执行的命令 | 命令行编译运行结果 | 命令行编译运行的命令 | 命令行直接运行结果 | 命令行直接运行的命令 |
|---|---|---|---|---|---|---|
| Go | 3.548秒 | go run main.go | 3.539秒 | time go run main.go | 3.068秒 | time ./main |
| C# | 10.295秒 | dotnet run | 10.075秒 | time dotnet run | 9.255秒 | time ./ConsoleApp |
| Java | 10.654秒 | javac t.java && java t | / | / | 9.814秒 | time java t |
| Python | 5.939秒 | set PYTHONIOENCODING=utf-8 && python.exe -u main.py | 5.943秒 | time python main.py | / | / |
| Erlang | 13.353秒 | escript app.erl | / | / | 12.836秒 | time erl -noshell -s app main 0 -s init stop |
| C语言迭代法使用mpz_set | / | / | / | / | 4.29秒 | time ./t |
| C语言迭代法使用mpz_swap | / | / | / | / | 1.932秒 | time ./t |
| C语言使用专用函数mpz_fib2_ui | / | / | / | / | 0.095秒 | time ./t |
总结：Go语言不愧是集快速开发与高性能平衡得最好的语言。
扫一扫
6594
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
