Quiche Lorraine [ditutup]

52

Karena itu adalah Pi hari baru-baru ini, saya telah melihat sebuah nomor dari tantangan yang meminta Anda untuk menghitung pi.

Tentu saja, quorum lorraine tidak terlalu banyak (Anda dapat mengklaim Skor Bonus¹ dari +1 jika Anda menebak tantangan di luar judul). Dengan demikian, tugas Anda adalah menulis algoritma atau metode yang kelihatannya mendekati Pi pada pandangan pertama, tetapi dijamin tidak akan menyatu ke arah Pi.

Ini adalah tantangan curang, jadi pastikan itu akan menghasilkan 3,14 ... untuk kasus uji sederhana, misalnya dengan 10 iterasi dari algoritma Anda. Ini juga merupakan tantangan popularitas, jadi jangan pergi untuk yang jelas echo(pi)dan mengatakan bahwa IEEE 754 floating point rounds beberapa digit naik atau turun.

Pemenang mendapat quiche lorraine².

¹ Peringatan: sebenarnya bukan skor bonus. Dengan mengklaim skor, Anda setuju untuk membuatkan saya pai sebelum Pi Day, 2016

² Peringatan: quiche lorraine digunakan sebagai metafora karena jawaban Anda ditandai sebagai 'diterima'

Sanchises
sumber
Terkait: tautan
Sp3000
2
Saya memilih untuk menutup pertanyaan ini sebagai di luar topik karena tantangan curang tidak lagi di-topik di sini. meta.codegolf.stackexchange.com/a/8326/20469
cat

Jawaban:

77

Algoritma

Menggunakan hasil yang terkenal:

masukkan deskripsi gambar di sini

kami mendefinisikan dengan 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

Pengujian

>>> 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 terpisahkan adalah ide matematika tentang lelucon praktis. Sementara identitas di atas berlaku hingga sinc (x / 13), nilai selanjutnya sebenarnya:

masukkan deskripsi gambar di sini

Uri Granta
sumber
12
Mungkin salah satu jawaban terbaik untuk pertanyaan curang belakangan ini.
Pengoptimal
14
"Gagasan matematika tentang lelucon praktis". +1
hapus klasifikasi
16
Itu bagus! IIRC, salah satu lelucon yang paling terkenal dengan integral ini adalah ketika seseorang mencatat hasilnya hingga yang aneh di Wolfram Alpha, dan mengirim laporan bug ... Yang WA menghabiskan waktu
bertahun-tahun
3
Referensi ini memberikan penjelasan yang baik tentang apa yang sedang terjadi.
TonioElGringo
59

Untuk menemukan pi, kami akan mengintegrasikan persamaan diferensial yang terkenal ini:

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

Dengan kondisi awal

> 0 <y0 <2 * pi

Sudah diketahui bahwa masalah nilai awal ini konvergen menjadi π karena t bertambah tanpa batas. Jadi, yang kita butuhkan hanyalah memulai dengan tebakan yang masuk akal untuk sesuatu antara 0 dan 2π, dan kita dapat melakukan integrasi numerik. 3 dekat dengan π, jadi kami akan memilih y = 3 untuk memulai.

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

Berikut adalah beberapa hasil di setiap langkah untuk jumlah langkah yang berbeda:

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

Bagaimana itu bekerja:

Persamaan diferensial itu terkenal karena sangat sulit untuk diintegrasikan dengan benar. Sementara untuk nilai t kecil, integrasi naif akan menghasilkan hasil yang dapat diterima, sebagian besar metode integrasi menunjukkan ketidakstabilan ekstrim karena t menjadi sangat besar.

AJMansfield
sumber
4
@UriZarfaty Artikel wikipedia ini menjelaskannya dengan cukup baik: en.wikipedia.org/wiki/Stiff_equation
AJMansfield
1
Apa itu n? ...
Cole Johnson
1
@ AJMansfield yang saya maksud: tidak dinyatakan di mana pun. forDeselerasi Anda digunakan t, tetapi loop Anda digunakan n.
Cole Johnson
1
@ ColeJohnson Saya baru saja memperbaikinya.
AJMansfield
2
Saya pikir persamaan diferensial Anda harus membaca dy / dt = sin (y) * exp (t).
David Zhang
6

Kode:

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

Saya pada dasarnya menemukan urutan ini secara tidak sengaja. Ini dimulai sebagai 1, 1dan setiap istilah setelah itu s(n)diberikan oleh s(n) = 2*s(n - 1) - s(n - 2) * (1 + m*m). Hasil akhirnya adalah yang terkecil nyang s(n) < 0dikalikan dengan 2m. Semakin mkecil, itu harus menjadi lebih dan lebih akurat.

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

Saya cukup yakin ini adalah kesalahan floating point karena (1 + m*m)mendekati satu, tapi saya tidak yakin. Seperti yang saya katakan, saya sengaja menemukan ini. Saya tidak yakin dengan nama resminya. Jangan coba ini dengan yang mterlalu kecil atau akan berjalan selamanya (jika 1 + m*m == 1karena mterlalu kecil).

Jika ada yang tahu nama urutan ini atau mengapa berperilaku seperti ini, saya akan sangat menghargainya.

soktinpk
sumber
Saya pikir ini karena pembatalan, yang merupakan kehilangan digit ketika mengurangi dua angka yang hampir sama. S1 dan s2 hampir sama setelah iterasi.
Sanchises
1
Saya belum mengetahui bagaimana cara kerjanya, tetapi itu mengingatkan saya pada sesuatu yang pernah saya lakukan: Saya berulang kali mengambil jumlah kumulatif dari sinyal berisik dan menormalkannya menjadi 0, nilai maks 1. Ini akan menyatu menjadi gelombang sinus, karena itulah satu-satunya sinyal yang merupakan anti-turunannya sendiri (dengan fase offset).
Sanchises
Saya menanyakannya di math.SE, dan mendapat jawaban ini .
Sanchises