## 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 $y$s, 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.