モンテカルロ法による円周率算出 (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;
}