Obliczanie liczby pi metodą Monte Carlo w C++

Algorytm rozwiązywania równania liniowego

Schemat algorytmu obliczania liczby pi metodą Monte Carlo

#include <iostream>
#include <iomanip>
#include <cstdlib>
#include <ctime>

using namespace std;

int main(int argc, char *argv[]) {
  int n;
  cout<<"Liczba losowanych punktow n = ";
  cin>>n;
  srand(time(NULL));
  float x,y;
  int k=0;
  for(int i=1; i<=n; i++){
    x=rand()/float(RAND_MAX);
    y=rand()/float(RAND_MAX);
    if(x*x+y*y<=1) k++;
  };
  float p=4.*k/n;
  cout<<"pi = "<<setprecision(10)<<p<<endl;
  char c;
  cin>>c;
}

 

 W pierwszych czterech liniach programu dołączamy pliki nagłówkowe bibliotek wykorzystywanych w programie: iostream (strumieni wejścia-wyjścia), iomanip (manipulowanie strumieniami i/o, wykorzystamy formatowanie strumienia wyjścia), cstdlib (biblioteka standardowa języka C zawierająca funkcje generatora liczb pseudolosowych), ctime (biblioteka obsługi czasu języka C). W szóstej linii informujemy kompilator, że będziemy używać w programie nazw zdefiniowanych w przestrzeni nazw std.

 W języku C++ program rozpoczyna się od wykonania funkcji main. W naszym przykładzie jest to jedyna funkcja programu. Linia 8 zawiera początek definicji funkcji main, jest to początek realizacji naszego algorytmu (blok START (1)). W kolejnej linii deklarujemy zmienną całkowitą n, którą za chwilę będziemy używali. W C++ zmienne można deklarować w dowolnym miejscu przed ich użyciem, dlatego nie deklarujemy tutaj od razu zmiennych i, k, p, które wystąpią w dalszej części programu.

  Linie 10 i 11 realizują operację wejścia/wyjścia (2), w której pobierane są dane wejściowe - liczba losowanych punktów.

  W linii 12 wywołujemy funkcję biblioteki standardowej srand. Jest to funkcja inicjująca wbudowany generator liczb pseudolosowych wartością podaną jako argument wywołania funkcji. Jako argument podajemy wynik działania funkcji time(NULL). Funkcja time() jest zawarta w bibliotece ctime, Wywołanie jej z argumentem NULL zwraca w wyniku bieżący czas systemowy, który oczywiście przy każdym uruchomieniu programu będzie inny. Gdybyśmy nie wywołali tej funkcji program przy każdym uruchomieniu generowałby ten sam ciąg liczb pseudolosowych i w rezultacie dla określonej liczby losowanych punktów otrzymywalibyśmy zawsze ten sam wynik. W każdym programie w C++, w którym wykorzystujemy generator liczb pseudolosowych powinniśmy jeszcze przed pierwszym losowaniem raz wywołać funkcję srand(time(NULL)).

 W linii 13 deklarujemy dwie zmienne zmiennoprzecinkowe (float), w których będziemy zapisywać współrzędne losowanych punktów, a w kolejnej linii definiujemy i inicjujemy wartością 0 zmienną k - licznik punktów leżących w kole o promieniu 1 (3).

  Linia 15 rozpoczyna instrukcję pętli for. Instrukcje zawarte pomiędzy nawiasami klamrowymi { } zostaną wykonane n razy (zmienna sterująca i pętli zmienia się od 1 do n, za każdym razem zwiększając swą wartość o 1). Linia 11 realizuje część punktu (3) algorytmu inicjując na początku zmienną i wartością 1 oraz punkty (7) i (8) algorytmu, zwiększając po każdym wykonaniu instrukcji pętli zmienną sterującą i o 1 i sprawdzając, czy wartość liczby i nie przekroczyła jeszcze górnej granicy pętli (i<=n). Wewnątrz pętli for wykonywana jest zawsze dokładnie jedna instrukcja. Jeżeli potrzebujemy w tym miejscu, tak jak ma to miejsce w naszym przykładzie, więcej instrukcji, musimy je zapisać pomiędzy nawiasami klamrowymi { } jako instrukcję złożoną.

 Instrukcje wykonywane w pętli zawarte są w liniach 16 - 18. Dla każdego losowanego punktu najpierw losujemy jego współrzędne x i y - punkt (4) algorytmu. Do losowania wykorzystujemy funkcję biblioteczną rand(), której wynikiem jest liczba losowa (pseudolosowa) typu int z zakresu od 0 do RAND_MAX. RAND_MAX jest stałą zdefiniowaną w bibliotece cstdlib i określa maksymalną liczbę całkowitą generowaną przez generator liczb pseudolosowych. Gdy wylosowaną liczbę całkowitą podzielimy przez jej maksymalną wartość (RAND_MAX) otrzymamy w wyniku liczbę losową z zakresu od 0 do 1. Jednak w języku C++ wynik dzielenia dwóch liczb całkowitych jest zawsze liczbą całkowitą. Jeśli chcemy otrzymać liczbę losową typu float z zakresu [0, 1] musimy na kompilatorze wymusić dzielenie zmiennoprzecinkowe, przekształcając (rzutując) przynajmniej jeden ze składników operacji dzielenia na liczbę zmiennoprzecinkową. float(RAND_MAX) rzutuje przed wykonaniem operacji dzielenia stałą całkowitą RAND_MAX na typ float.

 W linii 18 sprawdzamy, czy punkt o wylosowanych współrzędnych x, y leży w kole o promieniu r=1 (5). Jeśli tak to dalej w tej samej linii zwiększamy zawartość zmiennej k o 1 (6).

 Po zakończeniu pętli zmienna n zawiera całkowitą liczbę wylosowanych punktów (zawartość zmiennej n nie uległa zmianie), a zmienna k zawiera informację o liczbie wylosowanych punktów leżących w kole o promieniu 1.

 W linii 20 deklarujemy zmienną p i przypisujemy do niej przybliżoną wartość liczby π (9), a następnie w linii 21 wypisujemy otrzymany wynik (blok wejścia-wyjścia (10)).

 Instrukcje w liniach 22 - 23 nie mają żadnego związku z wykonaniem algorytmu dzielenia. Służą do zatrzymaniu programu do czasu wprowadzenia dowolnego znaku i naciśnięcia klawisza enter, co zapobiega natychmiastowemu zamknięciu okienka z programem.

 Nawias klamrowy } z linii 24 kończy funkcję main i zarazem działanie całego programu (blok końca algorytmu (11)).