Quiche Lorraine [cerrado]

52

Como fue el día de Pi recientemente, he notado una serie de desafíos que le piden que calcule pi.

Por supuesto, una quiche lorraine no es del todo un pastel (puede reclamar un Bonus Score¹ de +1 si adivinó el desafío fuera del título). Como tal, su trabajo es escribir un algoritmo o método que parezca que se aproxima a Pi a primera vista, pero se garantiza que no converja hacia Pi.

Este es un desafío poco claro, así que asegúrese de que salga 3.14 ... para un caso de prueba simple, por ejemplo, con 10 iteraciones de su algoritmo. Este también es un desafío de popularidad, así que no vaya por lo obvio echo(pi)y diga que el punto flotante IEEE 754 redondea algunos dígitos hacia arriba o hacia abajo.

El ganador recibe un quiche lorraine².

¹ Advertencia: en realidad no es un puntaje de bonificación. Al reclamar el puntaje, acepta hornearme un pastel antes del Día de Pi, 2016

² Advertencia: la quiche lorraine se usa como una metáfora para marcar su respuesta como 'aceptada'

Sanchises
fuente
Relacionado: enlace
Sp3000
2
Estoy votando para cerrar esta pregunta como fuera de tema porque los desafíos poco claros ya no están en el tema aquí. meta.codegolf.stackexchange.com/a/8326/20469
gato

Respuestas:

77

Algoritmo

Usando el resultado bien conocido:

ingrese la descripción de la imagen aquí

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

Pruebas

>>> 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  # ?!!

Revelación

La integral de Borwein es la idea matemática de una broma práctica. Si bien la identidad anterior se sostiene hasta sinc (x / 13), el siguiente valor es en realidad:

ingrese la descripción de la imagen aquí

Uri Granta
fuente
12
Probablemente una de las mejores respuestas a una pregunta encubierta en los últimos tiempos.
Optimizador
14
"idea matemática de una broma práctica". +1
unclemeat
16
¡Esa es buena! IIRC, una de las bromas más conocidas con esta integral fue cuando alguien grabó los resultados hasta el extraño en Wolfram Alpha, y envió un informe de error ... Que los desarrolladores de WA pasaron años tratando de arreglar =)
Mints97
3
Esta referencia da una buena explicación de lo que está sucediendo.
TonioElGringo
59

Para encontrar pi, integraremos esta conocida ecuación diferencial:

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

Con una condición inicial

> 0 <y0 <2 * pi

Es bien sabido que este problema de valor inicial converge a π a medida que t aumenta sin límite. Entonces, todo lo que necesitamos es comenzar con una conjetura razonable para algo entre 0 y 2π, y podemos realizar una integración numérica. 3 está cerca de π, por lo que elegiremos y = 3 para comenzar.

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

Aquí hay algunos resultados en cada uno para diferentes números de pasos:

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

Cómo funciona:

Esa ecuación diferencial es bien conocida porque es extremadamente difícil de integrar correctamente. Mientras que para valores t pequeños la integración ingenua producirá resultados aceptables, la mayoría de los métodos de integración exhiben una inestabilidad extrema a medida que t se vuelve muy grande.

AJMansfield
fuente
44
@UriZarfaty Este artículo de Wikipedia lo explica bastante bien: en.wikipedia.org/wiki/Stiff_equation
AJMansfield
1
¿Qué es n? ...
Cole Johnson
1
@AJMansfield quise decir: no está declarado en ninguna parte. Tu fordesaceleración usa t, pero tu ciclo usa n.
Cole Johnson
1
@ColeJohnson Lo acabo de arreglar.
AJMansfield
2
Creo que su ecuación diferencial debería leer dy / dt = sin (y) * exp (t).
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;
};

Básicamente descubrí esta secuencia por accidente. Comienza como 1, 1y cada término después de eso s(n)viene dado por s(n) = 2*s(n - 1) - s(n - 2) * (1 + m*m). El resultado final es el más pequeño n, s(n) < 0multiplicado por 2m. A medida que se mhace más pequeño, debería ser cada vez más 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

Estoy bastante seguro de que se trata de errores de coma flotante, ya que (1 + m*m)se acerca a uno, pero no estoy seguro. Como dije, me topé con esto por accidente. No estoy seguro de su nombre oficial. No intente esto con un mtamaño demasiado pequeño o se ejecutará para siempre (si es 1 + m*m == 1debido a que mes muy pequeño).

Si alguien sabe el nombre de esta secuencia o por qué se comporta así, lo agradecería.

soktinpk
fuente
Creo que esto se debe a la cancelación, que es una pérdida de dígitos al restar dos números casi iguales. S1 y s2 son casi iguales después de una iteración.
Sanchises
1
Todavía tengo que descubrir cómo funciona, pero me recuerda algo que hice una vez: repetidamente tomé la suma acumulativa de una señal ruidosa y la normalicé para que significara 0, valor máximo 1. Esto convergería en una onda sinusoidal, ya que esa es la única señal que es su propia anti-derivada (con un desplazamiento de fase).
Sanchises
Lo pregunté en matemáticas. SE, y obtuve esta respuesta.
Sanchises