其实上一篇「计算 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