Nenad Božidarević

Računanje broja Pi Monte Karlo metodom


Pre nekoliko dana na predavanju iz Verovatnoće i statistike dotakli smo jednu zanimljivu temu – problem Bifonove igle (Buffon’s needle).

U ravni je dat skup paralelnih pravih, pri čemu je rastojanje između susednih pravih a, a > 0. Na ravan se baca igla dužine l, l < a. Odrediti verovatnoću da igla seče neku od pravih. ("Verovatnoća i statistika za inženjere i studente tehnike", Milan Merkle)

Bacanje igle kod Bifonovog problemaU pitanju je običan zadatak iz verovatnoće, ali ono što ga čini posebno zanimljivim jeste činjenica da se na osnovu baš ovakvog bacanja igle može odrediti broj π. Naime, zbog različitih uglova pod kojima igla može da se nađe u odnosu na pravu (kao na slici), u formuli za verovatnoću figuriše broj π. Ako znamo njegovu vrednost, možemo odrediti verovatnoću — ali ako je ne znamo, možemo baciti iglu veliki broj puta da bismo dobili verovatnoću, a onda preko nje izračunati π.

Naravno, jasno je da “veliki broj bacanja” uopšte nije lako izvesti. Ako želimo iole preciznije rezultate, nekoliko hiljada teško da će biti dovoljno, pa se mora stremiti ka daleko većim brojevima. Italijanski matematičar Mario Lazzarini je početkom 20. veka izveo ovaj eksperiment sa 3408 bacanja i dobio vrednost 355/113, što je 3.14159292. Za poređenje, vrednost broja π je 3.14159265, odnosno razlikuju se tek na sedmoj decimali. Članak na Wikipedii objašanjava zašto je dobijena ovako precizna vrednost sa malim brojem bacanja, ali nas to trenutno interesuje. Zanemarićemo njegove rezultate, i pokušati lično da izračunamo π uz pomoć kompjuterske simulacije koja nam omogućava da izvršimo milione “bacanja” u samo jednoj sekundi.

Pre nego što krenemo sa programiranjem, malo ćemo izmeniti problem da bismo olakšali njegovu implementaciju, ne umanjujući preciznost. Umesto beskonačne ravni i Bifonove igle, koristićemo kvadrat određenih dimenzija i krug upisan u njega.

Kvadrat i krugKao što nam je poznato, površina kvadrata stranice 2r je 4r2, dok je površina kruga upisanog u isti r2π. Ukoliko znamo obe površine, jednostavnom matematikom možemo izračunati π. Međutim, mi ih ne znamo, pa ćemo morati da ih aproksimiramo. Kao kod Bifonove igle, to ćemo učiniti nasumičnim biranjem tačaka koje pripadaju kvadratu, s tim što ćemo pamtiti koliko je tačaka pripadalo i krugu. Posle dovoljno iteracija, odnosi ukupnog broja tačaka i tačaka u krugu će biti približno jednaki odnosu površina, što ćemo onda upotrebiti da bismo izračunali π.

Kako je ovde u pitanju aproksimacija, logično je da će više bacanja davati bolje rezultate, i poništiti negativan efekat kompjuterski generisanih random brojeva koji nisu zaista nasumični, već pseudo-nasumični. Računanje broja π na ovaj način se zasniva na tome da je podjednaka verovatnoća izbora bilo koje tačke u okviru kvadrata, pa ćemo mi pretpostaviti da nam naš program daje tačke sa uniformnom raspodelom. Iako radimo samo sa celobrojnim koordinatama, pa zapravo ne analiziramo kompletnu površinu kvadrata, zbog velikih vrednost brojeva rezultati bi trebalo da budu dovoljno precizni.

Ovog načina računanja broja π Monte Karlovom metodom sam se setio zbog nekoliko ljudi koji su pohađali seminare u Petnici i kojima je zadatak bio da na ovaj način izračunaju π — jedina razlika je bila u tome što oni nisu koristili računare već papir i pirinač koji su posle prosipanja morali precizno da prebroje. Ovaj način preporučujem samo najupornijima.

Realizacija algoritma, sa druge strane, je veoma jednostavna. U svakoj iteraciji generišemo dve nasumične vrednosti koje nam predstavljaju X i Y koordinatu tačke, i njih koristimo da bismo izračunali da li ona pripada krugu ili ne. U svom kôdu sam radio sa velikim nenegativnim brojevima da bih izbegao korišćenje manje preciznih realnih brojeva, odnosno ostavio ga za sam kraj. Kvadrat je postavljen u donji levi ćošak prvog kvadranta koordinatnog sistema, odnosno obuhvata tačke od (0,0) do (a,a). Formula na kraju algoritma je dobijena preko odnosa površina kruga i kvadrata.

Formula za dobijanje broja Pi

Pošto koristimo random generator, a želimo što veću preciznost, ograničio sam stranicu kvadrata na najveći broj koji možemo dobiti od random generatora, što je u mom slučaju bilo 231-1, odnosno nešto preko dve milijarde. Takođe, da bi se centar kruga nalazio na celobrojnim koordinatama, stranicu kvadrata sam smanjivao na paran broj ukoliko je konstanta RAND_MAX neparna, što bi uglavnom i trebalo da bude slučaj.

Kôd

#include <stdio.h>
#include <stdlib.h>
#include <time.h>

int main()
{
	// Duzina stranice kvadrata

	unsigned long long a = RAND_MAX;

	if (RAND_MAX & 1) // Ukoliko je A neparno
		a--;

	// Poluprecnik kruga

	unsigned long long r = a/2;

	// Brojac tacaka u krugu

	unsigned long long pogodaka = 0;

	// Broj iteracija algoritma, tj. ukupan broj tacaka

	unsigned long long n;

	// Seed za random generator

	srand(time(0));

	// Ucitavamo broj iteracija

	printf("Uneti broj iteracija: ");
	scanf("%lld", &n);

	// Izvrsavamo algoritam

	unsigned long long x, y;
	unsigned long long rr = r*r;

	unsigned long long i;
	for (i = 1; i <= n; i++)
	{
		x = rand();
		y = rand();

		if (x > a || y > a) // Da ne bismo dobijali tacke van kvadrata
		{
			i--; // Ponavljamo iteraciju
			continue;
		}

		if ((x-r)*(x-r) + (y-r)*(y-r) <= rr) // Jednacina kruga
			pogodaka++;
	}

	// Ispisujemo dobijenu vrednost

	double pi = 4 * (double)pogodaka / n;

	printf("Aproksimacija broja Pi: %lf\n", pi);

	return 0;
}

Rezultati

Slede rezultati koje sam dobio testiranjem za različite brojeve iteracija, sve do jednog biliona — za ostalo nisam imao živaca.

1.000 iteracija -> 3,108000
10.000 iteracija -> 3,153600
100.000 iteracija -> 3,144000
1.000.000 iteracija -> 3,140188
10.000.000 iteracija -> 3,141856
100.000.000 iteracija -> 3,141841
1.000.000.000 iteracija -> 3,141715 (50 sekundi)
10.000.000.000 iteracija -> 3,141595 (6 minuta)
100.000.000.000 iteracija -> 3,141595 (sat vremena - bez poboljšanja)
1.000.000.000.000 iteracija -> 3,141594 (deset sati - malo poboljšanje)

Kao što možete videti, što je više iteracija to je precizniji, s tim što broj potrebnih iteracija za povećanje preciznosti stalno raste. Naravno, ponovno pokretanje programa ne bi dalo iste vrednosti, a postoji mogućnost i da se vrednost u nekom trenutku udalji od broja π umesto da mu se približi, ali takva je verovatnoća: ništa nije sigurno, samo verovatno.

Komentari (38)

Рајко Велимировић 6 years ago

Било би ми драго да број ПИ=3,141… срачунате преко обрасца који сам пронашао а који гласи
ПИ=180*m*sin(1/m)
m-врло велики број пример 1.0Е+100000000
sin-синус
180-број
Ја сам добио 25 милиона децимала на 12 хиљада страна А4
Рајко Велимировић

Odgovori

Nenad Božidarević 6 years ago

Odgovorio sam Vam na e-mail koji ste mi poslali — problem je u limesu izraza m*sin(1/m) koji teži jedinici, pa je vrednost izraza 180.

Odgovori

Рајко Велимировић 6 years ago

Нема нигде грешке 180 степени је једнако ПИ радиана
један степен је једнако ПИ кроз 180 радиана
Уосталом ако тражиш површину круга или обим, ПИ се добија преко горе објашњеног претварањем степена у радијане
Да је синус икс кроз икс једнако један када икс тежи нули и врапци знају.

Odgovori

Nenad Božidarević 6 years ago

Dakle broj 180 u formuli predstavlja stepene. U redu, krajnji rezultat od 180 stepeni jeste Pi, ali mi broj 180 ne možemo pretvoriti u radijane ukoliko nam nije već poznat broj Pi.

Odgovori

Rajko Velimirović 6 years ago

Једноставно ћемо овако јер ме ниси разумео
узми да је m= 11111111111111111111111 и убаци у обрасцу
ПИ=180*m*sin(1/m)
ево како то изгледа у обрасцу
ПИ=180*11111111111111111111111 *sin(1/11111111111111111111111 )
добићеш 3.1415926535897932384626433832795
Па ако нађеш краћи и елегантнији образац за ирационални број ПИ молим те јави ми

Odgovori

Nenad Božidarević 6 years ago

Ne možemo se osloniti na kalkulatore kada je ovo u pitanju, jer oni u svojim proračunima već koriste poznatu vrednost 3,14… Da, rezultat za veliko M će biti Pi (ukoliko su vrednosti brojeva stepeni, a ne radijani), ali samo zato što će kalkulator na kraju izvršiti konverziju iz 180 stepeni u Pi radijana. Sa druge strane, Pi ipak možemo dobiti bez ovih konverzija, i to upotrebom Tejlorovog polinoma za arccos.

Odgovori

Trololo 6 years ago

Kraci obrazac bi bio
arccos(-1)

Odgovori

Rajko Velimirović 6 years ago

arccos(-1)=180
видиш ја множим ПИ/180 само у граничној вредности док у обрасцу не
и то треба да разјасниш

Odgovori

Рајко Велимировић 6 years ago

Ја се нисам ослонио само на калкулатор већ имам и граничну вредност коју сам ти послао.А да калкулатор врши конверзију до данас нисам знао.Ипак се држи Монте Карлове методе она је најбоља.

Odgovori

Vladimir Dimić 6 years ago

Gospodine Rajko, Vaša metoda nije matematički korektna. Korišćenjem Vašeg obrasca ПИ=180*m*sin(1/m) dobijemo neku vrednost izraženu u stepenima. Za neko veliko m, dobije se da je izraz približno jednak 180 stepeni. E sada, da biste to pretvorili u radijane, treba da dobijenu vrednost pomnožite sa Pi/180°. To, naravno, ne možete da uradite, jer vi u računanju broja Pi koristite isti taj broj Pi čiju vrednost ne znate. Zapravo, dobijete izraz Pi=180°*Pi/180°, što predstavlja identitet, a ne izraz za računanje broja Pi.

Odgovori

Рајко Велимировић 6 years ago

Апсолутно ме нисте разумели ево овако ћемо
образац Пи=180*m*sin(1/m) не треба множити са ПИ/180
узмите m=22222222222222222222222 замените у образац без множења са ПИ/180
или то изгледа овако нумерички
ПИ=180*22222222222222222222222*sin(1/22222222222222222222222)
Добићете ПИ у једном потезу шаховски речено
за добијање већих резултата преко мог обрасца препоручујем
http://harry-j-smith-memorial.com/index.html
калкулатор је бесплатан .Ја сам добио 25 милиона децимала тачно али за 12 сати рада рачунара.

Odgovori

Trololo 6 years ago

“ПИ=180*22222222222222222222222*sin(1/22222222222222222222222)”
180*22222222222222222222222*sin(1/22222222222222222222222) = 180.
Znaci Pi = 180.
Posto je ovo vrednost u stepenima, kako izracunate vrednost u radijanima? Kako od 180 dobijete 3.14…. ?

Odgovori

Рајко Велимировић 6 years ago

180*22222222222222222222222*sin(1/22222222222222222222222)=3.1415926535897932384626433832795 Јуначе како доби 180 објасни ми молим те

Odgovori

Vladimir Dimić 6 years ago

Prvo, čitajući ove komentare, došao sam do zaključka da Vi, uz svo dužno poštovanje, protivrečite sami sebi. Najpre u Vašem drugom komentaru kažete “ПИ се добија преко горе објашњеног претварањем степена у радијане“. Zatim, deo Vašeg odgovora na moj komentar glasi “Пи=180*m*sin(1/m) не треба множити са ПИ/180“. Jedno pitanje za Vas: Šta za Vas znači pretvaranje iz stepena u radijane?
Drugo, sam izraz nije validan. ПИ=180 * m * sin(1/m) jeste jednako 180 kad m teži beskonačnosti, a složićete se da Pi nije jednako 180. Ono što fali u Vašem izrazu je oznaka za stepen, tj. ПИ=180° * m * sin(1/m). Tada vrednost izraza jeste jednaka 180°, što jeste jednako Pi. Ali taj izraz nikako nije formula pomoću koje se može izračunati Pi, jer se on na kraju svodi na Pi=Pi.

Odgovori

Рајко Велимировић 6 years ago

Тако је уверили сте ме сад ме је срамота

Odgovori

Рајко Велимировић 6 years ago

За господина Tropolo рачун је тачан и јесте 180 али синус угла срачунај преко степена добио си 180 јер си угао рачунао у радиане ето изгледа да се ту неразумемо надам се сада да је све јасно.

Odgovori

Nenad Božidarević 6 years ago

Radijani ili ne, m*sin(1/m) je, kada M teži beskonačnosti, u svakom slučaju 1 i samim tim ni ne utiče na krajnju vrednost.

Odgovori

Рајко 6 years ago

Нити си меродаван да дајеш овакав коментар нити ти је то неко тражио.

Odgovori

Nenad Božidarević 6 years ago

Je li ni oni “doktori matematike koji sumnjaju u graničnu vrednost” nisu merodavni?

Odgovori

Trololo 6 years ago

Gospodine Rajko, da li mozete pokazete vasu formulu bez koriscenja stepena kao merne jedinice?

Odgovori

Рајко 6 years ago

Шта мислиш шта сам ја радио пуних шест година од кад сам пронашао не само овај образац већ још пуно других ствари.Твоје мишљење апсолутно ми ништа незначи било оно позитивно или негативно .Да погодим доктори са коима си ти разговарао су Wolfram|Alpha .Интересантно да те уопште није занимао поступак добијања обрасца .Знаш која је разлика између мене и тебе ако те занима рећи ћу ти

Odgovori

Trololo 6 years ago

Gospodine Rajko, molim vas da nam kazete razliku.

Rado bih video i nacin na koji ste dosli do obrasca.

Takodje, ja bih voleo da mi pokazete vas obrazac u formi koja se ne oslanja na racunanje sinusa ugla preko stepena.

Odgovori

Vladimir Dimić 6 years ago

Gospodine Rajko, želeo bih da Vam kažem nekoliko stvari. Voleo bih da pročitate ovaj dugačak komentar.

  1. Probao sam Vaš obrazac u kalkulatoru koji ste naveli. U pravu ste da se za neko veliko m dobije broj koji teži Pi. Ali jednu stvar ste prevideli – ako odete u podešavanja kalkulatora, videćete da on podrazumevano radi sa stepenima. Dakle, konvertuje vrednost iz stepena u radijane, a znamo da se to radi množenjem sa Pi/180. Takođe znamo da ovde Vaš obrazac pada u vodu (već sam obrazložio zašto).
  2. Guglajući Vaše cenjeno ime, našao sam pozamašan broj članaka na temu broja Pi na različitim sajtovima (na srpskom i ruskom jeziku), na kojima ste Vi ostavljali komentare o tome kako ste pronašli novi obrazac za broj Pi. Obrazac se razlikuje od ovog ovde. Nažalost, moraću da Vas razočaram – ni taj obrazac ne valja iz istog razloga kao i ovaj koji ste naveli ovde.
  3. Jasno mi je zašto ne možete da prihvatite da je Vaš obrazac netačan. Potrošili ste šest godina na njega (i na ostale stvari, kao što ste i sami rekli). Razumljivo je da ćete Vaše čedo koje ste gajili šest godina štititi svim snagama.
  4. Vi, kao akademski građanin, bi trebalo da znate osnove pravopisa srpskog jezika. Na primer, da se negacija uz glagol piše odvojeno. I ne budite škrti na znacima interpunkcije, umnogome pomažu da se Vaš cenjeni komentar bolje razume. Posle tačke za kraj rečenice se stavlja razmak, a ne pre. Takođe, ton kojim izlažete svoj stav, naročito u par komentara iznad ovog, svakako ne priliči odraslom čoveku.
  5. Kažete da ste dobili 25 miliona decimala broja Pi. Svaka Vam čast. Ali, moraću da Vas razočaram. Poznato je mnogo više od toga, po nekim podacima – 2,699,999,990,000 decimala. Tako da ne vidim razlog za izmišljanjem novih obrazaca za nalaženje ovog broja (pri tom, netačnih). Obrasci već postoje, problem je da se nađe dovoljno brz računar koji će u nekom razumnom roku izračunati Pi sa zadatom tačnošću.

I za kraj Vam želim sve najbolje u Vašem životu i radu, i preporučujem da Vaše slobodno vreme trošite na konstruktivnije stvari.

Odgovori

Рајко 6 years ago

Врло сам задовољан и захвалан Вашим одговором заиста.
Пробаћу да одговорим на овај начин
број m третирајте као променљиву
за m=10
3.14159
за m=100
3.1415926
за m=1000
3.141592653
за m=10 000
3.14159265358
за m=100 000
3.1415926535897
за m=1000 000
3.141592653589793
за m=10000000
3.14159265358979323
за m=100000000
3.1415926535897932384
за m=1000000000
3.141592653589793238462
Закључак за различите вредности променљиве m добијам различиту тачност ПИ броја
где је онда ту ПИ/180 како упорно тврдите сво време.
Осим тога када кажете степени постоје још минути и секунде.
Сви уређаји за мерење углова барем код нас су у степенима,минутама и секундама.Поново Вам кажем Ви нисте меродавни за исправност мога чеда.

Odgovori

Nenad Božidarević 6 years ago

Ono što mene zanima jeste kako biste izračunali vrednost vašeg obrasca na papiru, bez upotrebe kalkulatora, za recimo m = 1000? Tačnije, ako nam je ugao dat u stepenima, što ste naglasili da jeste slučaj, kako ćemo izračunati vrednost sinusa?

Odgovori

Trololo 6 years ago

PI/180 je u kalkulatoru koji koristite jer ste podesili Degree mode.

Odgovori

Vladimir Dimić 6 years ago

Kada u kalkulatoru računate sinus, kalkulator koristi Tejlorov red za sinus da bi našao njegovu vrednost. Ukoliko mu zadate ugao u stepenima, on automatski vrši konverziju u radijane, pa tek onda tako dobijenu vrednost ugla koristi u Tejlorovom redu. Eto tu se krije to množenje sa Pi/180.

Odgovori

Рајко 6 years ago

Господо честитам али признајте да нисте помињали Тајлеров ред тако сам мислио на друге ствари и стога се извињавам.
Ако интегралним рачуном хоћемо да добијемо обим или површину круга
резултат је у степенима ништа ново значи нисмо добили број ПИ
већ смо га увезли претварањем ПИ/180
ако поставимо себи задатак да нађемо површину и обим круга без увожења броја ПИ и успемо у томе онда је то у мом скромном раду урађено.Имам два идентична низа са заједничким константама и елиптички интеграл.А онај други образац Пи=n/2(sin(360/n)) где је n{3,4,5,6……….бесконачно} исто мој
образац и тиме хоћу да кажем да има бесконачно ПИ бројева зависно од променљиве n на пример за n=3 ради се о једнакостраничном троуглу чија је површина R на квадрат пута ПИ3 при чему је ПИ3=3/2(sin360/3)) R је полупречник описаног круга и тако за читав низ геометријских за квадрат ПИ4=2 цео број и за дванаестоугао Пи12=3 за што је библија у праву када каже да је ПИ=3 али не потпуно.Сличан низ имам код елипсе итд. Момци ипак се човек може ослонити на ВАС.

Odgovori

Рајко 6 years ago

Јуче сам био много уморан и извините на грешкама.Мислим да је најбоље да Вам пошаљем рад у потпуности .

Odgovori

Vladimir Dimić 6 years ago

Da, to bi bilo najbolje.

Odgovori

Рајко Велимировић 6 years ago

Ја једино знам Email гдина Божидаревића па ако ће он проследити Вама онда је проблем решен.

Odgovori

Vladimir Dimić 6 years ago

Imate moju adresu kad kliknete na moje ime, pa na skroz dnu strane.

Odgovori

Рајко 6 years ago

Послато
Користим прилику да се захвалим свим учесницима у дискусији и жао ми је што није било више учесника.

Odgovori

Milojko Pantic 6 years ago

Rajko mi te volimo, Rajko mi te volimo! Vrati se Rajko, vrati se Rajko!

Odgovori

kizami 4 years ago

Kolika ti je granulacija kod ovog algoritma?

Odgovori

Nenad Božidarević 4 years ago

Ne znam na šta misliš kada kažeš “granulacija algoritma”, ali koristim celobrojnu rešetku za random tačke i kvadrat širine 2^32.

Odgovori

marko 4 years ago

Dragi Rajko,

Formula vam je netacna!!! ALI za svako postovanje je sto rekreativno pokusavate da se bavite matematikom i mozda stvarno nekad ,ako dopunite svoje znanje, uspete nesto da dokazete.

Imajte na umu da se naucnici nikad nisu bavili naukom da bi postali slavni vec da bi uposlili mozak jer im je dosadjivala svakodnevna prazna prica i jer je najbolji osecaj kapirati stvari.

Takodje imajte na umu da mnogi brilijantni umovi ,iako su se dali celog sebe, nikad nista nisu dokazali jer nisu imali srece,mozda je samo 100 ljudi u celokupnoj istoriju covecanstva otkrilo 99% matematike.

Preporucicu vam ovaj sajt za besplatno skidanje knjiga http://gen.lib.rus.ec/

Malo progutajte ponosa i zaboravite ovu formulu jer je netacna, morate prvo da se ozbiljnije pozabavite limesima =D preporucujem sve od james stewart =D

Pozdrav, sve najbolje :) !!!

Odgovori

Postavi komentar