質問をすることでしか得られない、回答やアドバイスがある。

15分調べてもわからないことは、質問しよう!

ただいまの
回答率

89.06%

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

解決済

回答 3

投稿

  • 評価
  • クリップ 0
  • VIEW 180

nanyakanya

score 16

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

y軸上に0~1の範囲をとって、そこにながさ1の針の始点yをランダムでとり、またθを針とx軸とのなす角として、θを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;
   }


よろしくお願いします。

  • 気になる質問をクリップする

    クリップした質問は、後からいつでもマイページで確認できます。

    またクリップした質問に回答があった際、通知やメールを受け取ることができます。

    クリップを取り消します

  • 良い質問の評価を上げる

    以下のような質問は評価を上げましょう

    • 質問内容が明確
    • 自分も答えを知りたい
    • 質問者以外のユーザにも役立つ

    評価が高い質問は、TOPページの「注目」タブのフィードに表示されやすくなります。

    質問の評価を上げたことを取り消します

  • 評価を下げられる数の上限に達しました

    評価を下げることができません

    • 1日5回まで評価を下げられます
    • 1日に1ユーザに対して2回まで評価を下げられます

    質問の評価を下げる

    teratailでは下記のような質問を「具体的に困っていることがない質問」、「サイトポリシーに違反する質問」と定義し、推奨していません。

    • プログラミングに関係のない質問
    • やってほしいことだけを記載した丸投げの質問
    • 問題・課題が含まれていない質問
    • 意図的に内容が抹消された質問
    • 過去に投稿した質問と同じ内容の質問
    • 広告と受け取られるような投稿

    評価が下がると、TOPページの「アクティブ」「注目」タブのフィードに表示されにくくなります。

    質問の評価を下げたことを取り消します

    この機能は開放されていません

    評価を下げる条件を満たしてません

    評価を下げる理由を選択してください

    詳細な説明はこちら

    上記に当てはまらず、質問内容が明確になっていない質問には「情報の追加・修正依頼」機能からコメントをしてください。

    質問の評価を下げる機能の利用条件

    この機能を利用するためには、以下の事項を行う必要があります。

質問への追記・修正、ベストアンサー選択の依頼

  • ikadzuchi

    2020/07/31 09:55

    再現しませんでした。
    当該コードを私の環境で実行したところ、一例として以下になりました。
    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などの結果はどうやって出たものですか?

    キャンセル

  • nanyakanya

    2020/07/31 10:08

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

    キャンセル

回答 3

+1

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


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

投稿

編集

  • 回答の評価を上げる

    以下のような回答は評価を上げましょう

    • 正しい回答
    • わかりやすい回答
    • ためになる回答

    評価が高い回答ほどページの上位に表示されます。

  • 回答の評価を下げる

    下記のような回答は推奨されていません。

    • 間違っている回答
    • 質問の回答になっていない投稿
    • スパムや攻撃的な表現を用いた投稿

    評価を下げる際はその理由を明確に伝え、適切な回答に修正してもらいましょう。

  • 2020/08/01 17:43

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

    キャンセル

  • 2020/08/01 17:51

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

    キャンセル

  • 2020/08/01 17:56

    確かに
    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
    似たようなオーダーになりますね

    キャンセル

checkベストアンサー

0

気になったところ

  • 配列は必要ありません。冗長です。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;
}

投稿

編集

  • 回答の評価を上げる

    以下のような回答は評価を上げましょう

    • 正しい回答
    • わかりやすい回答
    • ためになる回答

    評価が高い回答ほどページの上位に表示されます。

  • 回答の評価を下げる

    下記のような回答は推奨されていません。

    • 間違っている回答
    • 質問の回答になっていない投稿
    • スパムや攻撃的な表現を用いた投稿

    評価を下げる際はその理由を明確に伝え、適切な回答に修正してもらいましょう。

  • 2020/07/31 09:15

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

    キャンセル

0

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

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

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

投稿

編集

  • 回答の評価を上げる

    以下のような回答は評価を上げましょう

    • 正しい回答
    • わかりやすい回答
    • ためになる回答

    評価が高い回答ほどページの上位に表示されます。

  • 回答の評価を下げる

    下記のような回答は推奨されていません。

    • 間違っている回答
    • 質問の回答になっていない投稿
    • スパムや攻撃的な表現を用いた投稿

    評価を下げる際はその理由を明確に伝え、適切な回答に修正してもらいましょう。

  • 2020/07/31 09:15

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

    キャンセル

15分調べてもわからないことは、teratailで質問しよう!

  • ただいまの回答率 89.06%
  • 質問をまとめることで、思考を整理して素早く解決
  • テンプレート機能で、簡単に質問をまとめられる

関連した質問

同じタグがついた質問を見る