2013/06/03

[C++]ビュフォンの針

,
ビュフォンの針のコード(C++)。課題でやらざるを得なかったので、メモ代わりに。

#include <iostream>
#include <cmath>
#include <random>
#include <algorithm>
#include <cstdint>

int main()
{
 const auto n = 10000000; //試行回数
 const auto d = 100.0; //平行線の間隔 //1.0
 const auto width = 100 * 255; //平行線の数 
 const auto l = 75; //針の長さ
 const auto M_PI = 6 * asin(0.5); //pi


 //乱数の種をまく
 std::random_device rnd;

 std::vector<std::uint_least32_t> v(10);

 std::generate(v.begin(), v.end(), std::ref(rnd));

 std::mt19937 engine(std::seed_seq(v.begin(), v.end()));

 std::uniform_real_distribution<double> distribution1(0.0, width);
 std::uniform_real_distribution<double> distribution2(0.0, 90.0);
/* std::uniform_real_distribution<double> dist(0.0, 1.0);

 //割ると(-1か1)を返す
 auto rev = [&dist, &engine] (double src) -> double {
  (dist(engine) < 0.5) ? src *= 1.0 : src *= -1.0;
  return src;
 };
*/
 //xの手前で一番近い平行線との距離
 auto func = [&d] (double x) -> double {
  while(x > (2.0 * d)){
   x -= 2.0 * d;
  }
  return x;
 };
 
 //ラジアンに直す
 auto radian = [&M_PI](double degree) -> double {
  return (degree / 180.0) * M_PI;
 };
 
 auto count = 0;

 //ビュフォンの針
 for(auto i = 0; i < n; i++){

  //針の中心位置
  auto x = distribution1(engine);

  //向いてる角度(-90度~90度)
  auto deg = distribution2(engine);
  
  //ラジアンに直す
  auto rad = radian(deg);

  if((func(x) + l * std::sin(rad)) >= (2.0 * d) || (func(x) - l * std::sin(rad)) <= 0.0){
   //円内にあるよ
   count++;
  }
 }

 //円周率を求める
 auto pi = 2.0 * (2.0 * l) * n / (2.0 * d * count);

 //出力する
 std::cout << "pi = " << pi << std::endl;

 return 0;
}
一応求まるけど収束性?がとても悪いです。

 一番困ったのは、最初定数を

0 コメント to “[C++]ビュフォンの針”

コメントを投稿