[9/17 12:15pm 更新: 增加位数范围检查, 修正 terms() 问题。] 这是为了完成一个小时候的梦想... 计算 e 到小数点下一亿位... ∞ 所使用的级数是 e = Σ(1/k!) = 1/0! + 1/1! + 1/2! + 1/3! + 1/4! + ..... k=0 这个级数看似简单, 其实很优秀, 计算不难, 且收敛速度非常快。 程式码: http://codepad.org/qtZaidZq 这次一样用二分法, 把 N 项的级数切成一半, 左右分别算完再合并结果, 要注意的是需要计算的总项数和分数通分加法数量和回圈法其实是一样多的, 那为什麽递回 Divide & Conquer 会比回圈法快这麽多呢? 是因为上次我打马虎眼被 Feis 抓包的大数运算复杂度并不是 O(1) ... 二分法切割到最後都是小小的运算, 两项两项小数字合并, 速度就很快, 回圈法跑没多久, 加总的数字开始变得很大, 就举步维艰了... 不要被 Divide & Conquer 这个名词骗了, 重点在 Conquer 之後的 merge, 只要能找到有效率的 merge 公式就成功了... 幸运地, 在这边是有的, 把计算结果用分数表示 (分子&分母分别用整数表示) 再通分即可合并。 千万不要全部转成浮点数计算, 每一次运算都用小数点以下一亿位的超高精度浮点 运算会非常慢... 整数的话只有最後的几次超级合并才有这麽大的数量级... 而且化成分数最後只进行一次超级除法, 可以把浮点误差降到最低。 一开始需要估计计算到第几项才能达到我们要的精度。 因为不知道怎麽找出 Stirling's approximation 的反函数, 我以二分逼近法求解。 在 64-bit Linux, 2.4GHz CPU 上面大约需要两分半钟, 所以我写了进度条, 其中有一分钟是将最後 120MB 左右的计算结果转换成十进位字串, 以及写到输出档。 如果一开始就使用十进位(十亿进位)来储存和运算, 最後的转换过程就可以省下来。 我在 Linux 和 MinGW 上都 compile 过, 32-bit 环境也可以跑, 没问题的。 输入位数最大可以到 2^64-1, 问题是 RAM 不够, 一亿(10^8)位就要 2GB RAM 了, 估计要计算十亿位需要 20GB RAM, 各位如果手边有暴力机台可以测试看看 :) // URL: http://codepad.org/qtZaidZq // -------- >8 -------- >8 -------- >8 -------- >8 -------- >8 -------- #include <stdio.h> #include <stdlib.h> #include <string.h> #include <limits.h> #include <gmp.h> #include <math.h> #define FILENAME "e.txt" #define PI 3.141592653589793238462643L #define E 2.718281828459045235360287L unsigned long long int count=0, total=0, old_p=0; long double logfactorial(unsigned long long int n) { return (logl(2.0*PI)/2.0+logl(n)/2.0+n*logl(n/E))/logl(10.0); } unsigned long long int terms(unsigned long long int digits) { long double f; unsigned long long int n, upper, lower; for (n=1;;n<<=1) { f = logfactorial(n); if (f>digits) break; } upper = n; lower = 0; while (upper-lower > 1) { n = (upper+lower)/2; f = logfactorial(n); if (f>digits) { upper = n; } else { lower = n; } } n = upper; return n; } void write_e(mpf_t e, unsigned long long int digits) { char *str; long int exp; int c; FILE *fp; fp = fopen(FILENAME, "w+"); if (fp == NULL) { return; } str = malloc(digits+2); mpf_get_str(str, &exp, 10, digits+1, e); fprintf(fp, " e = "); for (c=0; c<digits; c++) { if (c!=1) { if (c%50 == 1) { fprintf(fp, "\t"); } } fprintf(fp, "%c", 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"); } free(str); fclose(fp); } 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.5f*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 calc(mpz_t a, mpz_t b, mpz_t p, mpz_t q) { mpz_t m, p2, q1; mpz_inits(m, p2, q1, NULL); mpz_sub(m, b, a); if (mpz_cmp_ui(m, 0) <= 0) { if (mpz_cmp_ui(a, 0) == 0) { mpz_set_ui(p, 1); } else { mpz_set_ui(p, 0); } mpz_set_ui(q, 1); } else if (mpz_cmp_ui(m, 1) == 0) { if (mpz_cmp_ui(a, 0) == 0) { mpz_set_ui(p, 2); mpz_set_ui(q, 1); } else { mpz_set_ui(p, 1); mpz_set(q, b); } } else { mpz_add(m, a, b); mpz_fdiv_q_ui(m, m, 2); calc(a, m, p, q1); calc(m, b, p2, q); // Merge mpz_mul(p, p, q); mpz_add(p, p, p2); mpz_mul(q, q, q1); } //gmp_printf("P(%d,%d)=%Zd, Q(%d,%d)=%Zd\n", a, b, p, a, b, q); count++; progress(); mpz_clears(m, p2, q1, NULL); } void calc_e(unsigned long long int digits) { mpz_t z, p, q, n; mpf_t pf, qf; unsigned long long int d, n_terms, precision; if (digits < ULLONG_MAX) { d = digits + 1; } else { d = ULLONG_MAX; } n_terms = terms(d); total = n_terms*2-1; precision = ceill(d*(logl(10.0L)/logl(2.0L)))+1; printf("d = %llu, n = %llu, precision = %llu\n", d, n_terms, precision); mpf_set_default_prec(precision); precision = mpf_get_default_prec(); printf("Real precision = %llu\n", precision); mpz_inits(z, p, q, n, NULL); mpf_inits(pf, qf, NULL); mpz_set_ui(z, 0); mpz_set_ui(n, n_terms); calc(z, n, p, q); mpf_set_z(pf, p); mpf_set_z(qf, q); mpf_div(pf, pf, qf); mpz_clears(z, p, q, n, NULL); printf("Calculation done. Writing result to file...\n"); write_e(pf, d); mpf_clears(pf, qf, NULL); } int main(int argc, char *argv[]) { unsigned long long int digits = 100000000ULL; unsigned long long int prec_max; // I don't know how GNU MP sets it's maximum precision, it's black magic. if (ULONG_MAX == 4294967295UL) { prec_max = 1292913975UL; // 32-bit } else { prec_max = 5553023288523357111ULL; // 64-bit } if (argc == 2) { digits = strtoull(argv[1], NULL, 10); // 1292913985 = MAX_PRECISION*log10(2)-1 , MAX_PRECISION = 2^32-1 if (digits == 0 || digits > prec_max) { fprintf(stderr, "Invalid parameter, " "digits must range from 1 to %llu\n", prec_max); return -1; } } else if (argc > 2) { fprintf(stderr, "Usage: %s #digits\n", argv[0]); return -1; } calc_e(digits); return 0; } --

※ 发信站: 批踢踢实业坊(ptt.cc) ◆ From: 220.137.9.131 感谢提醒, terms() 的确是疏忽了。 这里用的是 GNU MP 不是 MPFR, 原因是 mpf_ 系列运算稍微快一点点, 而且都已经用了 gmp, 我不想再多写一份 mpfr 安装说明... 事实上 gmp/mpfr 两个版本我都有写, 所以差不多可以另外写一篇评比了 =.=" 「硬干」的意思是? 如果是指 implement 大数运算, 那整个程式码会膨胀很多, 但想要突破 RAM 的限制把 HD 拿来当计算纸, gmp 是办不到的, 只能自己写。 目前世界纪录 10^13 位数就是这样做的, HD 当 RAM 用, RAM 当 cache 用。 ※ 编辑: Schottky 来自: 220.137.7.154 (09/17 12:33)