[分享] 計算π到小數點下一億位

看板C_and_CPP (C/C++)作者Schottky (Schottky)時間12年前 (2013/09/18 18:18)推噓4(4推 0噓 0→)

留言4則, 4人參與討論串1/2 (看更多)

其實上一篇「計算 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