Quiche Lorraine [fermé]

52

Depuis que c'était la journée du pi récemment, j'ai remarqué un certain nombre de défis qui vous demandent de calculer le pi.

Bien sûr, une quiche lorraine n’est pas une tarte (vous pouvez prétendre à un Bonus Bonus¹ de +1 si vous avez deviné le défi du titre). En tant que tel, votre travail consiste à écrire un algorithme ou une méthode qui ressemble au premier abord à Pi, mais dont la garantie est de ne pas converger vers Pi.

C'est un défi sournois, alors assurez-vous qu'il générera 3.14 ... pour un cas de test simple, par exemple avec 10 itérations de votre algorithme. C’est aussi un défi de popularité, alors ne vous laissez pas aller à l’évidence echo(pi)et dites que la virgule flottante IEEE 754 arrondit les chiffres.

Le gagnant reçoit une quiche lorraine².

¹ Attention: pas réellement un score bonus. En réclamant le score, vous acceptez de me préparer une tarte avant la Journée Pi, 2016

² Avertissement: la quiche lorraine est utilisée comme une métaphore pour faire marquer votre réponse comme "acceptée"

Sanchises
la source
Connexes: lien
Sp3000
2
Je vote pour clore cette question hors sujet parce que les défis sournois ne sont plus au sujet ici. meta.codegolf.stackexchange.com/a/8326/20469
cat.

Réponses:

77

Algorithme

En utilisant le résultat bien connu:

entrez la description de l'image ici

on définit 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

Essai

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

Divulgacher

L' intégrale de Borwein est l'idée mathématique d'une blague. Alors que l'identité ci-dessus tient jusqu'à sinc (x / 13), la valeur suivante est en réalité:

entrez la description de l'image ici

Uri Granta
la source
12
Probablement l'une des meilleures réponses à une question sournoise de ces derniers temps.
Optimiseur
14
"L'idée mathématique d'une blague". +1
nourriture
16
C'en est une bonne! IIRC, l’une des blagues les plus connues avec cette intégrale a été lorsque l’on a enregistré les résultats jusqu’à l’étrange sur Wolfram Alpha et envoyé un rapport de bogue ... Ce que les développeurs de WA ont passé des siècles à essayer de corriger =)
Mints97
3
Cette référence donne une bonne explication de ce qui se passe.
TonioElGringo
59

Pour trouver pi, nous allons intégrer cette équation différentielle bien connue:

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

Avec une condition initiale

> 0 <y0 <2 * pi

Il est bien connu que ce problème de valeur initiale converge vers π lorsque t augmente sans limite. Donc, tout ce dont nous avons besoin est de commencer avec une estimation raisonnable pour quelque chose entre 0 et 2π, et nous pouvons effectuer une intégration numérique. 3 est proche de π, nous allons donc choisir y = 3 pour commencer.

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

Voici quelques résultats à chacun pour différents nombres d'étapes:

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

Comment ça fonctionne:

Cette équation différentielle est bien connue car il est extrêmement difficile de s’intégrer correctement. Alors que pour des valeurs t faibles, une intégration naïve produira des résultats acceptables, la plupart des méthodes d'intégration présentent une instabilité extrême car t devient très grand.

AJMansfield
la source
4
@UriZarfaty Cet article de Wikipédia l'explique plutôt bien: en.wikipedia.org/wiki/Stiff_equation
AJMansfield
1
Qu'est-ce que n? ...
Cole Johnson
1
@AJMansfield Je voulais dire: ce n'est déclaré nulle part. Votre fordécélération utilise t, mais votre boucle utilise n.
Cole Johnson
1
@ColeJohnson je viens de le réparer.
AJMansfield
2
Je pense que votre équation différentielle devrait se lire dy / dt = sin (y) * exp (t).
David Zhang
6

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

J'ai essentiellement découvert cette séquence par accident. Il commence au fur 1, 1et à mesure que chaque terme suivant s(n)est donné par s(n) = 2*s(n - 1) - s(n - 2) * (1 + m*m). Le résultat final est le plus petit ntel que s(n) < 0multiplié par 2m. Au fur et à mesure mque la taille diminue, la précision devrait être croissante.

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

Je suis à peu près sûr que ce sont des erreurs en virgule flottante (1 + m*m), mais je ne suis pas sûr. Comme je l'ai dit, je suis tombé par hasard dessus. Je ne suis pas sûr de son nom officiel. N'essayez pas cela avec un mtrop petit ou il fonctionnera pour toujours (si en 1 + m*m == 1raison d' mêtre si petit).

Si quelqu'un connaît le nom de cette séquence ou pourquoi elle se comporte comme ça, je l'apprécierais.

soktinpk
la source
Je pense que cela est dû à l'annulation, qui est une perte de chiffres lorsque l'on soustrait deux nombres presque égaux. S1 et s2 sont presque égaux après une itération.
Sanchises
1
Je ne sais pas encore comment cela fonctionne, mais cela me rappelle quelque chose que j'ai déjà fait: j'ai pris à plusieurs reprises la somme cumulée d'un signal bruyant et je l'ai normalisée en signifiant 0, valeur maximale 1. Cela convergerait vers une onde sinusoïdale, puisque c’est le seul signal qui possède sa propre anti-dérivée (avec un décalage de phase).
Sanchises
Je l'ai demandé à math.SE et j'ai eu cette réponse.
Sanchises