乱数 (復習)
rand()関数を使うと疑似乱数を生成できる。
#include <stdlib.h>
void
srand(unsigned seed);
int
rand(void);
The rand() function computes a sequence of pseudo-random integers in the range of 0 to RAND_MAX (as defined by the header file <stdlib.h>).
The srand() function sets its argument seed as the seed for a new sequence of pseudo-random numbers to be returned by rand(). These sequences are repeatable by calling srand() with the same seed value.
If no seed value is provided, the functions are automatically seeded with a value of 1.
seedとして適当な値を選び、srand(seed)を呼ぶことで乱数の系列が決まる。
rand()を呼ぶと順次疑似乱数が返される。srand()を何度も呼ぶのは間違い。
http://www9.plala.or.jp/sgwr-t/lib/srand.html
課題1 モンテカルロ法により円周率を求めよ
正方形内部に点(x,y)をMAXCT個作り、半径1の扇の中に(x,y)が入ってる場合の数をctとして数えることで,円周率を求めることができる.
より具体的には0<x<1, 0<y<1の正方形内部に乱数を使って点(x,y)をMAXCT個作る。
その中で、半径1の扇型の中に(x,y)が入ってる場合の数をctとして数える。
扇型の面積はπ/4なので、ct/MAXCTがπ/4になるはず。
#include <stdio.h>
#include <stdlib.h>
#define MAXCT 30000
double random_float(double, double);
int main()
{
int seed;
int k,ct=0;
float x,y;
scanf("%d",&seed); /* シードの入力 */
printf("seed = %d\n",seed);
srand(seed); /* 乱数の初期化 */
// 0<x<1, 0<y<1の正方形内部に乱数を使って点(x,y)をMAXCT個作り、
// 半径1の扇型の中に(x,y)が入ってる場合の数をctとして数える操作をここに書く.
printf("Approx. PI = %f\n",4.0*ct/MAXCT);
} /* minとmaxの間の実数乱数を発生する関数 */
double random_float(double min, double max) {
return min + ((max-min)*rand())/RAND_MAX;
}
浮動小数点の計算
0.1を100回足し算しても10にはならない.
#include<stdio.h>
int main(){
float sum;
int i;
sum = 0;
for (i = 1; i <= 100; i++){
sum += 0.1;
}
printf("%f\n", sum);
return 0;
}
0.1 を二進数にすると0.00011001100...という循環小数になる.
これを回避するには
1. 有効数字を考慮してあきらめる
2. 一時的に整数を使って計算し,その結果を小数点数で表示する
(例:1を100回足して,最後に10で割る.)
課題2 パイこね変換
(カオスとは、数的な誤差の蓄積により予測できないような複雑な振る舞いのことを呼ぶ。ランダムにも見えるが、ランダムではなく決定論的に決まるものである。パイコネ変換はカオスを生み出す典型的な仕組みの一つである。)
- 「パイ」をこねるような変換操作を計算機で実行することを考える。上図に従って,パイこね変換を数式化する。パイの初期状態では,パイは0<=x<=1の領域に広がっている。
- 図には,これを二倍に引き延ばし,折り畳む様子が示してある。このパイこね操作を繰り返した時に,初期状態時の点 X の座標(つまりx)はどこに動くかを求めたい。当然,x の初期値は 0<=x<=1 を満たす任意の実数となる。点 x は「引き延ばし+折りたたみ」によって,どこに移動するかと言えば,x<=0.5 の場合は,2x となり,x>0.5 の場合は,2-2x となる。この 2x,あるいは,2-2x を再度初期値として「引き延ばし+折りたたみ」を延々と繰り返す演算をする。上記の演算を繰り返し,点 x がどう動くのか(どうパイがこねられるのか)を追跡します。
- 課題2-1:x=x0(xの初期値。引数として与える)から開始して,N回(これも,引数として与える)「引き延ばし+折りたたみ」操作を繰り返した後の x を戻す関数を作成しなさい。
- double pie( double x0, int N ); のような関数を作成しなさい。
- 課題2-2:x0 = 0.1 とします。N = 1, 2, 3,...,10 の時に,上記関数の返す値を求め,x がどのように動くのかを確認し,動作結果について説明せよ。N = 100 やN = 1000 の場合も求めよ。
- 何を書いていいかわからない場合は、「初期値鋭敏性」などでWebで調べてみると良い。