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.

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!

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.

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.