Quiche Lorraine [fechado]

Desde que foi o dia do Pi recentemente, notei uma série de desafios que solicitam que você calcule o pi.

Obviamente, uma quiche lorraine não é exatamente uma torta (você pode reivindicar uma pontuação de bônus¹ de +1 se adivinhou o desafio do título). Como tal, seu trabalho é escrever um algoritmo ou método que pareça com o Pi à primeira vista, mas é garantido que não convergirá para o Pi.

Este é um desafio secreto, portanto, verifique se ele produzirá 3,14 ... para um caso de teste simples, por exemplo, com 10 iterações do seu algoritmo. Esse também é um desafio de popularidade; portanto, não vá pelo óbvio echo(pi)e diga que o ponto flutuante IEEE 754 arredonda alguns dígitos para cima ou para baixo.

O vencedor ganha uma quiche lorraine².

¹ Aviso: na verdade não é uma pontuação de bônus. Ao reivindicar a pontuação, você concorda em me fazer uma torta antes do Dia do Pi, 2016

² Atenção: quiche lorraine é usada como uma metáfora para ter sua resposta marcada como 'aceita'

Algoritmo

Usando o resultado conhecido:

nós definimos em 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

Testando

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

A integral de Borwein é a idéia matemática de uma piada prática. Embora a identidade acima seja sinc (x / 13), o próximo valor é realmente:

Para encontrar pi, integraremos esta conhecida equação diferencial:

Com uma condição inicial

É sabido que esse problema de valor inicial converge para π à medida que t aumenta sem limite. Portanto, tudo o que precisamos é começar com um palpite razoável para algo entre 0 e 2π, e podemos realizar a integração numérica. 3 está próximo de π, então escolheremos y = 3 para começar.

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;
        }
    }
}

Aqui estão alguns resultados em cada um para diferentes números de etapas:

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

Como funciona:

Essa equação diferencial é bem conhecida porque é extremamente difícil de integrar corretamente. Enquanto para pequenos valores de t a integração ingênua produzirá resultados aceitáveis, a maioria dos métodos de integração exibe extrema instabilidade à medida que t se torna muito grande.

Código:

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;
};

Eu basicamente descobri essa sequência por acidente. Começa como 1, 1todos os termos após o que s(n)é dado por s(n) = 2*s(n - 1) - s(n - 2) * (1 + m*m). O resultado final é o menor n, s(n) < 0multiplicado por 2m. À medida que mdiminui, deve ficar cada vez mais preciso.

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

Tenho certeza de que esses erros de ponto flutuante (1 + m*m)se aproximam de um, mas não tenho certeza. Como eu disse, eu tropecei nisso por acidente. Não tenho certeza do nome oficial. Não tente fazer isso com uma mmuito pequena ou ele vai correr para sempre (se 1 + m*m == 1devido a mser tão pequeno).

Se alguém souber o nome dessa sequência ou por que ela se comporta assim, eu agradeceria.