Quiche Lorraine [fechado]

52

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'

Sanchises
fonte
Relacionado: link
Sp3000
2
Estou votando para encerrar esta questão como fora do tópico, porque os desafios secretos não estão mais no tópico aqui. meta.codegolf.stackexchange.com/a/8326/20469
cat

Respostas:

77

Algoritmo

Usando o resultado conhecido:

insira a descrição da imagem aqui

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:

insira a descrição da imagem aqui

Uri Granta
fonte
12
Provavelmente, uma das melhores respostas para uma pergunta secreta nos últimos tempos.
Optimizer
14
"idéia matemática de uma piada prática". +1
unclemeat 19/03/2015
16
Essa é boa! IIRC, uma das piadas mais conhecidos com este integrante foi quando alguém gravou os resultados até aquele estranho no Wolfram Alpha, e enviou um relatório de erro ... Que o WA devs passado as idades tentando consertar =)
Mints97
3
Esta referência fornece uma boa explicação do que está acontecendo.
TonioElGringo
59

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

> dy / dt = sin (y) * exp (t)

Com uma condição inicial

> 0 <y0 <2 * pi

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

AJMansfield
fonte
4
@UriZarfaty Este artigo da Wikipedia explica muito bem: en.wikipedia.org/wiki/Stiff_equation
AJMansfield 19/15/15
11
O que é n...?
Cole Johnson
11
@AJMansfield eu quis dizer: não é declarado em lugar nenhum. Sua fordesaceleração usa t, mas seu loop usa n.
Cole Johnson
11
@ColeJohnson Acabei de corrigir.
precisa saber é o seguinte
2
Eu acho que sua equação diferencial deve ser dy / dt = sin (y) * exp (t).
21815 David Zhang
6

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.

soktinpk
fonte
Eu acho que isso se deve ao cancelamento, que é uma perda de dígitos ao subtrair dois números quase iguais. S1 e s2 são quase iguais após uma iteração.
Sanchises
11
Ainda estou para descobrir como funciona, mas isso me lembra algo que fiz uma vez: peguei repetidamente a soma cumulativa de um sinal barulhento e normalizei para significar 0, valor máximo 1. Isso convergiria para uma onda senoidal, já que esse é o único sinal que é seu próprio anti-derivado (com um deslocamento de fase).
Sanchises
Eu perguntei no math.SE, e obtive esta resposta.
Sanchises