# Quiche Lorraine [closed]

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'

• Related: link – Sp3000 Mar 19 '15 at 13:06
• I'm voting to close this question as off-topic because underhanded challenges are no-longer on-topic here. meta.codegolf.stackexchange.com/a/8326/20469 – cat Apr 15 '16 at 14:20

## 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:

• Probably one of the best answer to a underhanded question in recent times. – Optimizer Mar 19 '15 at 15:12
• "math's idea of a practical joke". +1 – unclemeat Mar 19 '15 at 22:31
• That's a good one! IIRC, one of the more well-known jokes with this integral was when someone recorded the results up to the weird one on Wolfram Alpha, and sent a bug report... Which the WA devs spent ages trying to fix =) – Mints97 Mar 20 '15 at 6:30
• This reference gives a good explanation of what is going on. – TonioElGringo Mar 20 '15 at 10:07

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.

• @UriZarfaty This wikipedia article explains it pretty well: en.wikipedia.org/wiki/Stiff_equation – AJMansfield Mar 19 '15 at 15:35
• What is n?... – Cole Johnson Mar 20 '15 at 23:30
• @AJMansfield I meant: it's not declared anywhere. Your for deceleration uses t, but your loop uses n. – Cole Johnson Mar 21 '15 at 4:59
• @ColeJohnson I just fixed it. – AJMansfield Mar 21 '15 at 10:40
• I think your differential equation should read dy/dt = sin(y)*exp(t). – David Zhang Mar 21 '15 at 14:24

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.

• I think this is due to cancellation, which is a loss of digits when subtracting two almost equal numbers. S1 and s2 are almost equal after an iteration. – Sanchises Oct 4 '15 at 20:14
• I'm yet to find out how it works, but it reminds me of something I once did: I repeatedly took the cumulative sum of a noisy signal and normalized it to mean 0, max value 1. This would converge to a sine wave, since that's the only signal that is its own anti-derivative (with a phase offset). – Sanchises Oct 4 '15 at 21:26
• I asked it over at math.SE, and got this answer. – Sanchises Jan 13 '16 at 15:04