モンテカルロ法

課題12-1 モンテカルロ法による円周率

モンテカルロ法により円周率を求めるプログラムです。

プログラム名 k1201.c

 1: /* モンテカルロ法による円周率 k1201.c */
 2: #include <stdio.h>
 3: #include <stdlib.h>
 4: #include <time.h>
 5: #define IMAX 10000000 
 6:
 7:int main()
 8: {
 9:     int naibu =0;
10:     int i=0;
11:     double x,y,r;
12:     srand( (unsigned int)time( NULL ) ); /*乱数初期化*/
13:     printf ("%d \n",RAND_MAX);
14:
15:     for ( i=1; IMAX>i; i++ ){
16:         x = rand()/(RAND_MAX+1.0); 
17:         y = rand()/(RAND_MAX+1.0); 
18:         r = x*x + y*y;
19:         if(1.0>r) naibu++;
20:         printf ("%f %f %d %d %18.16f \n",x,y,naibu,i,(double)naibu/i*4 );
21:     }
22: }

プログラム解説 k1201.c

3: #include <stdlib.h>

乱数のrand(),srand()を使うためのヘッダです。

4: #include <time.h>

乱数の種として時刻を使います。time()のためのヘッダです。

12: srand( (unsigned int)time( NULL ) );

時刻を種にすることで、毎回違う乱数になるようにします。

16: rand()

乱数をひとつ発生させます。C言語では[0,RAND_MAX] つまり 0以上 RAND_MAX 以下の整数です。RAND_MAXはシステムにより異なりますが、ここでは 2147483647 でした。

(RAND_MAX+1.0)で割ることで[0,1) つまり 0以上 1未満の数にしています。1でなく1.0なのは浮動小数点数にしないと整数のままで計算してしまい、[0,1)内の整数は0しかなくなってしまうからです。

(double)RAND_MAXで割ることで[0,1] つまり 0以上 1以下の数にできます。なぜこちらを使わないのかは考えてみてください。

18: r = x*x + y*y;

rは中心からの距離の二乗です。これが1.0より小さければ距離も1.0より小さいのでルートを計算しません。

実行結果 k1201.c

実行結果はこうなります。

~/c$ ./k1201.c
2147483647    ←RAND_MAX
0.437073 0.545812 1 1 4.0000000000000000 
0.841015 0.125503 2 2 4.0000000000000000 
0.559205 0.641552 3 3 4.0000000000000000 
0.458749 0.594224 4 4 4.0000000000000000 
0.498889 0.492063 5 5 4.0000000000000000 
0.047432 0.534532 6 6 4.0000000000000000 
0.618161 0.338228 7 7 4.0000000000000000 
0.343270 0.866057 8 8 4.0000000000000000 
0.522176 0.470137 9 9 4.0000000000000000 
0.419448 0.243027 10 10 4.0000000000000000 
0.000667 0.725147 11 11 4.0000000000000000 
0.237668 0.219389 12 12 4.0000000000000000 
0.038276 0.012017 13 13 4.0000000000000000 
0.671066 0.926532 13 14 3.7142857142857144 
0.028983 0.740724 14 15 3.7333333333333334 
0.755007 0.466057 15 16 3.7500000000000000 
0.286537 0.596022 16 17 3.7647058823529411 
0.591560 0.845742 16 18 3.5555555555555554 
0.237574 0.050309 17 19 3.5789473684210527 
0.439966 0.736463 18 20 3.6000000000000001 
0.542372 0.487398 19 21 3.6190476190476191 
0.270995 0.160533 20 22 3.6363636363636362 
0.825626 0.614265 20 23 3.4782608695652173 
0.026590 0.347801 21 24 3.5000000000000000 
0.084401 0.446038 22 25 3.5200000000000000 
0.590829 0.085069 23 26 3.5384615384615383 
0.171185 0.828497 24 27 3.5555555555555554 
0.304458 0.209461 25 28 3.5714285714285716 
0.840513 0.975524 25 29 3.4482758620689653 
0.135993 0.869496 26 30 3.4666666666666668 
0.716248 0.891000 26 31 3.3548387096774195 
0.335553 0.002784 27 32 3.3750000000000000 
0.487022 0.927113 27 33 3.2727272727272729 
0.848526 0.724595 27 34 3.1764705882352939 
0.977422 0.288493 27 35 3.0857142857142859 
0.461058 0.519793 28 36 3.1111111111111112 
0.775891 0.732053 28 37 3.0270270270270272 
0.680326 0.601516 29 38 3.0526315789473686 
0.346318 0.706916 30 39 3.0769230769230771 
0.949318 0.430719 30 40 3.0000000000000000 
0.152955 0.540147 31 41 3.0243902439024390 
0.515788 0.324140 32 42 3.0476190476190474 
0.368643 0.820246 33 43 3.0697674418604652 
0.533601 0.209156 34 44 3.0909090909090908 
0.795769 0.669594 34 45 3.0222222222222221 
0.078653 0.512017 35 46 3.0434782608695654 
0.560593 0.414206 36 47 3.0638297872340425 
0.514801 0.047615 37 48 3.0833333333333335 
0.341319 0.363328 38 49 3.1020408163265305 
0.772210 0.318741 39 50 3.1200000000000001 
0.651820 0.233268 40 51 3.1372549019607843 
0.838534 0.427711 41 52 3.1538461538461537 
0.965321 0.518860 41 53 3.0943396226415096 
0.029227 0.311639 42 54 3.1111111111111112 
0.225777 0.978545 42 55 3.0545454545454547 
0.742358 0.378731 43 56 3.0714285714285716 
0.518692 0.258146 44 57 3.0877192982456139 
0.702871 0.887335 44 58 3.0344827586206895 
0.078391 0.236472 45 59 3.0508474576271185 
0.096492 0.874161 46 60 3.0666666666666669 
0.906066 0.175145 47 61 3.0819672131147540 
0.386178 0.466659 48 62 3.0967741935483870 
0.589351 0.900979 48 63 3.0476190476190474 
0.514274 0.930670 48 64 3.0000000000000000 
0.264307 0.286484 49 65 3.0153846153846153 
0.249411 0.916127 50 66 3.0303030303030303 
0.519752 0.087945 51 67 3.0447761194029850 
0.343838 0.485074 52 68 3.0588235294117645 
0.606806 0.373066 53 69 3.0724637681159419 
0.796713 0.832582 53 70 3.0285714285714285 
0.351611 0.539071 54 71 3.0422535211267605 
0.211314 0.870303 55 72 3.0555555555555554 
0.797217 0.914185 55 73 3.0136986301369864 
0.757638 0.875608 55 74 2.9729729729729728 
0.150657 0.854130 56 75 2.9866666666666668 
0.749769 0.056723 57 76 3.0000000000000000 
0.029275 0.135946 58 77 3.0129870129870131 
0.523382 0.618626 59 78 3.0256410256410255 
0.036925 0.037656 60 79 3.0379746835443040 
0.549296 0.301232 61 80 3.0499999999999998 
0.324140 0.798707 62 81 3.0617283950617282 
0.217359 0.843892 63 82 3.0731707317073171 
0.886652 0.561197 63 83 3.0361445783132530 
0.328966 0.493458 64 84 3.0476190476190474 
0.934263 0.125679 65 85 3.0588235294117645 
0.326040 0.285874 66 86 3.0697674418604652 
0.664750 0.537354 67 87 3.0804597701149423 
0.156177 0.461966 68 88 3.0909090909090908 
0.451539 0.913815 68 89 3.0561797752808988 
0.337575 0.602196 69 90 3.0666666666666669 
0.767945 0.087343 70 91 3.0769230769230771 
0.658919 0.797220 70 92 3.0434782608695654 
0.223290 0.182301 71 93 3.0537634408602150 
0.415846 0.260215 72 94 3.0638297872340425 
0.219957 0.965141 73 95 3.0736842105263156 
0.561447 0.544097 74 96 3.0833333333333335 
0.763848 0.778807 74 97 3.0515463917525771 
0.387989 0.650500 75 98 3.0612244897959182 
0.340004 0.716955 76 99 3.0707070707070705 

点の数が100程度では、3.14まで出ない。

100000程度では3.14までは出てくる。

0.955767 0.364354 78590 99990 3.1439143914391439 
0.527298 0.749567 78591 99991 3.1439229530657760 
0.066892 0.856223 78592 99992 3.1439315145211615 
0.640014 0.923527 78592 99993 3.1439000730051103 
0.868690 0.540803 78592 99994 3.1438686321179272 
0.184084 0.428770 78593 99995 3.1438771938596930 
0.094987 0.893457 78594 99996 3.1438857554302171 
0.950835 0.397854 78594 99997 3.1438543156294689 
0.543792 0.468446 78595 99998 3.1438628772575452 
0.559280 0.790408 78596 99999 3.1438714387143873 

さらに10 000 000(一千万)程度では3.1415に近くなるように見えるが、かえって離れていくこともあり、これ以上は多く計算するほど正確になるというものではないようだ。

0.925574 0.527996 7854056 9999988 3.1416261699514041 
0.675092 0.279533 7854057 9999989 3.1416262557888812 
0.810590 0.219175 7854058 9999990 3.1416263416263415 
0.438512 0.324193 7854059 9999991 3.1416264274637848 
0.570633 0.392716 7854060 9999992 3.1416265133012105 
0.144585 0.418910 7854061 9999993 3.1416265991386192 
0.700746 0.257210 7854062 9999994 3.1416266849760111 
0.624001 0.169051 7854063 9999995 3.1416267708133856 
0.367598 0.677451 7854064 9999996 3.1416268566507428 
0.658572 0.160272 7854065 9999997 3.1416269424880827 
0.013899 0.285122 7854066 9999998 3.1416270283254057 
0.216488 0.829218 7854067 9999999 3.1416271141627115 

コンピュータの浮動小数点数は有限桁なので不連続な値です。本物の平面は連続な値でてきている(つまりどんなに近い点でも2つの点の間にまだ点を取ることができる)ので正確には表現できないのです。

また、コンピュータの出す乱数の品質の問題もあります。ほんとうにデタラメかという問題です。コンピュータの出す乱数は擬似乱数と呼ばれます。