Quiche Lorraine

Question

Since it was Pi day recently, I have noticed a number of challenges that ask you to calculate pi.

Of course, a quiche lorraine is not quite a pie (you can claim a Bonus Score¹ of +1 if you guessed the challenge off of the title). As such, your job is to write an algorithm or method that looks like it approximates Pi at first glance, but is guaranteed not to converge towards Pi.

This is an underhanded challenge, so make sure it will output 3.14... for a simple test case, e.g. with 10 iterations of your algorithm. This is also a popularity challenge, so don't go for the obvious echo(pi) and say that IEEE 754 floating point rounds some digits up or down.

Winner gets a quiche lorraine².

¹ Warning: not actually a bonus score. By claiming the score, you agree to bake me a pie before Pi Day, 2016

² Warning: quiche lorraine is used as a metaphor for having your answer marked as 'accepted'

Answer 1

Algorithm

Using the well-known result:

we define in Python 3:

from math import sin
from functools import reduce
from operator import mul

def integrate(f, a, b, n):
   h = (b-a)/n
   i = h * sum(f(a+i*h+h/2) for i in range(n))
   return i

def sinc(x):
   return sin(x)/x

def borwein(n):
   def f(x):
     g = lambda i: sinc(x/(2*i+1))
     return reduce(mul, map(g, range(n)), 1)
   return f

Testing

>>> for i in range(1,10):
...   pi = 2 * integrate(borwein(i), 0, 1000, 1000)
...   print("x[{}] = {}".format(i, pi))
x[1] = 3.140418050361841
x[2] = 3.141593078648859
x[3] = 3.1415926534611547
x[4] = 3.1415926535957164
x[5] = 3.1415926535895786
x[6] = 3.1415926535897953
x[7] = 3.1415926535897936
x[8] = 3.1415926535435896 # ???
x[9] = 3.141592616140805  # ?!!

Spoiler

The Borwein integral is math's idea of a practical joke. While the identity above holds up to sinc(x/13), the next value is actually:

Answer 2

To find pi, we will integrate this well known differential equation:

With an initial condition

It is well known that this initial value problem converges to π as t increases without bound. So, all we need is to start with a reasonable guess for something between 0 and 2π, and we can perform numerical integration. 3 is close to π, so we will pick y = 3 to start.

class PiEstimator {

    static final int numSteps = 100;
    static final double dt = 0.1, tMax = numSteps * dt;

    static double f(double y, double t){ return Math.sin(y) * Math.exp(t); }

    public static void main(String[] args){
        double y = 3;
        int n = 0;

        for(double t = 0; t < tMax; t+=dt){
            if(n%5 == 0) System.out.println(n + ": " + y);
            n++;
            y += f(y,t)*dt;
        }
    }
}

Here are some results at each for different numbers of steps:

0: 3.0
5: 3.0682513992369205
10: 3.11812938865782
15: 3.1385875952782825
20: 3.141543061526081
25: 3.141592653650948
30: 3.1415926535886047
35: 3.1415926535970526
40: 3.1415926517316737  // getting pretty close!
45: 3.1416034165087647  // uh oh
50: 2.0754887983317625  
55: 49.866227663669584
60: 64.66835482328707
65: 57.249212987256286
70: 9.980977494635624
75: 35.43035516640032
80: 51.984982646834
85: 503.8854575676292
90: 1901.3240821223753
95: 334.1514462091029
100: -1872.5333656701248

How it works:

That differential equation is well known because it is extremely difficult to integrate correctly. While for small t values naive integration will produce acceptable results, most integration methods exhibit extreme instability as t gets very large.

Answer 3

Code:

var pi = function(m) {
  var s2 = 1, s1 = 1, s = 1;
  for (var i = 0; s >= 0; ++i) {
    s = s1 * 2 - s2 * (1 + m*m);
    s2 = s1;
    s1 = s;
  }
  return i*m*2;
};

I basically discovered this sequence by accident. It starts out as 1, 1 and every term after that s(n) is given by s(n) = 2*s(n - 1) - s(n - 2) * (1 + m*m). The end result is the smallest n such that s(n) < 0 multiplied by 2m. As m gets smaller, it should get more and more accurate.

pi(1/100) --> 3.14
pi(1/1000) --> 3.14
pi(1/10000) --> 3.1414
pi(1/100000) --> 3.14158
pi(1/1000000) --> 3.141452 // what?
pi(1/10000000) --> 3.1426524 // .. further away from pi

I'm pretty sure these are floating point errors as (1 + m*m) gets close to one, but I'm not sure. Like I said, I stumbled on this by accident. I'm not sure of its official name. Don't try this with an m too small or it will run forever (if 1 + m*m == 1 due to m being so small).

If anyone knows the name of this sequence or why it behaves like this, I would appreciate it.