Bresenham’s Pie

Approximating \pi is always fun. Some approximations are well known, like Zu Chongzhi’s, \frac{355}{113}, 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 \pi. 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 O(r^2) process because a circle of (integer) radius r is contained in a 4r^2 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 r^2 area in r steps! Bresenham’s is actually a clever way of walking around the circle using an error-correction scheme. As we walk the circle along (x,y) points, incrementing x as we go, we sum the ys, because they represent the height, and therefore the number of points in the circle between (x,y) and the axis at (x,0).

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

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 O(r), and that its error is proportional to \frac{1}{r}. Indeed, above we see that r=10^9 gives an error of \approx 8\times{}10^{-10}, which is not too bad.

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s

%d bloggers like this: