Расчет PI с использованием метода Монте-Карло дает неточный ответ

Question

Я пытаюсь вычислить PI с использованием метода Монте-Карло. Мой код дает результат 3.000 независимо от того, насколько большой MAXLEN. После многократного отладки я не мог получить то, что я делаю неправильно.

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

#define sqr2(x) ((x)*(x)) 
#define frand() ((double) rand()/(RAND_MAX)) 
#define MAXLEN 1000 


int circumscribed(int radius){ 
    float xcoord = frand(); 
    float ycoord = frand(); 
    float coord = sqr2(xcoord) + sqr2(ycoord); 

    if(coord <= radius) 
     return 1; 
    return -1;  
} 

int main() 
{ 
    int i; 
    int circles = 0, rect = 0;; 
    for(i = 0; i < MAXLEN; i++) 
    { 
     if(circumscribed(1) > 0) // if(circumscribed(1)) shoul be enough but it doesn't work. Very odd in my opinion. 
      circles++; 
     rect++; //this is equal to MAXLEN, I just used it for debugging 
    } 

    float PI = 4 * circles/rect; 
    printf("PI is %2.4f: \n", PI); 
    return 0; 
} 

Answer 1

С circles и rect оба int, результатом 4 * circles/rect будет int. Вместо этого используйте числа с плавающей запятой.

float PI = 4.0 * (float)circles/rect; 

Answer 2

Это выражение все целые числа:

4 * circles/rect; 

Таким образом, результатом является целое число (3 в этом случае).

(как аналогичный пример: 10/3 == 3, но 10.0/3.0 == 3.333333)

Попробуйте вместо этого:

4.0 * circles/rect; 

Просто изменяя (int)4 к (double)4, обращаясь к ней, как 4.0 или даже 4. должно быть достаточно.

---

Другие Наблюдения Разное

Эта линия имеет дополнительный запятой:

int circles = 0, rect = 0;; 

Ваша функция circumscribed использует float. Ваша переменная PI также является float.
Если вы используете double, вы получите большую точность.

Answer 3

Вы целочисленной математикой здесь:

float PI = 4 * circles/rect; 

Изменение 4 для 4.0 достаточно, чтобы решить эту проблему:

float PI = 4.0 * circles/rect; 

Рабочая версия here

Answer 4

Вы делаете integer Math здесь. circles и rect являются int в вашем коде, поэтому результат 4 * circles/rect также является int. Вместо этого используйте числа с плавающей запятой. Используйте double для лучшей точности.

double PI = 4.0 * (double)circles/rect; 

Answer 5

float PI = 4 * circles/rect; 

Это выполняет целочисленную арифметику справа от = таким образом, ограничивая свой результат до 1 значащей цифры.

Вместо:

double PI = 4.0 * circles/rect; // Best 

---

Подробности

float PI = 4.0 * (float)circles/rect; // OK 

рекомендуют избегать использования float. Используйте вместо этого double, чтобы избежать ограничения этого Монте-Карло примерно на 7 цифр. Поскольку circles становится большим, он не может точно конвертироваться в то же значение float. Вместо этого он преобразуется в округленный float. Это происходит примерно на circles> 8 000 000 на многих машинах. При округлении вы без необходимости ограничиваете достижимую точность. Используя double, это округление не происходит до тех пор, пока circles не будет около 9e15.

float PI = 4.0 * (double) circles/rect; // Better 

Явный слепок (double) circles может быть полезным для читателя, но код будет работать так же, как и без него. 4.0 является double и приведет к продвижению double до circles, прежде чем произойдет множественность, даже без приведения.

double PI = 4.0 * circles/rect; // Better 

Как 4.0 * ... результат является двойным, лучше точность сохраняется за счет экономии на двойной. Использование float PI вызывает умножение/деление, которое было выполнено с двойной точностью, чтобы уменьшить его точность при сохранении в поплавке.

float or double PI = ... 
printf("PI is %2.4f: \n", PI); 

Примечание: Здесь, PI преобразуется в двойной передается printf() независимо, если PI был объявлен float или double. Таким образом, можно получить доступ к точности в PI, объявив его двойным.

---

Примечание: int, float и double range и precision зависят от машины. Вышеизложенное отражает общую реализацию.

Answer 6

Проверьте этот код:

#include <stdlib.h> 
#include <stdio.h> 
#include <math.h> 
#include <string.h> 
#define SEED 35791246 

int main(int argc, char** argv) 
{ 
    int niter=0; 
    double x,y; 
    int i,count=0; /* # of points in the 1st quadrant of unit circle */ 
    double z; 
    double pi; 

    printf("Enter the number of iterations used to estimate pi: "); 
    scanf("%d",&niter); 

    /* initialize random numbers */ 
    srand(SEED); 
    count=0; 
    for (i=0; i<niter; i++) { 
     x = (double)rand()/RAND_MAX; 
     y = (double)rand()/RAND_MAX; 
     z = x*x+y*y; 
     if (z<=1) count++; 
     } 
    pi=(double)count/niter*4; 
    printf("# of trials= %d , estimate of pi is %g \n",niter,pi); 

return 0; 
}