-PR-
解決
済み

高精度乱数関数

  • 暇なときにでも
  • 質問No.94094
  • 閲覧数756
  • ありがとう数12
  • 気になる数0
  • 回答数9
  • コメント数0

お礼率 77% (530/685)

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)区間の乱数関数を作ることはできますか?
この論理に間違いがあれば指摘してください。
通報する
  • 回答数9
  • 気になる
    質問をブックマークします。
    マイページでまとめて確認できます。

質問者が選んだベストアンサー

  • 回答No.9
レベル9

ベストアンサー率 71% (59/82)

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

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

というのはうそですね。

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

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

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

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

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

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

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

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

お礼率 77% (530/685)

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

その他の回答 (全8件)

  • 回答No.3
レベル11

ベストアンサー率 55% (155/280)

えっと、Fooky さんにちょっとした誤解があるような気がしますの で、追加です。質問者の意図は、2進数で考えたとき、 0.<1回目のrand()のbitパターン><2回目のrand()のbitパターン> としたいという意味です。(よね?) 10進数で考えれば、0x8000 とか 0x80000000 とかで割ると、非常 に桁数が増えて丸めが必要になるような気がします ...続きを読む
えっと、Fooky さんにちょっとした誤解があるような気がしますの
で、追加です。質問者の意図は、2進数で考えたとき、
0.<1回目のrand()のbitパターン><2回目のrand()のbitパターン>
としたいという意味です。(よね?)

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

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

お礼率 77% (530/685)

補足ありがとうございます。
問題なのはrand()の一様性だけなのでしょうか。
投稿日時 - 2001-06-23 11:00:59
  • 回答No.2
レベル9

ベストアンサー率 71% (59/82)

分母が「2の乗数なら、浮動小数において結果は 丸め誤差を含んでいない」というのはなぜでしょう? RAND_MAXが本当に32767であれば、1/(RAND_MAX+1) は0.000030517578125となって、桁数は11桁で 4倍精度なら確かに丸めは「たまたま」発生しません。 しかし、もしRAND_MAX=0x7fffffff(2147483647) であれば、(現に私のところのLi ...続きを読む
分母が「2の乗数なら、浮動小数において結果は
丸め誤差を含んでいない」というのはなぜでしょう?
RAND_MAXが本当に32767であれば、1/(RAND_MAX+1)
は0.000030517578125となって、桁数は11桁で
4倍精度なら確かに丸めは「たまたま」発生しません。

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

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

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

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

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

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

2つのさいころの和が2と12になる確率は
それぞれ1/36、3と11は1/18、というように、
足してその数になる2数の組合わせの数に
比例した確率になります。

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

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

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

お礼率 77% (530/685)

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

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

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

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

やはり、random()はMSDNホームページで見つけられませんでした。
これは、unix固有の関数なのでしょうか。
投稿日時 - 2001-06-23 10:58:54
  • 回答No.1
レベル11

ベストアンサー率 55% (155/280)

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

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

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

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

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

お礼率 77% (530/685)

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

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ライブラリのサイトで見つけられませんでした。
今は手元にローカルライブラリがないので。
投稿日時 - 2001-06-23 10:44:36
  • 回答No.4
レベル9

ベストアンサー率 71% (59/82)

> 0.<1回目のrand()のbitパターン><2回目のrand()のbitパターン> > としたいという意味です。(よね?) その解釈は成り立たないでしょう。 その理由は、 1.除算前にrand()をlong doubleにキャストしている    したがって、除算は4倍精度の浮動小数点形式で行われます。    2の乗数による除算が単純な右シフ ...続きを読む
> 0.<1回目のrand()のbitパターン><2回目のrand()のbitパターン>
> としたいという意味です。(よね?)

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

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

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

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

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

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

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()の出力が7ビット以上、longが8バイト
と仮定しています。
お礼コメント
haporun

お礼率 77% (530/685)

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

確かに、ビット合成すれば正確な値は得られそうですが、それにはlong doubleの内部表現を完全に理解し、ビットフィールドを用いなければ、ならないような気がするのですが。
投稿日時 - 2001-06-23 11:08:26
  • 回答No.5
レベル11

ベストアンサー率 55% (155/280)

えっと、試しにプログラム作りました。 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; ...続きを読む
えっと、試しにプログラム作りました。
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パターン>
になりそうであることは、期待できるわけです。
お礼コメント
haporun

お礼率 77% (530/685)

証明してくださってどうもありがとうございます。
ビット合成についても考えてみようと思うのですが、どちらがよいアルゴリズムと言えるでしょう?
投稿日時 - 2001-06-23 11:11:01
  • 回答No.7
レベル8

ベストアンサー率 51% (21/41)

乱数に関しては、下記のURLがお勧めです。 ...続きを読む
乱数に関しては、下記のURLがお勧めです。
お礼コメント
haporun

お礼率 77% (530/685)

詳しいページのご紹介どうもありがとうございます。
投稿日時 - 2001-06-26 12:47:41
  • 回答No.6
レベル9

ベストアンサー率 71% (59/82)

えー、通してみてみると、大体、3者の認識は 合っているように思います。私が杞憂に過ぎた のかもしれません。 「2の乗数で割れば丸めがない」に対する疑問は、 単に、指数部の限界を超えれば丸めが起こるんだから、 一般的にそういうことは言えないですよね、という ことです。 (指数部のビット数をよく調べずに書いてたんですが、 Linuxではdoubleでも10ビットあるので、1000桁分以上 ...続きを読む
えー、通してみてみると、大体、3者の認識は
合っているように思います。私が杞憂に過ぎた
のかもしれません。

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

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

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

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でもそうなりました。
お礼コメント
haporun

お礼率 77% (530/685)

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

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

ベストアンサー率 71% (59/82)

どうもしつこいようですが…。 nakashiさんのあげられてるようなページを見て、 一般的なアルゴリズムで乱数発生器をつくるのが、 最もお勧めだと私も思います。 > long doubleが10進で35桁、32768が5桁だから、 > 35/5=7と出してしまったのですが、内部では2進なので、 > 実際には何次くらいまで計算するべきかを考えると、 > 方法がわ ...続きを読む
どうもしつこいようですが…。

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

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

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

Mは10進数で35桁なので、
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)は15。
よって、両辺を15で割ると、
log_2(M) / log_2(32768) < 7.75

となりますので、最大の整数を取ると7。
よって、long doubleが本当に35桁であるなら、
7回で正解だと思います。
お礼コメント
haporun

お礼率 77% (530/685)

そうか!
対数をとればよかったのか。
2進と10進の変換については、あまり頭が回らなくて。
どうもありがとうございました。
投稿日時 - 2001-06-26 12:49:37
このQ&Aのテーマ
このQ&Aで解決しましたか?
関連するQ&A
-PR-
-PR-
こんな書き方もあるよ!この情報は知ってる?あなたの知識を教えて!
このQ&Aにはまだコメントがありません。
あなたの思ったこと、知っていることをここにコメントしてみましょう。

その他の関連するQ&A、テーマをキーワードで探す

キーワードでQ&A、テーマを検索する
-PR-
-PR-
-PR-

特集


-PR-

ピックアップ

-PR-
ページ先頭へ