A Different Kind of Pi

In my last post I left my readers (Hi Mom!) with a teaser about a possible improvement to my pi-calculating code. Let’s see how it panned out.

My naïve arc-drawing code used two square-root calculations per point – pretty slow. I figured that there must be a quicker version out there, and sure enough, the gold standard in translating curves to pixels is something called Bresenham’s Algorithm. (Edit: the Wikipedia article is awful. Here’s a better one.) Here’s my code, cribbed from this page:

function result = bresarc(r)
% BRESARC(R) - draw a quarter-arc of a circle using Bresenham's algorithm

r = abs(round(r));
x = -r;
y = 0;
err = 2 - 2 * r;

idx = 1;
result = zeros(round(sqrt(2) * r ), 2);

while x <= 0
   result(idx, [1 2]) = [-x y];
   idx = idx + 1;
   t = err;
   if t > x
       x = x + 1;
       err = err + 2 * x + 1;
   end
   if t <= y
       y = y + 1;
       err = err + 2 * y + 1;
   end
end
end

It’s about four times faster than my arcwalk function, and produces a better approximation to the circle, with fewer points.

Bresenham arc, compared to arc and simplified Bresenham arc

Bresenham arc, compared to arc and simplified Bresenham arc

Here’s an oddity. Look at line 10, where I pre-calculate the number of points to return in the result. It turns out (I discovered by experimentation) that this algorithm always generates radius*\sqrt2 points for some reason. Something this obvious, in an algorithm that’s very well known, should probably already be common knowledge – but I couldn’t find any evidence of it by Googling for a few minutes. Perhaps my 15 minutes of fame are on the way!

Let’s see how it does in my pi-calculating algorithm. The original form, slightly updated from last time (better memory management):

clc
clear

maxpow = 5;
epsilon = sqrt(2);

tic

for r = 10 .^ (1:maxpow) 
    pts = arcwalk(r);                       % make the arc
    pts = dpsimp(pts, epsilon);             % simplify it
    pts = pts(2:end,:) - pts(1:end-1,:);    % vector subtraction
    s = sum(hypot(pts(:,1), pts(:,2)));     % add up the lengths
    fprintf('For radius %d, pi is about %f\n', r, 2*s/r)
end

toc

The output:

For radius 10, pi is about 3.164823
For radius 100, pi is about 3.152301
For radius 1000, pi is about 3.142237
For radius 10000, pi is about 3.141706
For radius 100000, pi is about 3.141607
Elapsed time is 41.948334 seconds.

Swapping out arcwalk for bresarc, we get:

For radius 10, pi is about 3.170979
For radius 100, pi is about 3.143267
For radius 1000, pi is about 3.141778
For radius 10000, pi is about 3.141648
For radius 100000, pi is about 3.141602
Elapsed time is 29.530892 seconds.

Not bad! Let’s see about speeding it up some more.

The killer that remains is the DP polyline simplification routine dpsimp. Hmmm…. the new arc algorithm is closer to the circle. How would it work if I just calculated the arc length of the output of the bresarc function directly, without attempting to simplify it? It will be longer than the true arc (I think; proof is left as an exercise for the reader), but… how much more? Let’s find out.

After commenting out line 11 above, we get:

For radius 10, pi is about 3.297056
For radius 100, pi is about 3.308772
For radius 1000, pi is about 3.313458
For radius 10000, pi is about 3.313693
For radius 100000, pi is about 3.313704
Elapsed time is 0.283709 seconds.

Holy crap! 100 times faster! Farewell, favorite recursive algorithm.

Hmmm. Is that value converging? Let’s see. Bumping up maxpower from 5 to 7 and prettying up the output:

For radius       10, pi is about 3.297056274848
For radius      100, pi is about 3.308772003600
For radius     1000, pi is about 3.313458295101
For radius    10000, pi is about 3.313692609677
For radius   100000, pi is about 3.313704325403
For radius  1000000, pi is about 3.313709011707
For radius 10000000, pi is about 3.313708543080
Elapsed time is 29.052004 seconds.

It kinda looks like it. So what have we discovered here? A new constant? It would be the ratio of the perimeter length of a Bresenham circle to its diameter. Call it pi-sub-b, or \pi_b. For now, let’s assume that \pi_b\approx3.31371. We know that this number is transcendental like pi, because pi_b is made up of line segments with two possible lengths: 1 and sqrt(2) (which is also transcendental).

An obvious hack: we can use the (presumed) fact that pi_b is a constant to calculate pi. We figure out the difference between pi and pi_b, then subtract it from the calculated value of pi_b:

EDU>> pib = 3.31371

pib =

                   3.31371

EDU>> pidiff = pib - pi

pidiff =

         0.172117346410207

EDU>> 

Let’s see how it works:

clc
clear

maxpow = 7;
pidiff = 0.172117346410207;

tic

for r = 10 .^ (1:maxpow) 
    pts = bresarc(r);                       % make the arc
    pts = pts(2:end,:) - pts(1:end-1,:);    % vector subtraction
    s = sum(hypot(pts(:,1), pts(:,2)));     % add up the lengths
    fprintf('For radius %8d, pi is about %.12f\n', r, 2*s/r - pidiff)
end

toc
For radius       10, pi is about 3.124938928438
For radius      100, pi is about 3.136654657190
For radius     1000, pi is about 3.141340948691
For radius    10000, pi is about 3.141575263266
For radius   100000, pi is about 3.141586978993
For radius  1000000, pi is about 3.141591665297
For radius 10000000, pi is about 3.141591196670
Elapsed time is 29.198759 seconds.

Of course, that’s absolutely cheating. I baked the value for pi into the program, as part of the pidiff fudge factor.

But the possibility exists that there’s a more direct route to this result. Pi, the square root of two, and pi_b may all be related in a deep way. Almost certainly the Bresenham algorithm has \sqrt2 baked into it, as seen at the beginning of this post.

Consider a Bresenham circle. It has a perimeter length. The length won’t change if we rearrange the line segments it’s made of. It is made of horizontal lines of length 1, vertical lines of length 1, and diagonal lines of length \sqrt2. It must be possible to rearrange these into a square with the corners cut off.

Ah, the ancient dream of squaring the circle! My 15 minutes of fame approaches!

About Doug

I grew up an Air Force brat and have visited every European country except those once behind the Iron Curtain (they wouldn't let my father in for some reason). Now I'm enrolled in the Aerospace Engineering program at NMSU in Las Cruces, NM.
This entry was posted in . Bookmark the permalink.

2 Responses to A Different Kind of Pi

  1. Steve VanDevender says:

    There’s a basic geometric reason that the midpoint circle algorithm plots r*sqrt(2) points in a quarter-circle. First notice that the Bresenham line algorithm plotting a line between (x1, y1) to (x2, y2) will plot max(abs(x2-x1), abs(y2-y1)) points, no matter what the slope of the line; that is, it plots a point at each integral value of x or y along the line.

    The midpoint circle algorithm (as pointed out, it’s not due to Bresenham, although it’s based on the ideas in the Bresenham line algorithm) is normally used to generate the points along one-eighth of the circumference of a circle with the rest of the points obtaimed by symmetry. So it may typically generate points at integral values of x and y from the point (x=r, y=0) to a point pi/4 (or 45 degrees) along the circle from there, which is the point (x=r*cos(pi/4), y=r*sin(pi/4)), or in numeric terms (x=r*sqrt(2)/2, y=r*sqrt(2)/2). So one-eighth of a circle requires plotting r*sqrt(2)/2 points, and by extension one-quarter of a circle requires plotting r*sqrt(2) points, or the entire circle 4*r*sqrt(2) points.

  2. Doug says:

    That makes sense to me, thanks.

    I convinced myself of it after staring at the “Bresenham circle” and figuring out it could be rearranged into an “integer-ized” octagon and then doing some geometry, which I always enjoy.

Leave a Reply

Your email address will not be published. Required fields are marked *

You may use these HTML tags and attributes: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <strike> <strong>