推 damody:同學你小時候的夢想好棒~ 09/18 20:56
推 lc85301:有點猛啊LOL 09/18 20:56
推 suhorng:推~ 現在世界紀錄也是用Chudnovsky formula算 09/18 21:22
推 isaacting:我覺得這比背圓周率到一萬位更有意義 !!!!推 09/18 21:46
其實上一篇「計算 e 到一億位」是為這一篇鋪路,
用的演算法大致相同, 只是公式較複雜。
原始碼: http://codepad.org/vGYg1YaD
Compile 這段程式需要 gnu mp 及 gnu mpfr 兩個程式庫,
Linux 可以直接安裝各 distribution 提供的程式開發套件, 這兩項都是標準配備,
Dev C++ 搭配 MinGW32 的版本可以下載 MinGW 官方網站的套件:
MinGW GMP: http://sourceforge.net/projects/mingw/files/MinGW/Base/gmp/
MinGW MPFR: http://sourceforge.net/projects/mingw/files/MinGW/Base/mpfr/
在 64-bit Linux, 2.4GHz CPU 上計算一億位大約需要 2GB RAM, 時間約六分鐘。
這次採用 GMP 的整數運算和 MPFR 的浮點運算,
MPFR 比起 GMP 的 mpf 浮點運算有幾個優點:
1. 精確度 (precision) 定義明確, 使用 32-bit 和 64-bit library 不會在尾數
有差異, 而且可以從 MPFR_PREC_MAX 知道 precision 最大值, 不必用猜的。
2. 尾數的取捨有明確定義, 是要無條件捨去還是四捨五入都要在運算時指明。
3. 轉換成小數字串時, GMP 會自動省略尾數的 0, MPFR 不會。當計算π到 50 位
時, 最後一位剛好就是 0, 少一位會讓人懷疑是不是偷工減料。
4. 豐富完整的數學運算函數, 如 log, 非整數次方, 三角函數等等。
MPFR 的缺點則是:
1. MPFR 效能比 GMP 慢 10% 左右, 但目前這個程式大部份時間並非耗在浮點運算,
而是在遞迴函數內的整數乘法, 所以影響不大。
2. 多一個 library 要 link, 多幾行安裝說明要寫, 有點不方便。 (還是寫了)
3. 最大精確度只有 GMP 的一半, 但沒有人有這麼大的 RAM 所以也沒差。
整體來說我還是比較喜歡 MPFR, 因為設計得很細心, 對極限細微操作的掌握較好。
以這篇來說, 不會有人跑來問我輸入 50 位為什麼只出現 49 位小數...
http://en.wikipedia.org/wiki/Pi#Rapidly_convergent_series
用來計算π的是 Chudnovsky 公式, 維基百科寫每多計算一項增加 14 位有效小數,
我推導出來是當 N 值越大, 每項增加的位數趨近 3*log10(53360) = 14.1816474627
所以可以直接用除法計算數列共需計算幾項, 才能到達使用者要求的精確度。
R(a,b)*P(a,b)
定義 S(a,b) = -------------
Q(a,b)
(6b)! (6a)!
R(a,b) = ----- ÷ ----- (第 a+1 項到第 b 項提出來的乘數)
(3b)! (3a)!
b!
Q(a,b) = (----)^3*(640320)^(3*(b-a)) (第 a+1 項到第 b 項的分母)
a!
P(a,b) = 第 a+1 項到第 b 項的 (13591409+545140134*k) 分子部份通分後的和
剩下的就和上一篇計算 e 的級數一樣, 使用遞迴函數做 Divide and Conquer 。
為什麼是第 a+1 項到第 b 項? 這有點難解釋, 寫過一遍會比較容易理解...
這樣定義後在 merge 時可以省掉一個小乘法運算, 銜接處比較簡潔。
這次使用一模一樣的進度條。進度條的作用並不是逗貓或讓無聊的使用者打發時間,
而是讓我自己知道程式有沒有在跑? 大概還要跑多久? 藉此判斷是不是該喊卡了。
所以只要程式執行時間會超過一秒, 大概都會需要個進度條或里程表。
進度條並不難做, 短短 20 行左右的一個小函式就行了, 而且很容易重複使用。
另一個計算π的公式是 Gauss-Legendre 法, 這個數列無法使用 Divide & Conquer,
只能依序計算。它的收斂能力非常高, 23 項就有一千萬位, 27 項就有一億位。
每多算一項, 就會產生兩倍的有效數字, 而且使用 MPFR 來寫, 程式可以很簡潔。
我對高斯法的細節還沒有深入研究, 各位如果有興趣可以 implement 看看 :)
// URL: http://codepad.org/vGYg1YaD
// -------- >8 -------- >8 -------- >8 -------- >8 -------- >8 --------
#include <stdio.h>
#include <stdlib.h>
#include <gmp.h>
#include <mpfr.h>
#include <math.h>
#define FILENAME "pi.txt"
unsigned long long int count, total, old_p;
void progress(void)
{
int p, g, c;
p = (int)floorf(1000.0f*count/total+0.5f);
if (p > old_p) {
old_p = p;
g = (int)floorf(72.0f*count/total+0.5f);
for (c=0; c<72; c++) {
if (c<g) printf("H");
else printf("-");
}
printf(" %5.1f%%\r", p/10.0);
fflush(stdout);
}
if (count == total) {
printf("\n");
}
}
void write_pi(mpfr_t pi, int digits)
{
char *pi_str;
mp_exp_t exp;
int c;
FILE *fp;
pi_str = malloc(digits+1);
fp = fopen(FILENAME, "w+");
mpfr_get_str(pi_str, &exp, 10, digits, pi, MPFR_RNDZ);
fprintf(fp, "pi=");
for (c=0; c<digits; c++) {
if (c!=1) {
if (c%50 == 1) {
fprintf(fp, " ");
}
}
fprintf(fp, "%c", pi_str[c]);
if (c==0) {
fprintf(fp, ".");
} else if (c%1000 == 0) {
fprintf(fp, " << %d\n", c);
} else if (c%500 == 0) {
fprintf(fp, " <\n");
} else if (c%50 == 0) {
fprintf(fp, "\n");
} else if (c%5 == 0) {
fprintf(fp, " ");
}
}
if (c%50 != 1) {
fprintf(fp, "\n");
}
fclose(fp);
free(pi_str);
}
void s(mpz_t a, mpz_t b, mpz_t max, mpz_t p, mpz_t q, mpz_t r)
{
mpz_t m, p1, q1, r1;
mpz_inits(m, p1, q1, r1, NULL);
mpz_sub(m, b, a);
if (mpz_cmp_ui(m, 0) <= 0) {
gmp_printf("Error! s(%Zd, %Zd).\n", a, b);
} else if (mpz_cmp_ui(m, 1) == 0) {
if (mpz_cmp_ui(a, 0) == 0) {
mpz_ui_pow_ui(q, 640320, 3);
mpz_set_ui(r, 120); // 6!
mpz_mul_ui(p, q, 13591409);
mpz_submul_ui(p, r, 13591409+545140134);
} else {
mpz_set_ui(r, 8);
mpz_mul_ui(m, a, 6);
mpz_add_ui(m, m, 1);
mpz_mul(r, r, m); // r *= 6a+1;
mpz_add_ui(m, m, 2);
mpz_mul(r, r, m); // r *= 6a+3;
mpz_add_ui(m, m, 2);
mpz_mul(r, r, m); // r *= 6a+5;
mpz_mul_ui(q, b, 640320);
mpz_pow_ui(q, q, 3);
mpz_set_ui(p, 13591409);
mpz_addmul_ui(p, b, 545140134);
mpz_mul(p, p, r);
if (mpz_tstbit(b, 0) == 1) { mpz_neg(p, p); }
}
} else {
mpz_add(m, a, b);
mpz_fdiv_q_ui(m, m, 2);
s(a, m, max, p, q, r);
s(m, b, max, p1, q1, r1);
// Merge
mpz_mul(p, p, q1);
mpz_addmul(p, p1, r);
mpz_mul(q, q, q1);
if (mpz_cmp(b, max) == 0) {
mpz_realloc2(r, 1);
} else {
mpz_mul(r, r, r1);
}
}
// gmp_printf("s(%Zd, %Zd) = %Zd / %Zd (%Zd)\n", a, b, p, q, r);
count++;
progress();
mpz_clears(m, p1, q1, r1, NULL);
}
void chudnovsky(unsigned long int digits)
{
unsigned long int d, n, precision;
mpz_t max, p, q, r, z;
mpfr_t pf, qf;
d = digits+1;
n = ceill(d*logl(10)/(3*logl(53360)));
precision = ceill(d*logl(10.0)/logl(2.0));
count = old_p = 0;
total = n*2-1;
printf("Digits = %lu, N = %lu, Precision = %lu\n",
d, n, precision);
mpz_inits(max, p, q, r, z, NULL);
mpfr_inits2(precision, pf, qf, NULL);
mpz_set_ui(max, n);
mpz_set_ui(z, 0);
s(z, max, max, p, q, r);
mpz_mul_ui(q, q, 426880);
mpfr_set_z(pf, p, MPFR_RNDN);
mpfr_set_z(qf, q, MPFR_RNDN);
mpfr_div(pf, qf, pf, MPFR_RNDN);
mpfr_sqrt_ui(qf, 10005, MPFR_RNDN);
mpfr_mul(pf, pf, qf, MPFR_RNDN);
mpz_clears(max, p, q, r, z, NULL);
printf("Calculation done. Writing result to file <%s>...\n", FILENAME);
write_pi(pf, d);
mpfr_clears(pf, qf, NULL);
}
int main(int argc, char *argv[])
{
mpfr_t digits, max_digit;
mpfr_inits2(70, digits, max_digit, NULL);
// max_gitit = MPFR_PREC_MAX*log10(2)
mpfr_set_ui(max_digit, 2, MPFR_RNDN);
mpfr_log10(max_digit, max_digit, MPFR_RNDN);
mpfr_mul_ui(max_digit, max_digit, MPFR_PREC_MAX, MPFR_RNDN);
mpfr_sub_ui(max_digit, max_digit, 1, MPFR_RNDN);
mpfr_floor(max_digit, max_digit);
if (argc == 2) {
mpfr_set_str(digits, argv[1], 10, MPFR_RNDZ);
mpfr_floor(digits, digits);
if (mpfr_cmp_ui(digits, 1) < 0 || mpfr_cmp(digits, max_digit) > 0) {
mpfr_fprintf(stderr, "Invalid parameter, "
"digits must range from 1 to %.0Rf\n", max_digit);
return -1;
}
} else if (argc > 2) {
fprintf(stderr, "Usage: %s #digits\n", argv[0]);
return -1;
} else {
mpfr_set_ui(digits, 100000000, MPFR_RNDN);
}
chudnovsky(mpfr_get_ui(digits, MPFR_RNDZ));
mpfr_clears(digits, max_digit, NULL);
return 0;
}
--
※ 發信站: 批踢踢實業坊(ptt.cc)
◆ From: 118.167.24.180