モンテカルロ法 円周率算出

3.141まで収束するのに9分30秒掛かりました。一回計算の度に表示しているので時間掛かります。

1. 0～1.0までの間の浮動小数のランダム値を二つ得て、1.0 x 1.0の正方形内に仮想的にプロット。
2. xx+yy=1の方程式により、xx+yy<=1なら円内として c++。
3. *4してそれを計算した回数(N)で割ると円周率の近似値になります。
   近似解π=4C/N。Nが多いほど円周率に収束し有効桁数が増えます。

Output()関数は毎回、除算と表示を行っているので遅いですが、
収束具合が見えるので面白いです。

NoOutput()関数は速度優先で処理します。
速度優先と言っても1000億回で30分位です。1兆回では5時間は掛かります。

起動して1を選ぶと毎回表示で2では速度優先で処理します。
WIN32 APIと比べるため、VisualC++6.0版で作って比べた所、C#に比べてC言語では1/3の時間でした。

C# 10億回ループ 21sec (.NET FRAMEWORK)
C++10億回ループ 6sec (WIN32 API)
3.5倍。流石のWIN32ネイティブ。
(約33分 Core i5-7600K,C#)

回数(N)	pi≒ (約33分 Core i5-7600K,C#)
1000億	3.141592
100億	3.141600
10億	3.141569
1億	3.141713
1000万	3.140934
100万	3.146928
10万	3.136
1万	3.132

C# .NETFRAME WORK , Console App版

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++版

//参考 C++ 最近、VS2021で作り直したもの。なんか遅い
#include <iostream>
#include <time.h>
#include <random>

int main(){

    long cnt = 0;
    double X, Y;
    long c = 0;

    std::random_device rd;
    std::default_random_engine eng(rd());
    std::uniform_real_distribution<double> distr(0.1, 1.0);

    double g = 0;
    long f = 100000000;
    for (int n = 0; n < f; ++n) {

        X = distr(eng);
        Y = distr(eng);

        if ( (X * X + Y * Y) < 1.0) c++;

       // g= (4.0 * c / n);
       // std::cout << n <<" " << g << "\n";

        
    }

    g= (4.0 * c / f);
    std::cout << " " << g << "\n";
  
    int a;
    std::cin >> a;

}