モンテカルロ法による円周率算出 (2020年10月)
0~1.0までの間の浮動小数のランダム値を二つ得て、1.0 x 1.0の正方形内に仮想的にプロット。
x*x+y*y=1の方程式により、x*x+y*y<=1なら円内としてカウント(C)。 *4してそれを回数(N)で割ると円周率の近似値になります。 近似解π=4C/N。Nが多いほど円周率に収束し有効桁数が増えます。 Output()関数は毎回、除算と表示を行っているので遅いですが、 収束具合が見えるので面白いです。 NoOutput()関数は速度優先で処理します。 速度優先と言っても1000億回で30分位です。1兆回では5時間は掛かります。 起動して1を選ぶと毎回表示で2では速度優先で処理します。 WIN32 APIと比べるため、VisualC++6.0版で作って比べた所、C#に比べてC言語では1/3の時間でした。 作るのは楽ですが.NETフレームワークは明らかに速度面で劣ります。 C# 10億回ループ 21sec (.NET FRAMEWORK) C++10億回ループ 6sec (WIN32 API) 3.5倍。流石のWIN32ネイティブ。 1000億回ループしたデータです。 3.1415926535..... ------------------------------ 1000億回 pi ≒ 3.141592 (約33分 Core i5-7600K,C#) 100億回 pi ≒ 3.141600 10億回 pi ≒ 3.141569 1億回 pi ≒ 3.141713 1000万回 pi ≒ 3.140934 100万回 pi ≒ 3.146928 10万回 pi ≒ 3.136 1万回 pi ≒ 3.132
///////////////////////////////////
/// Simulation MonteCarloMethod.
/// By S.H VisualStdio Win-Forms
///////////////////////////////////
using System;
namespace ConsoleApp2 {
class Program {
static void Main(string[] args) {
Random Ra = new System.Random();
Console.WriteLine("Simulation MonteCarloMethod.\r\nMode Select Output...1 , NoOutput...2");
int mode = int.Parse(Console.ReadLine());
if (mode == 2) {
NoOut();
}
else Output();
}
static void Output() {
Random Ra = new System.Random();
ulong cnt = 0;
double X, Y;
ulong c = 0;
while (true) {
X = Ra.NextDouble();
Y = Ra.NextDouble();
if (X * X + Y * Y < 1.0) c++;
Console.WriteLine(cnt + " " + 4.0 * c / cnt + "\r\n\r\n\r\n");
cnt++;
}
}
static void NoOut() {
Random Ra = new System.Random();
ulong max = 10L;
ulong cnt = 0;
ulong c = 0;
double X, Y;
DateTime dt = DateTime.Now;
while (true) {
dt = DateTime.Now;
Console.WriteLine("StartTime " + dt);
Console.WriteLine("Loop " + max);
while (true) {
X = Ra.NextDouble();
Y = Ra.NextDouble();
if (X * X + Y * Y < 1.0) c++;
if (cnt == max) {
max *= 10;
Console.WriteLine(4.0 * c / cnt + "\r\n");
break;
}
cnt++;
}
} // MAINLOOP
}
}
}
// C++
#include "stdafx.h"
#include <time.h>
#include <stdlib.h>
#include <stdio.h>
int main(int argc, char* argv[]){
double x,y;
double c;
double long in=0L;
char tmp[1];
double max=10L;
double cnt2=0L;
srand((unsigned int)time(NULL));
time_t timer = time(NULL);
struct tm *date = localtime(&timer);
printf("Start ... %d.%02d.%02d %2d:%02d:%02d\r\n", date->tm_year+1900, date->tm_mon+1, date->tm_mday, date->tm_hour, date->tm_min, date->tm_sec);
while (1){
printf ("LOOP %f\r\n",max);
while (1){
x= rand() / (double)RAND_MAX;
y= rand() / (double)RAND_MAX;
if ( (x*x + y*y) < 1.0) in++;
cnt2++;
if (cnt2 == max) break;
}
printf("%lf\r\n\r\n",4.0*in/cnt2);
max*=10;
timer = time(NULL);
date = localtime(&timer);
printf("Start ..... %d.%02d.%02d %2d:%02d:%02d\r\n", date->tm_year+1900, date->tm_mon+1, date->tm_mday, date->tm_hour, date->tm_min, date->tm_sec);
}
return 0;
}