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ベストアンサー
- 回答日時:
もうここまで来たら、「しつこいよう」
じゃなくて明らかにしつこいですが、
答えは変わらないにしても間違いがあったので…。
> ですが、今回は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回と
すべきである訳です。
だんだん難しい話になってきましたね。
コンピュータと言うよりは、数学の話みたいです。
これ以上はもうちょっと時間をかけて考えなければならなそうです。
自分でもうちょっと考えてみます。
ありがとうございました。
No.8
- 回答日時:
どうもしつこいようですが…。
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回で正解だと思います。
No.7
- 回答日時:
No.6
- 回答日時:
えー、通してみてみると、大体、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でもそうなりました。
たびたびありがとうございます。
一様性が失われる部分について理解しました。
これなら、long doubleの有効桁数内でやる限りは、大丈夫でしょう。
ところで、簡単に7次と言いましたが、後で考えてみると不安になりました。
long doubleが10進で35桁、32768が5桁だから、35/5=7と出してしまったのですが、内部では2進なので、実際には何次くらいまで計算するべきかを考えると、方法がわかりません。
浮動小数ではごみを加減しても、ほとんど何も起こりませんが、ごみになる前に止めるのがセオリーでしょう。
10進←→2進は色々と複雑なので、難しいですね。
No.5
- 回答日時:
えっと、試しにプログラム作りました。
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
- 回答日時:
> 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バイト
と仮定しています。
右シフトと言うよりは、指数部分の減算を期待しているのですが、これは間違いなのでしょうか。
確かに、ビット合成すれば正確な値は得られそうですが、それにはlong doubleの内部表現を完全に理解し、ビットフィールドを用いなければ、ならないような気がするのですが。
No.3
- 回答日時:
えっと、Fooky さんにちょっとした誤解があるような気がしますの
で、追加です。質問者の意図は、2進数で考えたとき、
0.<1回目のrand()のbitパターン><2回目のrand()のbitパターン>
としたいという意味です。(よね?)
10進数で考えれば、0x8000 とか 0x80000000 とかで割ると、非常
に桁数が増えて丸めが必要になるような気がしますが、内部表現は
2進数ですから、単に 15bit あるいは 31bit シフトされるだけで
す(というか、それを期待してのアルゴリズムです)。もちろん、
この内部表現を10進に変換して表示するときには丸めて表示されま
すが、計算途上では丸めは起きません。
あと、もし rand() が理想的な一様乱数であれば、それぞれ重なら
ないようにシフトしたものを加算しても、何も問題なく、結果は一
様です。サイコロの問題は同じ桁で加算しているから起きるのです。
No.2
- 回答日時:
分母が「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()では
改善されていると思います。
指摘ありがとうございます。
最大で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固有の関数なのでしょうか。
No.1
- 回答日時:
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ライブラリのサイトで見つけられませんでした。
今は手元にローカルライブラリがないので。
お探しのQ&Aが見つからない時は、教えて!gooで質問しましょう!
似たような質問が見つかりました
- 高校 レピュニット数の性質についてです。 レピュニット数とは、各桁が1のみの数で、以下1がk桁の数をRkと 3 2023/07/21 19:58
- Visual Basic(VBA) ExcelのVBAコードについて教えてください。 6 2022/06/08 12:55
- Visual Basic(VBA) ExcelのVBAコードについて教えてください。 1 2022/04/01 12:11
- Android(アンドロイド) 通電はするが画面がたちらないスマホのデータをどうやって再現するか? 4 2022/05/29 06:08
- Mac OS rmしたドットファイルのある場所 1 2023/07/18 22:45
- ブルーレイ・プレーヤー・レコーダー 学習リモコンSONY RM-PLZ430Dについて 4 2023/05/07 11:26
- UNIX・Linux Linuxについて質問です。 以下のhistoryの出力結果から、sedコマンドのファイル名tmp1 1 2023/02/03 20:11
- Visual Basic(VBA) 乱数の表示 4 2022/10/12 10:43
- K-POP これは何年前の写真なのか知りたいです! BTSに最近ハマって、色んな画像集めたりグッズを集めたりして 1 2023/03/14 07:31
- 統計学 Excelによるサンプルの拡大について 6 2023/08/22 16:03
関連するカテゴリからQ&Aを探す
おすすめ情報
デイリーランキングこのカテゴリの人気デイリーQ&Aランキング
-
16進数 加算 減算 C言語
-
C言語プログラミングにて、arct...
-
VB6.0での小数点の扱いについて
-
O(n log n)について2
-
VB6のFIX関数での誤差について
-
ExcelでPC(パソコン)によって...
-
C言語でセルオートマトンを作成...
-
時刻の比較
-
パソコンで階乗を計算
-
”/”を使わずに割り算したいんで...
-
floatの有効桁数
-
c languageで 簡単な質問があ...
-
2進数の0.2?
-
浮動小数点演算を固定小数点演...
-
VBAのINT関数について
-
ラズベリーパイ>MM-TXS03で温度...
-
EXCELの関数"STDEV(標準偏差)"...
-
Double型について
-
浮動小数演算は実行環境の変化...
-
VBAでの割り算の余りの求め方
マンスリーランキングこのカテゴリの人気マンスリーQ&Aランキング
-
O(n log n)について2
-
ExcelでPC(パソコン)によって...
-
ExcelのINT関数の計算結果がお...
-
16進数 加算 減算 C言語
-
VB.net Double と...
-
floatの有効桁数
-
三菱シーケンサ(Aシリーズ)で...
-
c languageで 簡単な質問があ...
-
除算を使わずに10で割りたい。
-
VBAでミリ秒まで出力する方法
-
VBAでの割り算の余りの求め方
-
VB6.0での小数点の扱いについて
-
VB6のFIX関数での誤差について
-
有効数字について 以前質問をし...
-
100桁の計算ができなくて困って...
-
浮動小数演算は実行環境の変化...
-
EXCELの関数"STDEV(標準偏差)"...
-
BCD・HEX・BINについて
-
コンピューターは指数関数をど...
-
乱数 なぜ剰余を使うのか
おすすめ情報