Question: I'm looking for the fastest way to obtain the value of π, as a personal
challenge. More specifically I'm using ways that don't involve using
Solution:
Bellard's formula, as used by PiHex, the now-completed distributed computing project, is used to calculate the nth digit of π in base 2. It is a faster version (about 43% faster) of the Bailey–Borwein–Plouffe formula.
#define
d constants like M_PI
, or hard-coding the number in.Solution:
Bellard's formula, as used by PiHex, the now-completed distributed computing project, is used to calculate the nth digit of π in base 2. It is a faster version (about 43% faster) of the Bailey–Borwein–Plouffe formula.
/** * Prints the nth number of pi followed by the next 8 numbers in base 10. * This program is based on Bellard's work. */ public class Bpp { final static int NUM = 990; // nth number of pi to print out /** * Returns the nth digit of pi followed by the next 8 numbers * @param n - nth number of pi to return * @return returns an integer value containing 8 digits after n */ public int getDecimal(long n) { long av, a, vmax, N, num, den, k, kq, kq2, t, v, s, i; double sum; N = (long) ((n + 20) * Math.log(10) / Math.log(2)); sum = 0; for (a = 3; a <= (2 * N); a = nextPrime(a)) { vmax = (long) (Math.log(2 * N) / Math.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 = mulMod(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 = mulMod(den, t, av); kq2 += 2; if (v > 0) { t = modInverse(den, av); t = mulMod(t, num, av); t = mulMod(t, k, av); for (i = v; i < vmax; i++) t = mulMod(t, a, av); s += t; if (s >= av) s -= av; } } t = powMod(10, n - 1, av); s = mulMod(s, t, av); sum = (sum + (double) s / (double) av) % 1; } return (int) (sum * 1e9); // 1e9 is 9 decimal places } private long mulMod(long a, long b, long m) { return (long) (a * b) % m; } private long modInverse(long a, long n) { long i = n, v = 0, d = 1; while (a > 0) { long t = i / a, x = a; a = i % x; i = x; x = d; d = v - t * x; v = x; } v %= n; if (v < 0) v = (v + n) % n; return v; } private long powMod(long a, long b, long m) { long tempo; if (b == 0) tempo = 1; else if (b == 1) tempo = a; else { long temp = powMod(a, b / 2, m); if (b % 2 == 0) tempo = (temp * temp) % m; else tempo = ((temp * temp) % m) * a % m; } return tempo; } private boolean isPrime(long n) { if (n == 2 || n == 3) return true; if (n % 2 == 0 || n % 3 == 0 || n < 2) return false; long sqrt = (long) Math.sqrt(n) + 1; for (long i = 6; i <= sqrt; i += 6) { if (n % (i - 1) == 0) return false; else if (n % (i + 1) == 0) return false; } return true; } private long nextPrime(long n) { if (n < 2) return 2; if (n == 9223372036854775783L) { System.err.println("Next prime number exceeds Long.MAX_VALUE: " + Long.MAX_VALUE); return -1; } for (long i = n + 1;; i++) if (isPrime(i)) return i; } /** * Runs the program * @param args */ public static void main(String args[]) { long duration = System.currentTimeMillis(); Bpp bpp = new Bpp(); System.out.println("Decimal digits of pi at position " + NUM + ": " + bpp.getDecimal(NUM) + "\n"); duration = System.currentTimeMillis() - duration; System.out.println("> " + duration + " ms"); } }
No comments:
Post a Comment