The Limber Lambda

This one got me thinking for a long time, mostly because I subconsciously eliminated “cheat” solutions.  The problem is to calculate an estimate of π by working out the perimeter of an n-sided circle-inscribing polygon.  Let’s start with some pictures:

Now, the obvious solution (focus on Figure A) is to use a trigonmetric function to solve the problem.  Using a convenient r = 0.5 (implies that the perimeter of the circle is π) one side of our polygon (d) can be determined as follows:

d = 2rsin θ

And since the perimeter of the polygon = p = nd, and r = 0.5, and θ = 2π/2n, then:

p = nsin (π/n)

All too easy, and well … we’re using π to find an approximation of π.  Somehow that doesn’t sit well with me.  Which is the reason why I took two days to come up with a solution that didn’t involve trigonometric functions, or π.  Enter Figure B.

We apply a method of virtual construction–let’s assume a cartesian plane, our circle having centre (0,0), and radius r = 0.5.  The method is iterative, so we’ll need to start with a guess of the length of a side of the polygon, R_guess =P_guess / n, where P_guess is our “perimeter guess”.  P_guess can’t be any larger than π, so let’s make the upper bound = 3.2.  R_guess then forms the radius of a “construction circle”, with centre (x,y), starting at (x₀,y₀) in the diagram.  The intersection of said construction circle and our original circle will determine a polygon corner; of course, since there are two intersection points, we eliminate the one that has already been “visited” (or occurs where y is negative in the case of the first attempt).  We progressively “construct” our way around the circle n times, each time using the previously determined “corner” as the centre of the next construction circle, finally arriving at the final intersection point: (x_n, y_n).  If R_guess was correct, then (x_n, y_n) would equal (x₀,y₀).  Adjusting R_guess for the next iteration then becomes a case of increasing P_guess if y_n turned out negative and decreasing P_guess if y_n turned out positive.  A simple zero-by-halving (aka Bisection) technique is employed to get a value for the perimeter that satisfies the accuracy constraint (1e-9).

The source code for my π-less solution is below:

using System;

class Archimedes
{
    public double approximatePi(int numSides)
    {
        double upperPi = 3.5;
        double lowerPi = 0.0;
        while (Math.Abs(upperPi - lowerPi) > 1e-10)
        {
            bool first = true;
            double cx = 0.5;
            double cy = 0.0;
            double piGuess = (upperPi + lowerPi) / 2;
            double lastCx = 0.0, lastCy = 0.0;
            for (int i = 0; i < numSides; i++)
            {
                double[][] intersections = GetCirclesIntersection(cx, cy, piGuess, numSides);
                if (intersections[0][0] == Double.NaN ||
                    intersections[0][1] == Double.NaN ||
                    intersections[1][0] == Double.NaN ||
                    intersections[1][1] == Double.NaN)
                {
                    cy = 1;
                    break;
                }
                if (first)
                {
                    if (intersections[0][1] > 0)
                    {
                        lastCx = cx;
                        lastCy = cy;
                        cx = intersections[0][0];
                        cy = intersections[0][1];
                    }
                    else
                    {
                        lastCx = cx;
                        lastCy = cy;
                        cx = intersections[1][0];
                        cy = intersections[1][1];
                    }
                    first = false;
                }
                else
                {
                    if (Math.Abs(intersections[0][0] - lastCx) < 1e-10 &&
                        Math.Abs(intersections[0][1] - lastCy) < 1e-10)
                    {
                        lastCx = cx;
                        lastCy = cy;
                        cx = intersections[1][0];
                        cy = intersections[1][1];
                    }
                    else
                    {
                        lastCx = cx;
                        lastCy = cy;
                        cx = intersections[0][0];
                        cy = intersections[0][1];
                    }
                }
            }
            if (cy < 0.0)
            {
                lowerPi = piGuess;
            }
            else
            {
                upperPi = piGuess;
            }
        }
        return upperPi;
    }

    private double[][] GetCirclesIntersection(double xo, double yo, double piGuess, int n)
    {
        double e = xo;
        double f = yo;
        double p = Math.Sqrt(e*e + f*f);
        double k = (p * p + 0.25 - (piGuess / n) * (piGuess / n)) / (2 * p);
        return new double[2][] {
            new double[] {
                (e*k)/p + (f/p) * Math.Sqrt(0.25-k*k),
                (f*k)/p - (e/p) * Math.Sqrt(0.25-k*k) },
            new double[] {
                (e*k)/p - (f/p) * Math.Sqrt(0.25-k*k),
                (f*k)/p + (e/p) * Math.Sqrt(0.25-k*k)
            }
        };
    }
}