2014年6月9日月曜日

C言語で一様乱数を作ろうとしてうっかりする

C言語で乱数を発生させるときによく使う「0以上1未満の一様乱数を作る」方法なのですが、適当にインターネットからコピペして使っていると、編集ミスなどをやらかしヒドイ目に遭うようです。型変換を理解せず、コンパイラの警告も適当に流して使ったりすると、そういうことになるので、学生さんは気をつけましょう、という話。

 まずは12個、パターンを列挙してみます。

double randnnA(){
 return rand() / (RAND_MAX + 1);
}

double randnnB(){
 return rand() / (RAND_MAX + 1.0);
}

double randndA(){
 return rand() / (double)(RAND_MAX + 1);
}

double randndB(){
 return rand() / (double)(RAND_MAX + 1.0);
}

double randndC(){
 return rand() / ((double)RAND_MAX + 1);
}

double randndD(){
 return rand() / ((double)RAND_MAX + 1.0);
}

double randdnA(){
 return (double)rand() / (RAND_MAX + 1);
}

double randdnB(){
 return (double)rand() / (RAND_MAX + 1.0);
}

double randddA(){
 return (double)rand() / (double)(RAND_MAX + 1);
}

double randddB(){
 return (double)rand() / (double)(RAND_MAX + 1.0);
}

double randddC(){
 return (double)rand() / ((double)RAND_MAX + 1);
}

double randddD(){
 return (double)rand() / ((double)RAND_MAX + 1.0);
}
コンパイル可能なソースコードはこちら。
https://gist.github.com/taraijpn/7789f5e5382a2111a6c4
(動作を確認するため若干いじってあります。)

おかしな動きをする関数が4つあるのですが、分かりますでしょうか? 

A で終わっている関数は全て挙動が不適切になります。

$ ./a.exe
set 0 and RAND_MAX
randnnA min:0    max:0
randnnB min:0    max:0.9999999995343387
randndA min:-0   max:-0.9999999995343387
randndB min:0    max:0.9999999995343387
randndC min:0    max:0.9999999995343387
randndD min:0    max:0.9999999995343387
randdnA min:-0   max:-0.9999999995343387
randdnB min:0    max:0.9999999995343387
randddA min:-0   max:-0.9999999995343387
randddB min:0    max:0.9999999995343387
randddC min:0    max:0.9999999995343387
randddD min:0    max:0.9999999995343387

generate random number 10000 times
randnnA min:0    max:0
randnnB min:0.0001787277869880199        max:0.9999335804022849
randndA min:-0.9998125517740846          max:-0.0001268410123884678
randndB min:6.128335371613503e-05        max:0.9999860003590584
randndC min:1.606764271855354e-05        max:0.99976617237553
randndD min:6.103422492742538e-06        max:0.9999661864712834
randdnA min:-0.9998871022835374          max:-3.686733543872833e-05
randdnB min:2.388842403888702e-07        max:0.9999870401807129
randddA min:-0.9999983948655427          max:-2.792663872241974e-05
randddB min:5.736947059631348e-05        max:0.9998592492192984
randddC min:3.99700365960598e-05         max:0.9999427762813866
randddD min:8.147209882736206e-05        max:0.999992452096194

randnnAがダメな理由、まあどちらもintだし、というところまでは想像つくと思いますが、何故(RAND_MAX)/(RAND_MAX+1)がdouble型に代入されると0.0になるのか、はちょっと説明に引っかかるかも? また、なにか演算した後にあわてて型変換しても手遅れ(randndAがダメ)なことも分かります。

それと、分子側が(double)で型変換されてれば問題ないだろ、と思ってしまうのもいけません(randdnA, randddAがダメ)です。 [0,1)の一様乱数を作ろうとしているのに負の値が出てる時点で何かおかしいと気づくだろうと思います。

CとDは分母の (double)RAND_MAX でdouble型になるので、これに1を足そうが1.0を足そうがいずれにせよdouble型であることは変わりません。書く意味はほとんどなかったかも。まあパターンで列挙した次第です。

gccのバージョンが新しめならコンパイラがWarningで教えてくれるようですので、 Warningは見落とさないようにしましょう。

$ gcc --version
gcc (GCC) 4.8.2
Copyright (C) 2013 Free Software Foundation, Inc.
This is free software; see the source for copying conditions.  There is NO
warranty; not even for MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
$ gcc urandtest.c
urandtest.c: In function ‘randnnA’:
urandtest.c:8:24: warning: integer overflow in expression [-Woverflow]
  return ir / (RAND_MAX + 1);
                        ^
urandtest.c: In function ‘randndA’:
urandtest.c:18:32: warning: integer overflow in expression [-Woverflow]
  return ir / (double)(RAND_MAX + 1);
                                ^
urandtest.c: In function ‘randdnA’:
urandtest.c:38:32: warning: integer overflow in expression [-Woverflow]
  return (double)ir / (RAND_MAX + 1);
                                ^
urandtest.c: In function ‘randddA’:
urandtest.c:48:40: warning: integer overflow in expression [-Woverflow]
  return (double)ir / (double)(RAND_MAX + 1);
                                        ^
4.7でもいけました。コンパイラのバージョンが古かったり32bit環境だったりすると、なんとなく正しく動くものが増えるところが厄介です。
個人的には、ちょっと使う時には Mersenne Twister with improved initialization (2002) をお勧めしています。 本格的に使う場合は改良版(SFMT)の導入を検討したほうがよいでしょう。

おまけ。

double lrandnnA(){
 return rand() / (RAND_MAX + 1L);
}

double lranddnA(){
 return (double)rand() / (RAND_MAX + 1L);
}

整数リテラルってやつですね。それぞれ希望した結果になるでしょうか? ならないでしょうか?

0 件のコメント: