The Bailey-Borwein-Plouffe formula provides a spigot algorithm for the computation of the nth binary digit of π. This summation formula was discovered in 1995 by Simon Plouffe. The formula is named after the authors of the paper in which the formula was first published, David H. Bailey, Peter Borwein, and Plouffe. The discovery of this formula came as a surprise. For centuries it had been assumed that there was no way to compute the nth digit of π without calculating all of the preceding n-1 digits.
The reason this pi formula is so interesting is because it can be used to calculate the N-th digit of Pi (in base 16) without having to calculate all of the previous digits!
Moreover, one can even do the calculation in a time that is essentially linear in N, with memory requirements only logarithmic in N. This is far better than previous algorithms for finding the N-th digit of Pi, which required keeping track of all the previous digits!
The algorithm is the fastest way to compute the nth digit (or a few digits in a neighborhood of the nth), but π-computing algorithms using large data types remain faster when the goal is to compute all the digits from 1 to n.
Original C Code For BBP.
/*
This program implements the BBP algorithm to generate a few hexadecimal
digits beginning immediately after a given position id, or in other words
beginning at position id + 1. On most systems using IEEE 64-bit floating-
point arithmetic, this code works correctly so long as d is less than
approximately 1.18 x 10^7. If 80-bit arithmetic can be employed, this limit
is significantly higher. Whatever arithmetic is used, results for a given
position id can be checked by repeating with id-1 or id+1, and verifying
that the hex digits perfectly overlap with an offset of one, except possibly
for a few trailing digits. The resulting fractions are typically accurate
to at least 11 decimal digits, and to at least 9 hex digits.
*/
/* David H. Bailey 2006-09-08 */
#include
#include
main()
{
double pid, s1, s2, s3, s4;
double series (int m, int n);
void ihex (double x, int m, char c[]);
int id = 1000000;
#define NHX 16
char chx[NHX];
/* id is the digit position. Digits generated follow immediately after id. */
s1 = series (1, id);
s2 = series (4, id);
s3 = series (5, id);
s4 = series (6, id);
pid = 4. * s1 - 2. * s2 - s3 - s4;
pid = pid - (int) pid + 1.;
ihex (pid, NHX, chx);
printf (" position = %i\n fraction = %.15f \n hex digits = %10.10s\n",
id, pid, chx);
}
void ihex (double x, int nhx, char chx[])
/* This returns, in chx, the first nhx hex digits of the fraction of x. */
{
int i;
double y;
char hx[] = "0123456789ABCDEF";
y = fabs (x);
for (i = 0; i < y =" 16." s =" 0.;" k =" 0;" ak =" 8" p =" id" t =" expm" s =" s" s =" s">= id. */
for (k = id; k <= id + 100; k++){ ak = 8 * k + m; t = pow (16., (double) (id - k)) / ak; if (t < s =" s" s =" s" expm =" 16^p" tp1 =" 0;" tp1 ="=" tp1 =" 1;" i =" 1;" ak ="=" i =" 0;"> p) break;
pt = tp[i-1];
p1 = p;
r = 1.;
/* Perform binary exponentiation algorithm modulo ak. */
for (j = 1; j <= i; j++){ if (p1 >= pt){
r = 16. * r;
r = r - (int) (r / ak) * ak;
p1 = p1 - pt;
}
pt = 0.5 * pt;
if (pt >= 1.){
r = r * r;
r = r - (int) (r / ak) * ak;
}
}
return r;
}
Results:
Position - Hex Digits At the Beginning
10^6 - 26C65E52CB4593
10^7 - 17AF5863EFED8D
10^8 - ECB840E21926EC
10^9 - 85895585A0428B
10^10 - 921C73C6838FB2
10^11 - 9C381872D27596
1.25×10^12 - 07E45733CC790B
2.5×10^14 - E6216B069CB6C1
Links:
http://mathworld.wolfram.com/BBPFormula.html
http://en.wikipedia.org/wiki/BBP-type_formula
No comments:
Post a Comment