Bresenham’s Pie

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.