π like an Egyptian

In the Rhind Mathematical Papyrus (RMP) we find that the Egyptians used the approximation

\displaystyle 4\times\left(\frac{8}{9}\right)^2 =4\times\frac{64}{81} =\frac{256}{81} \approx{3.16}

For \pi. The RMP shows how to arrive at this result: you construct a 9\times{9} grid and draw an octagon on it that approximates the circle. Turns out this irregular octagon as \frac{7}{9} of the area of the bounding square. But the ancient Egyptians rounded that value, for reasons unknown, from 63 to 64 (I have a theory on why; but that may be a later post). This gives the above approximation.

egyptian-pie-crop-small

But what if we increase the grid density? 9\times{9} is convenient for the computationally-challenge ancient Egyptians, but can we do better? In previous posts, we have examined how to approximate \pi by rejection and by approximating the area of a circle on a discrete grid. The rejection method showed very slow convergence and the scanning method was computationally expensive.

Let’s revisit the scanning method. In the previous post, I used floating points numbers, which, admittedly, are in themselves a source of numerical error. Let’s rewrite it using integer-only arithmetic.

template <typename T> T sqr(const T & x) { return x*x; }

uint64_t pi_like_an_egyptian(uint64_t steps)
 {
  const uint64_t radius=sqr(steps-0.5);

  uint64_t inside=0;

  for (uint64_t x=0;x<steps;x++)
   {
    const uint64_t sx=sqr(x);
    for (uint64_t y=0;y<steps;y++)
     inside+=(sx+sqr(y)<=radius);
   }

  return inside;
 }

This method is not any better, complexity-wise, than the previous implementation using floating point numbers. It scans every point on the grid, checks if the midpoint of the square falls inside the circle or not, and increments a counter accordingly. It is a O(n^2) algorithm for a grid density of n\times{n}—unsurprisingly.

What now if we just try to find where the circle’s edge is on the grid? If we somehow only followed the circle, we would not have to scan the entire grid. Well, we could use the Pythagorean theorem and find, from a known point, which neighboring points lie on the edge of the circle. That sounds complicated. Fortunately, some good fellow named Bresenham devised an algorithm to do just that!

The midpoint circle algorithm uses the derivatives of the circle equation to find, without extracting square roots, the next point to draw. But this next point to draw lies on the edge of the circle, and that’s what we want. A possible implementation of this algorithm on our n\times{n} grid is as follows:

uint64_t bresenham_pi(uint64_t radius)
 {
  int64_t x=0;
  int64_t y=radius;
  int64_t d=3-2*radius;
  uint64_t inside=0;

  while (x&lt;y)
   {
    x++;
    inside+=y;

    if (x==y) break;

    if (d&lt;0)
     d+=4*x+6;
    else
     {
      d+=4*(x-y)+10;
      y--;
     }
   }

  return 2*inside-sqr(y);
 }

It only goes one eighth of a circle then uses inclusion-exclusion principle to compute the correct number of squares inside the circle.

*
* *

Let’s compare the two algorithms. The following figure shows the (log-scale) absolute errors of approximation given the grid size. The green is the integer scanning algorithm, the red, the algorithm using Bresenham’s midpoint circle algorithm.

egyptian-h1

Bresenham’s method (let’s call it that for short) does much worse. Well, that might not be surprising as it uses an approximation to compute the circle. Visually, it will look nice, but in fact, sometimes the limit is drawn a bit too much inside the circle, sometimes, and more often, just a bit outside. In need, we have to a 68000\times{}68000 grid to reach the precision the brute-force scanning algorithm gets with a 50000\times{50000} grid, around 4.3\times{10}^{-6}:

egyptian-h2

But while the first algorithm has a O(n^2) complexity, Bresenham method’s has O(n) complexity. It is much faster. It will allow us to go into million-sized grids in an insignificant fraction of a second.

egyptian-h3

Pushing the grid to 2000000\times{2000000}, we get an error of 1.4\times{}10^{-7}, which isn’t all that bad!

*
* *

Bresenham’s circle algorithm is finicky and uses constants that depends on the wanted precision. The constants I used seem to be giving a bit more precise circle than Wikipedia’s; but that more or less just numerical tinkering. Maybe a Pythagorean-based method that corrects the path and allows jumps further than one pixel would do better. Questions for later.

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 )

Google+ photo

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

Connecting to %s

%d bloggers like this: