ビュッホンの針を用いたpiの近似

Question

ビュッフォンの針を用いてpiを導くプログラムを作成しようと思います。

y軸上に0～1の範囲をとって、そこにながさ1の針の始点yをランダムでとり、またθを針とｘ軸とのなす角として、θを0～2piまでの乱数とする。ここで針の終点がy+sinθであるから、その値が0以下か1以上であるならばnumを+1して最終的にm(1000)階の施行からnum/1000として確率を求める
悩んでいること

result1 :2.100000e-002
result2 :4.300000e-002
result3 :7.300000e-002
result4 :9.200000e-002
result5 :1.220000e-001
以上のような結果となってしまった。
あまりにも理論値とかけ離れているのがどこが原因なのかわからない
コード

#include<stdio.h>
#include<math.h>
#include<stdlib.h>
#include<time.h>
#define pi 3.141592
 int main(void){
   int m = 1000;
   int i,j,k,num = 0;
   srand ((unsigned int)time(NULL));
   int randam[5][1000] = {0};
   for(i = 0; i < 5 ; i++ ){
     for(j = 0; j < m; j++){
       randam[i][j] =rand();
     }
   }

   for(i = 0; i < 5; i++ ){
     for(j = 0 ; j < m ; j++ ){
     double x,y;

     x = (float)randam[i][j] / RAND_MAX;//0～1の範囲の小数点乱数
     y = x + sin( ((float)randam[i][m-j])* 2 * pi);

     if(y <= 0 || 1 <= y){
       num++;
     }

     //printf("%lf     %lf\n",x,y);
   }
    printf("result%d :%e\n",i+1,(double)num/m);
     }
     return 0;
   }

よろしくお願いします。

Comments

[ikadzuchi]

再現しませんでした。
当該コードを私の環境で実行したところ、一例として以下になりました。
result1 :6.430000e-01
result2 :1.301000e+00
result3 :1.925000e+00
result4 :2.553000e+00
result5 :3.224000e+00
初期化忘れでresult2以降結果が誤っている点を除き、これは想定通りの値です。(線の間隔の違いで)2で割り、逆数をとればπになります。
お示しの2.1e-2などの結果はどうやって出たものですか?

[naka]

atomを用いてwindowsのコマンドプロンプトでの結果となります。
再現性がないのはなんでなのかわからりません。
いまreple.itで試したところikadzuchiさんのようになりました。

Answer 1

ちょっと気づいたのですが、sin()の引数をラジアンでなく度だと解釈して計算すると、お示しのものによく似た出力が出ます。
考えにくいことではありますが、まさか引数を度で解釈する特殊なsin関数が用意されているということはないでしょうか。
atomは使ったことが無いのですが(そもそもatomはテキストエディタであってコンパイラは別にあるのではないかと思いますが)、その環境でsinに様々な値を入れるとどのような結果が出ますか?

---

(追記)
上記、表計算で値をいじって確かめていた際の気づきなのですが、この際コードを読み違えており、実際のコードでは以上のような結果は出ません。すみませんでした。(詳細はコメント参照)
ただ、そのような出力が出る理由が分からないことには変わりありません。
sinと、あとrandの出力もどのようになるか気になりますので試していただけると幸いです。

Comments

[Penpen7]

度だと解釈するとというのは、質問のコードの一文を
y = x + sin( ((float)randam[i][m-j])* 2 * 180)
という風にするという意味ですかね？
result1 :6.380000e-01
result2 :1.292000e+00
result3 :1.923000e+00
result4 :2.562000e+00
result5 :3.210000e+00
となってオーダーが一桁違うので、同じとは言えないとおもいます。

[ikadzuchi]

それは「『引数を度と解釈するsin関数』があったとして、それに正しい結果を出させるための入力」ですね。
引数をラジアンと解釈する(普通の)sin関数にその入力を入れてもほぼ正常な出力が出て参考になりません。
「引数を度と解釈するsin関数」を普通のsin関数で再現したいなら、
sin( ((float)randam[i][m-j])* 2 * pi * pi/180)
です。

[Penpen7]

なるほど。引数を度ととるために、現在の引数をラジアンに変更するという意味ですね。
しかし、y = x + sin( ((float)randam[i][m-j])* 2 * pi * pi/180);としても、
result1 :6.370000e-01
result2 :1.260000e+00
result3 :1.893000e+00
result4 :2.531000e+00
result5 :3.156000e+00
となりますが...

[ikadzuchi]

おや、本当ですね。
計算は表計算で行っていたのですが何か勘違いしていたかもしれません。すいません、考え直します。

[Penpen7]

承知しました
返信ありがとうございます

[ikadzuchi]

分かりました。sin関数に入れている乱数は0～1に正規化していない非常に大きい数なのを見落としていました。
これではsinの引数が度でもラジアンでも実質的に何の違いも出ません。

[Penpen7]

確かに
y = x + sin( ((float)randam[i][m-j])/RAND_MAX* 2 * pi * pi/180);
とすれば、
result1 :5.600000e-02
result2 :1.050000e-01
result3 :1.580000e-01
result4 :2.120000e-01
result5 :2.650000e-01
似たようなオーダーになりますね

Answer 2

気になったところ

• 配列は必要ありません。冗長です。0-1の値さえ出てくればいいので関数を作りました。
• 棒の長さなどが与えられていないです。(こちらで陽に与えました)
• math.hをインクルードしているのであれば、円周率は自分で定義せずM_PIを使いましょう。
• 答えの分数は逆数にしなければなりません。
• numは試行するたびに0に初期化しなければなりません。

コード

#define _USE_MATH_DEFINES
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
// 0から1の一様分布に従う乱数を返す
double calculateUniformRandom()
{
  return rand() / (double)(RAND_MAX + 1.0);
}

int main(void)
{
  // 棒の長さ
  const double l = 1.0;
  // 平行線の距離
  const double d = 2 * l;
  // 試行回数
  const int m = 1000;

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

  for (int i = 0; i < 5; i++) {
    int num = 0;
    for (int j = 0; j < m; j++) {
      double x = calculateUniformRandom() * 0.5 * d; // [0:d/2]の乱数
      double y = l * 0.5 * sin(calculateUniformRandom() * M_PI * 0.5);

      if (x - y <= 0 || x + y >= d) {
        num++;
      }
    }
    printf("result%d :%lf\n", i + 1, m / (double)num);
  }
  return 0;
}

結果

result1 :3.067485
result2 :3.215434
result3 :3.154574
result4 :3.215434
result5 :3.012048

円周率を使用しない例

#define _USE_MATH_DEFINES
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>

// [min, max]の範囲の一様分布に従う乱数を返す
double calculateUniformRandom(double min, double max)
{
  return min + (double)(rand() * (max - min + 1.0) / (1.0 + RAND_MAX));
}

// 円周率を使用せず大きさrのベクトルを計算する(一様分布)
void calculeteUniformRandomVectorNoPI(double* x, double* y, double r)
{
  double xtemp, ytemp, rtemp2;

  // 正方形の中から, 半径1の円内部の座標を得るまで繰り返す。
  do {
    xtemp = calculateUniformRandom(0, 1.0), ytemp = calculateUniformRandom(0, 1.0);
    rtemp2 = xtemp * xtemp + ytemp * ytemp;
  } while (rtemp2 > 1.0);

  // 正規化し, 大きさrにする。
  double rtemp = sqrt(rtemp2);
  *x = xtemp / rtemp * r;
  *y = xtemp / rtemp * r;
}

// 円周率を使用して大きさrのベクトルを計算(一様分布)
void calculeteUniformRandomVectorPI(double* x, double* y, double r)
{
  double theta = M_PI * 0.5 * calculateUniformRandom(0, 1.0);
  *x = r * cos(theta);
  *y = r * sin(theta);
}

double calculatePi(int tryNumber, double l)
{
  int num = 0;
  double d = l * 2;

  for (int i = 0; i < tryNumber; i++) {
    double centerY = calculateUniformRandom(0, d * 0.5);
    double directionX, directionY;

    calculeteUniformRandomVectorNoPI(&directionX, &directionY, l * 0.5);
    if (centerY - directionY <= 0.0 || centerY + directionY >= d)
      num++;
  }

  return tryNumber / (double)num;
}

int main(void)
{
  srand((unsigned int)time(NULL));

  for (int i = 0; i < 5; i++) {
    printf("result%d :%lf\n", i + 1, calculatePi(1000, 1.0));
  }
  return 0;
}

追記(最小限の修正)

ご質問のコードだと、平行線の間隔が1, 棒の長さが0.5になりますね。

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#define pi 3.141592
int main(void)
{
  int m = 1000;
  int i, j;
  srand((unsigned int)time(NULL));
  int randam[5][1000] = { 0 };
  for (i = 0; i < 5; i++) {
    for (j = 0; j < m; j++) {
      randam[i][j] = rand();
    }
  }

  for (i = 0; i < 5; i++) {
    int num = 0;
    for (j = 0; j < m; j++) {
      double x, y;
      x = (float)randam[i][j] / (RAND_MAX + 1.0); //0～1の範囲の小数点乱数
      y = x + 0.5 * sin(((float)randam[i][m - j]) / (RAND_MAX + 1.0) * 2 * pi);

      if (y <= 0 || 1 <= y) {
        num++;
      }
    }
    printf("result%d :%e\n", i + 1, m / (double)num);
  }
  return 0;
}

Comments

[naka]

ビュッフォンの針の理解自体が甘いことに気づきました。
別海まで示していただきありがとうございました！

Answer 3

そもそも円周率を求めたいのに

θを0～2piまでの乱数とする

と、前提に円周率を使っては意味がありません。

Comments

[naka]

ご指摘ありがとうございます！
改善します