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.

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

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
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.