大貓咪的新玩具GMP Library: C/C++ 任意精度的整數和浮點數
目錄
在算排列的時候發現 64 位元的整數不夠用
於是換成了 128 位元的整數
但發現還是不夠用…
GMP
GMP 是一個 Library，他支援任意精度的運算，裡面包括整數、有理數和無理數的型態，最初在 1991 年的時候釋出，到現今已經 32 年了，最新的版本是 6.2.1，Library 撰寫是用 C 語言，因此介面提供的 C 語言，但到現在有 C++ 的介面可以使用，另外也有其他語言的包裝介面可以使用。
編譯 Compile
進入資料夾後，在 bash 上運行以下指令：
- 預設位置在
/usr/local，可以加上--prefix=<Path>更改安裝位置 - 若要支援 C++，加上
--enable-cxx
./configureNo usable m4 in $PATH
如果中途有發生找不到 m4 的話，可以下指令安裝 m4，並重新執行一次
sudo apt-get install m4生成好 Makefile 後就開始編譯：
make編譯完可以跑內建的測試檢查一下：
make check測試完沒問題可以安裝：
make installCross Compile
由於 Windows 上沒有一些所需的程式可以使用，因此在檢查環境的時候就會抱錯，在 Windows 上官方建議使用 cygwin，MSYS2 應該也是可以使用，然後在上面的 bash 安裝所需的程式再編譯即可。但因為我不想安裝一堆東西，所以我選擇了 Cross Compile，在 Ubuntu 裡面編譯 Windows 的 library。
這邊使用 MinGW 的 Cross Compiler，所以先安裝他：
sudo apt-get install gcc-mingw-w64再來設定參數，因為要在 Linux 中 build，所以設定 CC_FOR_BUILD 的數值為 x86_64-linux-gnu-gcc，build 的參數為 x86_64-linux-gnu，另外因為我們最終要在 Windows 上使用，因此 CC 設為 x86_64-w64-mingw32-gcc，host 的參數設為 x86_64-w64-mingw32。
export CC=x86_64-w64-mingw32-gcc
export CC_FOR_BUILD=x86_64-linux-gnu-gcc
./configure --build=x86_64-linux-gnu --host=x86_64-w64-mingw32最後一樣直接 make 開始編譯，然後檢查跟安裝
make
make check
make install使用
在使用的時候要先引入 header
- C:
gmp.h - C++:
gmpxx.h
#include "gmp.h" // C
#include "gmpxx.h" // C++編譯的時候最後要連結 library
- C: -lgmp
- C++: -lgmpxx
記得設定 library 的路徑和 include 的路徑到 GMP 安裝的地方
# C
gcc main.c -o main.exe -L <GMP Install Folder>/lib -I <GMP Install Folder>/include -lgmp
# C++
g++ main.cpp -o main.exe -L <GMP Install Folder>/lib -I <GMP Install Folder>/include -lgmpxx整數 mpz_t
整數數字的型態為 mpz_t，下面列出 C 的 function
如果想要用 C++ 介面的可以看 C++ Interface Integers
// 宣告整數變數
mpz_t num, num1, num2;
// 初始變數，為變數分配記憶體空間
mpz_init(num); // 初始一個
mpz_inits(num1, num2, NULL); // 初始很多個，最後一個一定要是 NULL
mpz_init_set_ui(num, 123456); // 使用 unsigned long int 初始
mpz_init_set_si(num, -123456); // 使用 signed long int 初始
mpz_init_set_d(num, 11111.0); // 使用 double 初始
mpz_init_set_str(num, "123456789123456789", 10); // 使用字串初始
// 設定數值
mpz_set(num, num1);
mpz_set_ui(num, 123456); // ui: unsigned long int
mpz_set_si(num, -123456); // si: signed long int
mpz_set_d(num, 11111.0); // d: double
// 轉換
char buf[1024];
mpz_get_ui(num); // ui: unsigned long int
mpz_get_si(num); // si: signed long int
mpz_get_d(num); // d: double
mpz_get_str(buf, 10, num); // 將 num 以 10 進位表示的方式轉成字串放到 buf
// 從檔案設定數值
mpz_inp_str(num, <FILE*>, 10)
mpz_inp_str(num, stdin, 10); // 從標準輸入
// 輸出數值到檔案
mpz_out_str(<FILE*>, 10, num);
mpz_out_str(stdout, 10, num); // 輸出到標準輸出
// 加法
// num = num1 + num2
mpz_add_ui(num, num1, 654321); // ui: unsigned long int
mpz_add(num, num1, num2);
// 減法
// num = num1 - num2
mpz_sub(num, num1, num2);
mpz_sub_ui(num, num1, 654321); // ui: unsigned long int
mpz_ui_sub(num, 654321, num1); // num = 654321 - num1
// 乘法
// num = num1 * num2
mpz_mul(num, num1, num2);
mpz_mul_si(num, num1, -654321); // si: signed long int
mpz_mul_ui(num, num1, 654321); // ui: unsigned long int
// 先乘後加
// num = num + num1 * num2
mpz_addmul(num, num1, num2);
mpz_addmul_ui(num, num1, 654321); // ui: unsigned long int
// 先乘後減
// num = num - num1 * num2
mpz_submul(num, num1, num2);
mpz_submul_ui(num, num1, 654321); // ui: unsigned long int
// 負數
// num = -num1
mpz_neg(num, num1);
// 絕對值
// num = abs(num1)
mpz_abs(num, num1);
// 交換
mpz_swap(num1, num2);
// 比較
// > 0: num1 > num2
// = 0: num1 = num2
// < 0: num1 < num2
mpz_cmp(num1, num2);
mpz_cmp_d(num, 22222.0); // d: double
mpz_cmp_si(num, -654321); // si: signed long int
mpz_cmp_ui(num, 654321); // ui: unsigned long int
// bitwise operation
// num = num1 <operation> num2
mpz_and(num, num1, num2); // &
mpz_ior(num, num1, num2); // | (inclusive or)
mpz_xor(num, num1, num2); // ^ (exclusive or)
mpz_com(num, num1); // ~ (one's complement)
mpz_popcount(num); // number of 1-bit
mpz_hamdist(num1, num2); // hamming distance
mpz_setbit(num, 0); // set bit
mpz_clrbit(num, 0); // clear bit
mpz_combit(num, 0); // complment bit
mpz_tstbit(num, 0); // test bit, 是 1 就回傳 1，是 0 回傳 0
// 釋放記憶體空間
// 不要使用後記得釋放空間
mpz_clear(num); // 釋放一個
mpz_clears(num1, num2, NULL); // 釋放多個，最後要加上 NULL
有理數 mpq_t
有理數是 mpq_t，所有操作都假設有理數都是 canonical form
canonical form:
- 最簡分數
- 分母是正數
- 數值是 0 的話分子是 0，分母是 1，也就是 (0/1)
// 宣告整數變數
mpq_t num, num1, num2;
mpq_init(num);
mpq_inits(num1, num2, NULL);
// 設定數值的時候，如果數值不是 canonical form 的話，要在設定完後，使用 mpq_canonicalize 做轉換再做其他運算
mpq_set(num, num1);
mpq_set_ui(num, 2, 654321); // 2/654321
mpq_set_si(num, 2, -654321);
mpq_set_d(num, 66666.0);
mpq_set_str(num, "41/8", 10); // base 10
mpq_swap(num1, num2);
// 轉換
char buf[1024];
mpq_get_d(num); // 轉成 double，必要時會截斷
mpq_get_str(buf, 10, num); // base 10
// 將 num 轉成 canonical form
mpq_canonicalize(num);
mpq_add(num, num1, num2); // num = num1 + num2
mpq_sub(num, num1, num2); // num = num1 - num2
mpq_mul(num, num1, num2); // num = num1 * num2
mpq_div(num, num1, num2); // num = num1 / num2
mpq_neg(num, num1); // num = -num1
mpq_abs(num, num1); // num = abs(num1)
mpq_inv(num, num1); // num = 1/num1;
// 比較
// > 0: num1 > num2
// = 0: num1 = num2
// < 0: num1 < num2
mpq_cmp(num1, num2); // 比較 num1 和 num2
mpq_cmp_ui(num, 2, 654321); // 比較 num 和 2/654321
mpq_cmp_si(num, -2, 654321); // 比較 num 和 -2/654321
// 比較 num1 和 num2 是否相等
// (要知道是否相等，用這個比 cmp 快)
// 0: 不相等
// 非0: 相等
mpq_equal(num1, num2);
// 取得正負號
// +1: num > 0
// 0: num = 0
// -1: num < 0
mpq_sgn(num);
// 取得分子和分母的 reference
// 類型是整數 mpz_t，可以直接使用 mpz_t 的函數
mpq_numref(num); // 分子
mpq_denref(num); // 分母
mpq_get_num(mpz_t numerator, num); // 取得 num 的分子，放到 numerator
mpq_get_den(mpz_t denominator, num); // 取得 num 的分母，放到 denominator
mpq_set_num(num, const mpz_t numerator); // 將 numerator 放到 num 的分子
mpq_set_den(num, const mpz_t denominator); // 將 denominator 放到 num 的分母
// 輸入 & 輸出
// base 10
mpq_inp_str(num, FILE*, 10);
mpq_inp_str(num, stdin, 10); // 標準輸入
mpq_out_str(FILE*, 10, num);
mpq_out_str(stdout, 10, num); // 標準輸出
mpq_clear(num);
mpq_clears(num1, num2, NULL);浮點數 mpf_t
在這個類別裡，使用者可以自行選擇有效位數 (mantissa) 的位數，可以在任何時候調整，來調整計算的精度。
指數 (exponent) 的部分則是固定的，在 32 位元的系統上，可以表示的範圍大概是 $2^{-68719476768}$ 到 $2^{68719476736}$。
GMP 的 mptf_t 並不是 IEEE 754 的擴展，因此在不同 word size 的機器上計算的結果會不一樣，所以建議使用 MPFR 代替，MPFR 就是 IEEE 754 的擴展。
mpf_t num, num1, num2;
mpf_set_default_prec(20); // 設定預設的精度為 20 bits
mpf_get_default_prec(); // 取得預設的精度
mpf_set_prec(num, 15); // 設定精度為 15 bits
mpf_get_prec(num);
mpf_init(num);
mpf_inits(num1, num2, NULL);
mpf_init2(num, 10); // 初始並設定精度為 10 bits
mpf_set(num1, num2);
mpf_set_ui(num, 123456);
mpf_set_si(num, -654321);
mpf_set_d(num, 123.321);
mpf_set_str(num, "123@-2", 10); // 123 的 -2 次方，10 進位
// 初始並設定
mpf_init_set(num, num1);
mpf_init_set_ui(num, 123456);
mpf_init_set_si(num, -654321);
mpf_init_set_d(num, 123.321);
mpf_init_set_str(num, "123@-2", 10);
mpf_add(num, num1, num2); // num = num1 + num2
mpf_add_ui(num, num1, 123456); // num = num1 + 123456
mpf_sub(num, num1, num2); // num = num1 - num2
mpf_ui_sub(num, 123456, num1); // num = 123456 - num1
mpf_sub_ui(num, num1, 123456); // num = num1 - 123456
mpf_mul(num, num1, num2); // num = num1 * num2
mpf_mul_ui(num, num1, 123456); // num = num1 * 123456
mpf_div(num, num1, num2); // num = num1 / num2
mpf_ui_div(num, 123456, num1); // num = 123456 / num1
mpf_div_ui(num, num1, 123456); // num = num1 / 123456
mpf_neg(num, num1); // num = -num1
mpf_abs(num, num1); // num = abs(num1)
// 比較
// > 0: num1 > num2
// = 0: num1 = num2
// < 0: num1 < num2
mpf_cmp(num1, num2);
mpf_cmp_d(num, 66666.6);
mpf_cmp_ui(num, 123456);
mpf_cmp_si(num, -654321);
mpf_sgn(num);
// 輸入輸出
mpf_inp_str(num, FILE*, 10); // 10 進位
mpf_inp_str(num, stdin, 10); // 10 進位
mpf_out_str(FILE*, 10, 20, num); // 10 進位，精度 20
mpf_out_str(stdout, 10, 20, num); // 10 進位，精度 20
mpf_ceil(num, num1);
mpf_floor(num, num1);
mpf_trunc(num, num1); // 直接捨去小數，所以當數字是正的時候，與 floor 一樣，當數字是負的時候，與 ceil 一樣
// 轉換
char buf[1024];
mp_exp_t exp;
mpf_get_d(num); // 轉成 double
mpf_get_si(num); // 去掉小數部分，轉成 signed long int
mpf_get_ui(num); // 去掉小數部分，轉成 unsigned long int
mpf_get_str(buf, &exp, 10, 500, num); // 10 進位，最多轉換 500 個 digit，設定為 0 則是轉換最多的 digit
mpf_swap(num1, num2);
mpf_clear(num);
mpf_clears(num1, num2, NULL);標準輸入 & 標準輸出
要正確使用 printf 和 scanf 的話要使用 gmp 版本的：
- gmp_scanf
- gmp_printf
標準輸入輸出符號
mpz_t: Z，%Zdmpq_t: Q，%Qdmpf_t: F，%Ff
輸入範例：
- 記得要先初始變數再接收輸入
- 不用加上
&，因為已經是 call by reference 了
mpz_t z; mpz_init(z);
gmp_printf("Input Z value: ");
gmp_scanf("%Zd", z);
gmp_printf("Z value: %Zd\n", z);
mpq_t q; mpq_init(q);
gmp_printf("Input Q value: ");
gmp_scanf("%Qd", q);
gmp_printf("Q value: %Qd\n", q);
mpf_t f; mpf_init(f);
gmp_printf("Input F value: ");
gmp_scanf("%Ff", f);
gmp_printf("F value: %.Ff\n", f);
########### Output ###########
Input Z value: 123123123123123123123123
Z value: 123123123123123123123123
Input Q value: 321321321321321/999999911111
Q value: 321321321321321/999999911111
Input F value: 1231231312312.32132132132132132
F value: 1231231312312.32132132輸出範例：
mpz_t z; mpz_init_set_d(z, 13579246810.0);
gmp_printf("Z value: %Zd\n", z);
mpq_t q; mpq_init(q); mpq_set_ui(q, 123456789, 321654987);
gmp_printf("Q value: %Qd\n", q);
mpf_t f; mpf_init_set_str(f, "1.111111111111111@-2", 10);
gmp_printf("F value: %.Ff\n", f);
########### Output ###########
Z value: 13579246810
Q value: 123456789/321654987
F value: 0.01111111111111111範例
#include <stdio.h>
#include <stdlib.h>
#include "gmp.h"
int main(){
mpz_t tmp;
mpz_init_set_ui(tmp, 1);
int limit = 100;
for(int i = 0; i <= limit; ++i){
gmp_printf("%3d! = %Zd\n", i, tmp);
if(i+1 <= limit) mpz_mul_ui(tmp, tmp, i+1);
}
return 0;
}
0! ~ 50!
50! ~ 100!
後記
上面只放一些常用的 function，另外有支援質數、最大公因數一些數學上的運算，C++ 的介紹和一些更深層的運用就留給大家慢慢發掘了。
GNU MP 6.2.1 manual
C++ Class Interface
Reference
如果你覺得這篇文章有用 可以考慮贊助飲料給大貓咪