BBP functions (Bailey-Borwein-Plouffe) =====Warren D Smith====Dec 2011=== I seem to have made another mistake (will be explained below). In this post, I will review "BBP functions" hopefully more carefully this time... BBP1: ln(3) = sum(n>=1) 16^(-n) * (16/(4*n-3) + 4/(4*n-1)) pi = sum(n>=0) (-4)^(-n) * (2/(4*n+1) + 1/(2*n+1) + 1/(4*n+3)) are examples of "BBP formulas." They enable computing the Nth digit of ln(3) in radix 16, or the Nth digit of pi in radix 4, WITHOUT computing the preceding N-1 digits... and the whole algorithm takes roughly O(N) operations using O(logN) bits of memory. This is essentially because there is a fast way [based on fast modular powering] to figure out the digits in the vicinity of the Nth in the decimal (actually radix=R) expansion of any rational number p/q with p,q small. BBP2: Simon Plouffe pointed out that he invented a weaker kind of BBP-like formula which then got explored & improved further by Fabrice Bellard & Xavier Gourdon: http://arxiv.org/abs/0912.0303 http://bellard.org/pi/pi_n2/pi_n2.html http://numbers.computation.free.fr/Constants/Algorithms/nthdecimaldigit.pdf Examples of this kind of formula are pi+3 = sum(n>=1) n * 2^n / binomial(2*n, n) zeta(3) = (5/2) * sum(n>=1) (-1)^(n-1) / (n^3 * binomial(2*n, n)). In these formulas, the rational numbers p/q involved have LARGE q, but this is ok because q is factorable into small primes with no large powers of a prime, which allows the whole rational p/q can be rewritten as a sum of rationals with small denominators. As a result, these authors showed how to evaluate the Nth digit of pi+3 now in ANY fixed radix in roughly O(N^2) time using O(logN)^2 bits of memory. The BBP2 formulas are of less interest than BBP1 because they actually run much SLOWER to compute the Nth digit than algorithms which compute ALL digits 1,2,...,N. (You could argue they still have some interest because of their tiny memory needs and since they now work in every radix...) BBP3: Another weaker kind of BBP formula would be exemplified by F(x) = -ln(1-x)/(x*(1-x)) = sum(n>=0) H[n+1] * x^n where H[n] = 1 + 1/2 + 1/3 + ... + 1/n is the nth "Harmonic number." By regarding the series as a DOUBLE sum we can compute the Nth decimal digit of F(1/10) in roughly O(N^2) operations and memory O(logN). Similarly one could consider triple sums, etc; then we would evaluate the Nth digit in polynomial time O(N^S) operations roughly, and O(logN) bits of memory, for an S-fold sum for any fixed S. Again, I consider these BBP3 formulas of less interest than BBP1. However, in at least some cases I think you can "cheat" to evaluate the Nth digits of these things effectively with a smaller value of S, i.e. faster runtime than the naive bound. If so, in such cases they would regain interest, so it would be good to figure out just when this happens. BBP4: I also pointed out to Plouffe that if you are willing to consume "quasipolynomial" compute time N^O(D) where D can now depend logarithmically on N -- or even worse like a power of log(N) -- then you can compute the Nth digit of a much wider class of reals using log(N)^O(1) bits of memory. This works in any radix and for a far larger class of real numbers. I consider this result to be of even less interest than BBP2 and BBP3 because it is even slower. Nevertheless (since Simon Plouffe told me he did not believe me), I will explain how this is done. Let X be ANY real number which is "D(N)-shallow" meaning it can be computed to N-digit accuracy starting from 0 and 1 by an D(N)-depth "arithmetic circuit" each of whose nodes adds, subtracts, or multiplies its two inputs and outputs the result. There are also certain other operations nodes could use, such as, the nodes could only look at certain subranges of the digits of their inputs, or could keep track of the number of digits the number has if it is an integer (these go beyond the 3 arithmetic operations strictly speaking, but still would be allowed). Then I claim the Nth digit of X can be computed by an algorithm taking roughly N^D(N) time and running in O(D(N)*logN) bits of memory. The method is simple. To compute the Nth digit of X -- actually a O(logN)-wide block of digits starting at position N, and we shall assume no "carries" propagate further than distance O(logN), an assumption which the previous BBP algorithms also needed to make: * if X=A+B then compute the digit blocks for A and B and add them. * If X=A-B ditto. * If X=A*B then simply scan over k computing the kth digit of A and the (N-k)th digit of B, recursively, and accumulate products. The recursion "stack" gets to be O(D(N))-deep and each recursion level you need to keep track of a counter (we called it k) using O(logN) bits of memory. Example: e=1+1/2!+1/3!+... We can compute the Nth digit of e in O(logN)^2 bits of memory and time N^O(logN) as follows. We need to show that e can be computed to N-digit accuracy by an O(logN)-deep arithmetic circuit. The divisions will be computed using +, -, * only via Newton iteration in O(logN) depth to get N digit accuracy. The sum with N terms will be computed in only O(logN) depth using binary-tree summation. The factorials will be computed in only O(logN) depth by using * if N odd: N! = N * (N-1)! * if N even: N! = (N/2)!^2 * binomial(N,N/2) and get digits of binomial(N,N/2) inside high powers by using Newton binomial formula and using binary powering of integers like 1000001^N [these numbers will be O(N^2) digits long but at we shall only employ at O(N)-long subranges of them; but even using the full length would be fine]. QED. Further examples: A. Every hypergeometric function value pFq([a1,a2,...,ap],[b1,b2,...,bq]; x) where x=1/radix or x=integer and q>p (or q=p for x=1/radix) is O(logN)-shallow. B. For general integer a,b every hypergeometric function value pFq([a1,a2,...,ap],[b1,b2,...,bq]; a/b) is O(logN)^2-shallow. C. Any rational power (a/b)^(p/q) of a rational, or any algebraic number, is O(logN)-shallow. D. Composing functions basically adds their depths. LARGE CLASSES OF FUNCTIONS THAT ARE BBP: I. If R(x) is any rational function of x with integer coefficients and whose poles lie on the unit circle |x|=1, then R(x) is BBP1 whenever x=1/radix. Proof sketch: expand R(x) using partial fractions, find the Maclaurin series for each term, and note its coefficients of x^n are rational(n). Also if R(x) has poles on circles |x|=r where r is a power of the radix, this also is ok -- but that specializes things to a particular radix, whereas with r=1 only, R(x) will preserve BBPness for general radix. II. Integration: integral F(x) dx is BBP1 if F was. Proof sketch: c[n]*x^n in the Maclaurin series becomes c[n]*x^(n+1)/(n+1). III. Derivative: F'(x) is BBP1 if F was. IV. Sum: (a/b)*F(x)+(c/d)*G(x) is BBP1 if F and G are and a,b,c,d are integers. V. Product: F(x)*G(x) is BBP3 if F and G were BBP1 or BBP3. Note there is a weakening here: we fall into the weaker BBP-like function class BBP3 not the most interesting class BBP1. (That was my error before.) VI. Powers of argument: F(x^Q) is BBP1 if Q is positive integer and F is BBP1. In the above it is ok to replace all "BBP1" by "BBP3." By performing operations II, III, IV, V, VI starting from functions in I we get a very large set of BBP3 functions F(x). If x=a/b is rational where b is a power of the radix, then the Nth digit of F(x) can be computed in polynomial time with O(logN) bits of memory. I think (but have not carefully checked) that every BBP formula so far compiled in Bailey's huge compendium, lies within a tiny subclass of the above. -- Warren D. Smith http://RangeVoting.org <-- add your endorsement (by clicking "endorse" as 1st step)