我一直想找到一个能做到这一点的算法.我不在乎它有多慢,只要它可以返回Pi的第n位:
例如:
size_t piAt(long long int n) { }
优选地,不使用无限系列.
如果有人有一个函数或类来做这个,在C或C++中,我真的很有兴趣看到它.
谢谢
这个非凡的解决方案展示了如何在O(N)时间和O(log·N)空间中计算π的第 N 位数,并且无需计算导致它的所有数字.
哦,它是十六进制的.
如果你不想这样做,你可以很容易地从shell中做到这一点:
% perl -Mbignum=bpi -wle 'print bpi(20)' 3.1415926535897932385 % perl -Mbignum=bpi -wle 'print bpi(50)' 3.1415926535897932384626433832795028841971693993751 % perl -Mbignum=bpi -wle 'print bpi(200)' 3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303820 % perl -Mbignum=bpi -wle 'print bpi(1000)' 3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117067982148086513282306647093844609550582231725359408128481117450284102701938521105559644622948954930381964428810975665933446128475648233786783165271201909145648566923460348610454326648213393607260249141273724587006606315588174881520920962829254091715364367892590360011330530548820466521384146951941511609433057270365759591953092186117381932611793105118548074462379962749567351885752724891227938183011949129833673362440656643086021394946395224737190702179860943702770539217176293176752384674818467669405132000568127145263560827785771342757789609173637178721468440901224953430146549585371050792279689258923542019956112129021960864034418159813629774771309960518707211349999998372978049951059731732816096318595024459455346908302642522308253344685035261931188171010003137838752886587533208381420617177669147303598253490428755468731159562863882353787593751957781857780532171226806613001927876611195909216420199
以下是Fabrice Bellard编写的Simon Plouffe解决方案:
/* * Computation of the n'th decimal digit of \pi with very little memory. * Written by Fabrice Bellard on January 8, 1997. * * We use a slightly modified version of the method described by Simon * Plouffe in "On the Computation of the n'th decimal digit of various * transcendental numbers" (November 1996). We have modified the algorithm * to get a running time of O(n^2) instead of O(n^3log(n)^3). * * This program uses mostly integer arithmetic. It may be slow on some * hardwares where integer multiplications and divisons must be done * by software. We have supposed that 'int' has a size of 32 bits. If * your compiler supports 'long long' integers of 64 bits, you may use * the integer version of 'mul_mod' (see HAS_LONG_LONG). */ #include#include #include /* uncomment the following line to use 'long long' integers */ /* #define HAS_LONG_LONG */ #ifdef HAS_LONG_LONG #define mul_mod(a,b,m) (( (long long) (a) * (long long) (b) ) % (m)) #else #define mul_mod(a,b,m) fmod( (double) a * (double) b, m) #endif /* return the inverse of x mod y */ int inv_mod(int x, int y) { int q, u, v, a, c, t; u = x; v = y; c = 1; a = 0; do { q = v / u; t = c; c = a - q * c; a = t; t = u; u = v - q * u; v = t; } while (u != 0); a = a % y; if (a < 0) a = y + a; return a; } /* return (a^b) mod m */ int pow_mod(int a, int b, int m) { int r, aa; r = 1; aa = a; while (1) { if (b & 1) r = mul_mod(r, aa, m); b = b >> 1; if (b == 0) break; aa = mul_mod(aa, aa, m); } return r; } /* return true if n is prime */ int is_prime(int n) { int r, i; if ((n % 2) == 0) return 0; r = (int) (sqrt(n)); for (i = 3; i <= r; i += 2) if ((n % i) == 0) return 0; return 1; } /* return the prime number immediatly after n */ int next_prime(int n) { do { n++; } while (!is_prime(n)); return n; } int main(int argc, char *argv[]) { int av, a, vmax, N, n, num, den, k, kq, kq2, t, v, s, i; double sum; if (argc < 2 || (n = atoi(argv[1])) <= 0) { printf("This program computes the n'th decimal digit of \\pi\n" "usage: pi n , where n is the digit you want\n"); exit(1); } N = (int) ((n + 20) * log(10) / log(2)); sum = 0; for (a = 3; a <= (2 * N); a = next_prime(a)) { vmax = (int) (log(2 * N) / log(a)); av = 1; for (i = 0; i < vmax; i++) av = av * a; s = 0; num = 1; den = 1; v = 0; kq = 1; kq2 = 1; for (k = 1; k <= N; k++) { t = k; if (kq >= a) { do { t = t / a; v--; } while ((t % a) == 0); kq = 0; } kq++; num = mul_mod(num, t, av); t = (2 * k - 1); if (kq2 >= a) { if (kq2 == a) { do { t = t / a; v++; } while ((t % a) == 0); } kq2 -= a; } den = mul_mod(den, t, av); kq2 += 2; if (v > 0) { t = inv_mod(den, av); t = mul_mod(t, num, av); t = mul_mod(t, k, av); for (i = v; i < vmax; i++) t = mul_mod(t, a, av); s += t; if (s >= av) s -= av; } } t = pow_mod(10, n - 1, av); s = mul_mod(s, t, av); sum = fmod(sum + (double) s / (double) av, 1.0); } printf("Decimal digits of pi at position %d: %09d\n", n, (int) (sum * 1e9)); return 0; }
它有效:
C:\tcc>tcc -I libtcc libtcc/libtcc.def -run examples/pi.c 1 Decimal digits of pi at position 1: 141592653 C:\tcc>tcc -I libtcc libtcc/libtcc.def -run examples/pi.c 1000 Decimal digits of pi at position 1000: 938095257
http://bellard.org/pi/pi_n2/pi_n2.html