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

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 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 . For now, let’s assume that . 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 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 . 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!