高精度乱数関数

質問

rand()の値をRAND_MAX+1で割ると、[0, 1)の範囲になりますよね。
もし、RAND_MAX+1が2の乗数なら、浮動小数において結果は丸め誤差を含んでいないはずですよね。
unixのgccでコンパイルしたのですが、RAND_MAXの値は32767だったので、

#define (RM RAND_MAX + 1)
#define RND ((long double)rand())
long double x;
x = RND / RM
+ RND / RM / RM
+ RND / RM / RM / RM
+ RND / RM / RM / RM / RM
・
・
・;

とすれば、精度の高い[0, 1)区間の乱数関数を作ることはできますか?
この論理に間違いがあれば指摘してください。

回答No.9

Fooky · ベストアンサー

もうここまで来たら、「しつこいよう」
じゃなくて明らかにしつこいですが、
答えは変わらないにしても間違いがあったので…。

> ですが、今回はlog_2(M) / log_2(32768)が取りうる
> 最大の整数を求めたいので、上限値だけ考えればよく、

というのはうそですね。

丸めを防ぐためには安全側に考えて、Ｍの下限値の方を
計算すべきでした。

Mをlong doubleの仮数最大値＋１とすると、
10^34 < M <= 10^35　　（細かい間違いですが、１足してるので
34 < log_10(M) <= 35　　等号は右辺につきます）
112.9 < log_2(M) <= 116.3

Mは仮数部の全ビットが立ったものに１を足しているので
２の乗数です（正確には仮数部の全ビットが立ったものに
正規化によるビット増の１ビットを付け加えたものに
１を足したものがM）。したがってlog_2(M)は整数になります。
log_2(M) = {113, 114, 115, 116}

これは、(M-1)のビット数は113ビット以上116ビット
以下であるということに相当します。

したがって、丸めが起こらないようにするには安全側に
考えて、(M-1)は113ビットの数であると考えるべきです。

32767は15ビットの数なので、
113 / 15 = 7.53

したがって、可能な繰り返し回数は最大の整数を取って、
（これと混同したわけです）７である、というのが
冗長ですが正しい計算過程です。

今回は上限をとっても下限をとっても同じ答えですが、
たとえば(M-1)のビット数の上限が134とかになった場合、
８回繰り返せる可能性がかなり高いわけですが、
涙を飲んで下限のビット数113を取って７回と
すべきである訳です。

お礼

だんだん難しい話になってきましたね。
コンピュータと言うよりは、数学の話みたいです。
これ以上はもうちょっと時間をかけて考えなければならなそうです。
自分でもうちょっと考えてみます。
ありがとうございました。

回答No.8

Fooky

どうもしつこいようですが…。

nakashiさんのあげられてるようなページを見て、
一般的なアルゴリズムで乱数発生器をつくるのが、
最もお勧めだと私も思います。

> long doubleが10進で35桁、32768が5桁だから、
> 35/5=7と出してしまったのですが、内部では2進なので、
> 実際には何次くらいまで計算するべきかを考えると、
> 方法がわかりません。

についてですが、long doubleの最大値＋１を
10進数Ｍとすると、可能な繰り返し回数は、
log_2(M) / log_2(32768) で表されます。

Ｍは１０進数で３５桁なので、
34 <= log_10(M) < 35

ですが、今回はlog_2(M) / log_2(32768)が取りうる
最大の整数を求めたいので、上限値だけ考えればよく、
log_10(M) < 35

log_p(q) / log_p(r) = log_r(q)より、
両辺をlog_10(2)で割ると、2進数の桁数に直り、
log_2(M) < 35 / log_10(2)
log_2(M) < 116.3

一方、log_2(32768)は１５。
よって、両辺を１５で割ると、
log_2(M) / log_2(32768) < 7.75

となりますので、最大の整数を取ると７。
よって、long doubleが本当に35桁であるなら、
７回で正解だと思います。

お礼

そうか!
対数をとればよかったのか。
2進と10進の変換については、あまり頭が回らなくて。
どうもありがとうございました。

回答No.7

nakashi

乱数に関しては、下記のURLがお勧めです。

参考URL：

http://www-6.ibm.com/jp/developerworks/security/000818/j_playing.html

お礼

詳しいページのご紹介どうもありがとうございます。

回答No.6

Fooky

えー、通してみてみると、大体、３者の認識は
合っているように思います。私が杞憂に過ぎた
のかもしれません。

「２の乗数で割れば丸めがない」に対する疑問は、
単に、指数部の限界を超えれば丸めが起こるんだから、
一般的にそういうことは言えないですよね、という
ことです。
（指数部のビット数をよく調べずに書いてたんですが、
Linuxではdoubleでも10ビットあるので、1000桁分以上
あるわけで、気にすることもないんでしょうね）

さいころの話は、Σ1/(RAND_MAX^n)の計算においては、
仮数部の有効桁数が問題になってくる、ということを
言ってるつもりなんですが、たとえが悪かったかもしれません。

ｎ回目の乱数発生と除算の結果をf(n)と記述するとし、
（ややこしいので10進数で話しますが、）有効桁数が４桁で
あるとすると、例えば、

f(1) = 0.11
f(3) = 0.00001
Σf(n) = 1.100e-1

f(1) = 0.11
f(3) = 0.00002
Σf(n) = 1.100e-1

f(1) = 0.01
f(3) = 0.00001
Σf(n) = 1.001e-2

というように、仮数部で丸めが発生すると、ある X = Σf(n) に対して
可能なf(n)の系列が複数存在する可能性が発生する、そうすると、
さいころと同様に、一様性が崩れるということです。（f(1)の仮数部桁数が
大きいXになる確率が高くなります。）

> だから、実際は7次くらいで止めるつもりです。

と、仰ってるので大丈夫なんでしょう。その辺の上限をどう考えられて
るのかな、ということが気になっていました。

また、整数を浮動小数点に変換する処理と浮動小数点計算の内部処理が、
どの程度理論的に期待される処理と合致するのか（変換過程や計算過程で
ゴミが入ったりしないのか）ということも気になっていましたが、
確かに、Σrand()/（RAND_MAX＾n） は期待する結果（各rand()の仮数部が
順に並んだ仮数部）が得られるんですね。手もとのLinuxでもそうなりました。

お礼

たびたびありがとうございます。
一様性が失われる部分について理解しました。
これなら、long doubleの有効桁数内でやる限りは、大丈夫でしょう。

ところで、簡単に7次と言いましたが、後で考えてみると不安になりました。
long doubleが10進で35桁、32768が5桁だから、35/5=7と出してしまったのですが、内部では2進なので、実際には何次くらいまで計算するべきかを考えると、方法がわかりません。
浮動小数ではごみを加減しても、ほとんど何も起こりませんが、ごみになる前に止めるのがセオリーでしょう。
10進←→2進は色々と複雑なので、難しいですね。

回答No.5

punchan_jp

えっと、試しにプログラム作りました。
main()
{
double a = 1;
double b = a/32768;
double c = b/32768;
double d = c/32768;
double e = d/32768;
double x = b + c + d + e;
printf("%08x:%08x\n", *((long*)&a+1), *((long*)&a+0));
printf("%08x:%08x\n", *((long*)&b+1), *((long*)&b+0));
printf("%08x:%08x\n", *((long*)&c+1), *((long*)&c+0));
printf("%08x:%08x\n", *((long*)&d+1), *((long*)&d+0));
printf("%08x:%08x\n", *((long*)&e+1), *((long*)&e+0));
printf("%08x:%08x\n", *((long*)&x+1), *((long*)&x+0));
}
を、FreeBSD の gcc でコンパイルし、実行すると
3ff00000:00000000
3f000000:00000000
3e100000:00000000
3d200000:00000000
3c300000:00000000
3f000020:00400080
となり、IEEE754の倍精度浮動小数点形式として、完全に誤差なく
表現されていることがわかります。a に 1 じゃなく、5 など別の
値を入れても同様です。

long double については試していません。

浮動小数点数についても、2の羃での除算は仮数部のシフト+正規化
でしかありませんから、
>0.<1回目のrand()のbitパターン><2回目のrand()のbitパターン>
になりそうであることは、期待できるわけです。

お礼

証明してくださってどうもありがとうございます。
ビット合成についても考えてみようと思うのですが、どちらがよいアルゴリズムと言えるでしょう?

回答No.4

Fooky

> 0.<1回目のrand()のbitパターン><2回目のrand()のbitパターン>
> としたいという意味です。（よね？）

その解釈は成り立たないでしょう。
その理由は、

１．除算前にrand()をlong doubleにキャストしている
　　　したがって、除算は４倍精度の浮動小数点形式で行われます。
　　　２の乗数による除算が単純な右シフトで表現されるデータ形式は、
　　　符合無し整数型だけです。

また、もし整数での除算であったとしても、

２．RAND_MAX+1での除算に相当する右シフトの結果はかならず０になる
　　　非除数は最大でもRAND_MAXなので、８ビット目より上位にビットが
　　　立つことはないんですから

> 0.<1回目のrand()のbitパターン><2回目のrand()のbitパターン>

このような意図であれば、bitパターンの生成は整数型で行ない、
また、除算ではなく、１回目のrand()にRAND_MAX+1を掛けて
左シフトし、2回目を足し合わせる。ということをした後に、
RAND_MAXの２乗で割る、という手順を踏まないと行けないでしょう。

long double x;
int i;
unsigned long k;

for( i = 0, k = 0 ; i < 4 ; i++ ){
　k = k << 15;
　k += (unsigned long)random() & 0x7fff;
}
x = (long double)k;
for( i = 0 ; i < 4 ; i++ )
　x /= (long double)(RAND_MAX + 1);

random()の出力が７ビット以上、longが８バイト
と仮定しています。

お礼

右シフトと言うよりは、指数部分の減算を期待しているのですが、これは間違いなのでしょうか。

確かに、ビット合成すれば正確な値は得られそうですが、それにはlong doubleの内部表現を完全に理解し、ビットフィールドを用いなければ、ならないような気がするのですが。

回答No.3

punchan_jp

えっと、Fooky さんにちょっとした誤解があるような気がしますの
で、追加です。質問者の意図は、2進数で考えたとき、
0.<1回目のrand()のbitパターン><2回目のrand()のbitパターン>
としたいという意味です。（よね？）

10進数で考えれば、0x8000 とか 0x80000000 とかで割ると、非常
に桁数が増えて丸めが必要になるような気がしますが、内部表現は
2進数ですから、単に 15bit あるいは 31bit シフトされるだけで
す（というか、それを期待してのアルゴリズムです）。もちろん、
この内部表現を10進に変換して表示するときには丸めて表示されま
すが、計算途上では丸めは起きません。

あと、もし rand() が理想的な一様乱数であれば、それぞれ重なら
ないようにシフトしたものを加算しても、何も問題なく、結果は一
様です。サイコロの問題は同じ桁で加算しているから起きるのです。

お礼

補足ありがとうございます。
問題なのはrand()の一様性だけなのでしょうか。

回答No.2

Fooky

分母が「2の乗数なら、浮動小数において結果は
丸め誤差を含んでいない」というのはなぜでしょう？
RAND_MAXが本当に32767であれば、1/(RAND_MAX+1)
は0.000030517578125となって、桁数は１１桁で
4倍精度なら確かに丸めは「たまたま」発生しません。

しかし、もしRAND_MAX=0x7fffffff（2147483647）
であれば、（現に私のところのLinuxはそうですが）
1/(RAND_MAX+1) =
0.0000000004656612873077392578125
で、桁数は22桁。4倍精度の有効桁数を正確には
把握してませんが、これもたまたまＯＫでしょうか？

２の乗数なら丸めが発生しない、なんて一般的には
言えなくて、有効桁数を上回ってしまえば当然、
丸めは発生するわけです。

では、1/(RAND_MAX+1)が４倍精度で収容できれば、
この手法はＯＫなのか、というとそうではなく、
さらに問題があります。

それは、結果の乱数ｘを1/(RAND_MAX+1)を公比とする
等比級数で表現している点です。これをやると、
桁数がどんどん増えるので、必ずどこかで丸めが
発生しますよね。

丸めが発生すると非常にまずい、というのは
ｘの計算の仕方から見てhaporunさんは認識されて
いるんだと思いますが、一応指摘しておくと、
一様性が崩れてしまいます。

さいころを２つ振ったときの和を考えてみれば
すぐに分かりますが、一様乱数の和は一様乱数には
なりませんもんね。

２つのさいころの和が２と１２になる確率は
それぞれ１／３６、３と１１は１／１８、というように、
足してその数になる２数の組合わせの数に
比例した確率になります。

ですから、丸めが発生すると、可能なrand()出力値の
組合わせが複数存在するようなｘが存在する可能性が
できるため、一様性が崩れてしまうわけです。

１つのrandom()で出力される値は、RAND_MAXで
割ろうが割るまいが、高々RAND_MAX+1通りしか
ないんで、複数の出力値を組合わせるという考えは、
より高精度の乱数を（アルゴリズムまで
下らずに）手軽につくる方法であるとは思いますが、
出力を足したり掛けたりするのは非常にまずく、
punchan_jpの仰るように、bit列を合成するなり、
桁ごとに0～9の乱数を振るなりした方が良いです。

そのときに、どうしてもrand()を使うなら、
rand() & 1とか、rand() % 10といった乱数の
作り方は止めましょう。理由は多分rand()の
manページに載ってますが、rand()は下位ビットの
乱雑性が低いらしいです。random()では
改善されていると思います。

お礼

指摘ありがとうございます。

最大で32767という数値は、10進数で5桁、2進数で15桁なので、long doubleの10進でおよそ35桁ある精度を満たしてみたいと思ったのです。
だから、実際は7次くらいで止めるつもりです。

乱数の値が0～0x7fffffffのように、32ビットフルに使っているとしても、0x80000000で割ることは指数部分を減算することでしかないので、精度は失われないと思ったのですが。

>さいころを２つ振ったときの和を考えてみれば～
これは場合が違います。
10面体のさいころで、1つを10の位、1つを1の位とすれば、0～99の整数を一様に出すことができます。
この場合は(RAND_MAX+1)面体というものを考えています。

やはり、random()はMSDNホームページで見つけられませんでした。
これは、unix固有の関数なのでしょうか。

回答No.1

punchan_jp

unix の RAND_MAX って、そんなに小さかったでしたっけ？手元の
FreeBSD では 0x7fffffff ですけど。

で、rand() がほんとうに乱数であるなら、その方法でいいように
見えます。でも2の羃の除算が、丸め誤差なくちゃんとシフトにな
るのか心配なので、安全のために、浮動少数点形式を調べて、
[0,1) となるビットパターンを合成してしまう方が気が楽ですけど。

ただ、rand() の呼び出し前後の独立性が弱いので、複数回続けて
rand() を呼んで数値の組を作ったり、それらをまとめて一つの数
値を得るというのは、乱雑性の点から問題があるでしょう。

あと、rand() は互換性のためにのみあるといってもよく、
random() を使うことが推奨されています。random() でも上の問題
は完全には消えませんが。

結局、浮動小数点数の最小の単位まで本当に十分乱雑な乱数を得た
いなら、必要な bit 数をちゃんと合成するルーチンを作る方がい
いでしょう。

お礼

詳しい説明ありがとうございます。

SunOS 5.6のようですが、RAND_MAXは32767でした。
まぁ、何にせよ、RAND_MAX+1が2のn乗であることは間違いないので。

UNIXの浮動小数が2進数で管理されているなら、2のn乗で割り切れないはずがないと思うので、丸め誤差が発生しないと思ったのですが。
10進数で 3257 / 10000 = 0.3257
2進数で 11 0111 / 100 0000 = 0.1101 11
のように、内部表現で必ず割り切れているはずです。
16進数で管理されているとしても、割り切れていると思うのですが。

>rand() の呼び出し前後の独立性が弱い
というのは、擬似乱数だからですか?

random()は、MSDNライブラリのサイトで見つけられませんでした。
今は手元にローカルライブラリがないので。