[Question]
  [
**Closed.** This question is [off-topic](/help/closed-questions). It is not currently accepting answers.
 
 
---
 This question does not appear to be about code golf or coding challenges within the scope defined in the [help center](https://codegolf.stackexchange.com/help/on-topic).
Closed 7 years ago.
**Locked**. This question and its answers are [locked](/help/locked-posts) because the question is off-topic but has historical significance. It is not currently accepting new answers or interactions.
 
 
Since it was Pi day recently, I have [noticed](https://codegolf.stackexchange.com/questions/47759/monte-carlo-estimator-of-pi) a [number](https://codegolf.stackexchange.com/questions/47808/transmit-pi-precisely) of [challenges](https://codegolf.stackexchange.com/questions/47912/calculate-%CF%80-with-quadratic-convergence) 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]
## Algorithm
Using the well-known result:
![enter image description here](https://i.stack.imgur.com/ERe8P.png)
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](http://en.wikipedia.org/wiki/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:  
>   
> 
> ![enter image description here](https://i.stack.imgur.com/pCTcT.png)
> 
> 
> 
[Answer]
To find pi, we will integrate this well known differential equation:
![> dy/dt = sin(y)*exp(t)](https://i.stack.imgur.com/UQVsV.gif)
With an initial condition
![> 0 < y0 < 2*pi](https://i.stack.imgur.com/zXjJ4.gif)
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]
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.
]