Approximating is always fun. Some approximations are well known, like Zu Chongzhi’s, , that is fairly precise (it gives 6 digits).

We can derive a good approximation using continue fractions, series, or some spigot algorithm. We can also use Monte Carlo methods, such as drawing uniformly distributed values inside a square, count those who falls in the circle, then use the ratio of inside points to outside points to evaluate . This converges slowly. How do we evaluate the area of the circle then? Well, we could do it exhaustively, but in a smart way.

If we count the number of (discrete) points inside a circle, one by one, that would be an process because a circle of (integer) radius is contained in a points square. But we can actually only walk on the edge of the disc, that is, on the circle itself in a much faster way.

Using Bresenham’s midpoint circle algorithm, and symmetries, we can cover an area in steps! Bresenham’s is actually a clever way of walking around the circle using an error-correction scheme. As we walk the circle along points, incrementing as we go, we sum the s, because they represent the height, and therefore the number of points in the circle between and the axis at .

We have to be a bit careful on the symmetries, as the algorithm works only on one eighth of the circle. If is on the circle, so will ; and we have to be careful to add only when changes. We compute the area of a quarter circle, times four (minus overlaps), and we get a fairly accurate estimation of .

A quick implementation would look like this:

// c++20
#include
#include
#include
template T sqr(const T & x) { return x*x; }
int main()
{
int64_t r=1'000'000'000;
int64_t x=0;
int64_t y=r;
int64_t oy=0; // old y
int64_t d=5-4*r; // everything times 4 to avoid fractions
size_t s1=0;
size_t s2=0;
while (x<y)
{
s1+=y;
if (d<0)
d+=4*x+6;
else
{
d+=8*(x-y)+20;
y--;
}
x++;
if (oy!=y) // symmetry
s2+=x;
oy=y;
}
double pi_est=4.0*(s1+s2-r/2)/sqr(r);
std::cout
<< (s1+s2) << '/' << sqr(r) << std::endl
<< pi_est << std::endl
<< pi_est-std::numbers::pi << std::endl;
return 0;
}

The output is:

time a.out
785398163690333774/1000000000000000000
3.14159
-8.28458e-10
real 0m0.412s
user 0m0.412s
sys 0m0.000s

The interesting thing about this algorithm is that it is , and that its error is proportional to . Indeed, above we see that gives an error of , which is not too bad.

This entry was posted on Tuesday, March 29th, 2022 at 11:35 am and is filed under Uncategorized. You can follow any responses to this entry through the RSS 2.0 feed.
You can leave a response, or trackback from your own site.

RT @FondsKheops: 📚 [Editions Khéops/ Editions du Louvre] "La voix des hiéroglyphes – Nouvelle édition revue et augmentée 2022"
par Christop… 33 minutes ago

RT @BLMedieval: Our Sherborne Missal episode of Moving Pictures radio programme is going to be repeated on Radio 4 today, 12 August at 11am… 33 minutes ago