RandomLib::InversePiProb Class Reference

Return true with probability 1/π. More...

#include <RandomLib/InversePiProb.hpp>

Public Member Functions
template<class Random >
bool	operator() (Random &r)

Detailed Description

Return true with probability 1/π.

InversePiProb p; p(Random& r) returns true with prob 1/π using the method of Flajolet et al. It consumes 9.6365 bits per call on average.

The method is given in Section 3.3 of

• P. Flajolet, M. Pelletier, and M. Soria,
  On Buffon Machines and Numbers,
  Proc. 22nd ACM-SIAM Symposium on Discrete Algorithms (SODA), Jan. 2011.
  http://www.siam.org/proceedings/soda/2011/SODA11_015_flajoletp.pdf
  

using the identity

\[ \frac 1\pi = \sum_{n=0}^\infty {{2n}\choose n}^3 \frac{6n+1}{2^{8n+2}} \]

It is based on the expression for 1/π given by Eq. (28) of

• S. Ramanujan,
  Modular Equations and Approximations to π,
  Quart. J. Pure App. Math. 45, 350–372 (1914);
  In Collected Papers, edited by G. H. Hardy, P. V. Seshu Aiyar, B. M. Wilson (Cambridge Univ. Press, 1927; reprinted AMS, 2000).
  http://books.google.com/books?id=oSioAM4wORMC&pg=PA36
  

\[\frac4\pi = 1 + \frac74 \biggl(\frac 12 \biggr)^3 + \frac{13}{4^2} \biggl(\frac {1\cdot3}{2\cdot4} \biggr)^3 + \frac{19}{4^3} \biggl(\frac {1\cdot3\cdot5}{2\cdot4\cdot6} \biggr)^3 + \ldots \]

The following is a description of how to carry out the algorithm "by hand" with a real coin, together with a worked example:

1. Perform three coin tossing experiments in which you toss a coin until you get tails, e.g., HHHHT; HHHT; HHT. Let h₁ = 4, h₂ = 3, h₃ = 2 be the numbers of heads tossed in each experiment.
2. Compute n = ⌊h₁/2⌋ + ⌊h₂/2⌋ + mod(⌊(h₃ − 1)/3⌋, 2) = 2 + 1 + 0 = 3. Here is a table of the 3 contributions to n:

   0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17  h
   0  0  1  1  2  2  3  3  4  4  5  5  6  6  7  7  8  8  floor(h1/2)
   0  0  1  1  2  2  3  3  4  4  5  5  6  6  7  7  8  8  floor(h2/2)
   1  0  0  0  1  1  1  0  0  0  1  1  1  0  0  0  1  1  mod(floor((h3-1)/3), 2)
   

3. Perform three additional coin tossing experiments in each of which you toss a coin 2n = 6 times, e.g., TTHHTH; HHTHH|H; THHHHH. Are the number of heads and tails equal in each experiment? yes and no and no → false. (Here, you can give up at the |.)

The final result in this example is false. The most common way a true result is obtained is with n = 0, in which case the last step vacuously returns true.

Proof of the algorithm: Flajolet et al. rearrange Ramanujan's identity as

\[ \frac 1\pi = \sum_{n=0}^\infty \biggl[{2n\choose n} \frac1{2^{2n}} \biggr]^3 \frac{6n+1}{2^{2n+2}}. \]

Noticing that

\[ \sum_{n=0}^\infty \frac{6n+1}{2^{2n+2}} = 1, \]

the algorithm becomes:

1. pick n ≥ 0 with prob (6n+1) / 2²ⁿ⁺² (mean n = 11/9);
2. return true with prob (binomial(2n, n) / 2²ⁿ)³.

Implement (1) as

• geom4(r) + geom4(r) returns n with probability 9(n + 1) / 2²ⁿ⁺⁴;
• geom4(r) + geom4(r) + 1 returns n with probability 36n / 2²ⁿ⁺⁴;
• combine these with probabilities [4/9, 5/9] to yield (6n + 1) / 2²ⁿ⁺², as required.

Implement (2) as the outcome of 3 coin tossing experiments of 2n tosses with success defined as equal numbers of heads and tails in each trial.

This class illustrates how to return an exact result using coin tosses only. A more efficient implementation (which is still exact) would replace prob59 by r.Prob(5,9) and geom4 by LeadingZeros z; z(r)/2.

Definition at line 104 of file InversePiProb.hpp.

Member Function Documentation

template<class Random >

bool RandomLib::InversePiProb::operator() ( Random &  r )

inline

Return true with probability 1/π.

Template Parameters

Random

the type of the random generator.

Parameters

[in,out]

r

a random generator.

Returns

true with probability 1/π.

Definition at line 140 of file InversePiProb.hpp.

---

The documentation for this class was generated from the following file:

• InversePiProb.hpp