7
$\begingroup$

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.

  • 3
    I take it when you write "digit" you really mean "bit". We are considering binary, not decimal, expansions, right?2012-02-06

1 Answers 1

8

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