# Using Bellard's formula for calculating $\pi$?

I don't have any experience with these types of formulae and am finding it difficult understanding how to use Bellards formula.

Say I wanted to get the $3^\text{rd}$ digit of $\pi$ (which equals 1), would I simply replace all the instances of $n$ in the formula for $1$, and then multiply by $\frac{1}{2^6}$ ?

If I wanted to get more than one digit, do I get them all without multiplying by $\frac{1}{ 2^6}$, sum and then multiply?

I have googled the problem extensively and the only resources that are available seem to explain how the formula was gotten and not its use.

-
I take it when you write "digit" you really mean "bit". We are considering binary, not decimal, expansions, right? –  Gerry Myerson Feb 6 '12 at 23:15

You should see how BBP formulas are used to extract the $p$-th digit.

The Bellard's formula, when truncated to $m$ terms will give you a rational number $r_m$, approximating $\pi$.

Suppose we are interested in finding $p^\text{th}$ binary digit of $\pi$. Choose $m$ is sufficiently large, so that $| \pi -r_m | \leqslant 2^{-p}$. In order to extract the digit, we need to compute $d_p = \lfloor 2^p r_m \rfloor \bmod 2$.

Consider $$2^{p-1} r_m = \sum_{n=0}^\infty (-1)^{n} 2^{p-1-10n} \Big( \frac{4}{10 n+1}-\frac{1}{10 n+3}-\frac{1}{2 (4 n+1)}-\frac{1}{2^4 (10 n+5)}-\frac{1}{2^4 (10 n+7)}+\frac{1}{2^6 (10 n+9)}-\frac{1}{2^6 (4 n+3)} \Big) = \Sigma_1 - \Sigma_2 - \Sigma_3 - \Sigma_4 - \Sigma_5 + \Sigma_6 - \Sigma_7$$ Now the trick is not to compute the digits which precede the $p$-th one. To this end one splits the sum into 7 pieces. We consider just one: $$\Sigma_1 = \sum_{n=0}^\infty \frac{(-1)^{n} 2^{1 +p-10n} }{10n+1}$$ Not computing the preceding digits is equivalent to only retaining the fractional part. This is done by computing the numerator modulo the numerator: $$\sum_{n=0}^\infty \frac{(-1)^{n} 2^{1 +p-10n} \bmod (10 n+1) }{10n+1}$$

Here is a little Mathematica code for Bellard's algorithm:

Sigma[p_Integer, a_Integer, b_Integer, c_Integer, prec_] :=
Module[{sum = 0, prev, den, n = 0, num, exp,
bprec = 1 + prec Log2[10.]},
While[True,
den = a n + b;
exp = c + p - 1 - 10 n;
If[exp > 0,
num = Mod[PowerMod[2, exp, den], den, 1],
num = 1; den = den*2^(-exp);
];
prev = sum;
sum += Divide[(-1)^n SetPrecision[num, prec], den];
If[sum - prev == 0 && exp < -bprec, Break[]];
n++
];
sum
]

PthPiDigitBellard[p_Integer?NonNegative, prec_: MachinePrecision] :=
Floor[2 Mod[Sigma[p, 10, 1, 2, prec] - Sigma[p, 10, 3, 0, prec] -
Sigma[p, 4, 1, -1, prec] - Sigma[p, 10, 5, -4, prec] -
Sigma[p, 10, 7, -4, prec] + Sigma[p, 10, 9, -6, prec] -
Sigma[p, 4, 3, -6, prec], 1]]

In[166]:= {Table[PthPiDigitBellard[k, 20], {k, 1, 60}], 0} ==
RealDigits[Pi, 2, 60, -1]

Out[166]= True
-