c - rand() を使用して、(a,b)、[a,b)、(a,b]、および [a,b] で均一に分散された浮動小数点数を生成します。

Question

4 種類の間隔すべてで乱数を生成する「最良の」方法を 1 か所に集めたいと考えています。私はこれをグーグルするのにうんざりしています。検索結果は多くのがらくたを出します。関連する結果でさえ、完全に間違っていることが多いページやブログであるか、自任の専門家が技術的な問題について互いに意見を異にする議論を行っているページやブログであり、多くの場合、彼らの「回答」は、さまざまなタイプについて知らないことを明らかにしているように見えます (非公開) 、オープン、セミオープン）間隔の。このような「単純な」質問に対して、C で乱数を生成することに関する悪い情報を読むのにうんざりしています。

一様に分散された浮動小数点数を生成する方法を教えてください。（a、b）、[a、b）、（a、b]、および[a、b]に対する私の典型的な方法（例として「ロングダブル」を使用）は次のとおりです。

long double a=VALUE1,b=VALUE2;
long double x1,x2,x3,x4;

srand((unsigned)time(NULL));

/* x1 will be an element of [a,b] */
x1=((long double)rand()/RAND_MAX)*(b-a) + a;

/* x2 will be an element of [a,b) */
x2=((long double)rand()/((long double)RAND_MAX+1))*(b-a) + a;

/* x3 will be an element of (a,b] */
x3=(((long double)rand()+1)/((long double)RAND_MAX+1))*(b-a) + a;

/* x4 will be an element of (a,b) */    
x4=(((long double)rand()+1)/((long double)RAND_MAX+2))*(b-a) + a;

単位間隔 (0,1)、[0,1)、(0,1]、および [0,1] の特殊なケースの場合:

long double x1,x2,x3,x4;

srand((unsigned)time(NULL));

/* x1 will be an element of [0,1] */
x1=((long double)rand()/RAND_MAX);

/* x2 will be an element of [0,1) */
x2=((long double)rand()/((long double)RAND_MAX+1));

/* x3 will be an element of (0,1] */
x3=(((long double)rand()+1)/((long double)RAND_MAX+1));

/* x4 will be an element of (0,1) */    
x4=(((long double)rand()+1)/((long double)RAND_MAX+2));

RAND_MAX と rand() の戻り値の両方のキャストが必要だと思います。これは、整数の除算を避けたいという理由だけでなく、それらが int であり、1 つ (または 2 つ) を追加するとオーバーフローする可能性があるためです。

「double」と「float」のバージョンはまったく同じですが、タイプを置き換えているだけだと思います。さまざまな浮動小数点型で発生する微妙な点はありますか?

上記の実装に問題はありますか? もしそうなら、何をどのように修正しますか？

編集：上記の実装は、それらが正しいために必要なテストに合格します（少なくとも64ビットLinuxを実行している64ビットIntel Core 2 Duoマシンでは）：x1は0と1の両方を生成でき、x2は0を生成できますが、されていませんx3 は 1 を生成できますが、0 を生成することは確認されておらず、x4 は 0 または 1 を生成することが確認されていません。

Answer 1

範囲内のすべての double を可能にし、その確率が隣接する double 値との差に比例するようにしたい場合、実際には非常に困難です。

範囲を検討してください[0, 1000]。範囲の非常に小さな最初の部分には、値の絶対的なバケットロードがあります: と の間に 100 万個の値があり0、1000000*DBL_MIN約DBL_MIN2 * 10 ^-308です。範囲内には複数の値があるため、すべてを生成するには2^32を 1 回呼び出すだけでは明らかに不十分です。rand()あなたがする必要があるのは、二重の仮数を一様に生成し、指数分布で指数を選択し、結果が範囲内にあることを確認するために少しファッジすることです。

範囲内のすべての double が可能である必要がない場合、「真の」連続一様ランダム分布では、正確な値が発生する確率はとにかく 0であるため、開いた範囲と閉じた範囲の違いはほとんど関係ありません。したがって、オープンレンジで数値を生成することもできます。

そうは言っても、あなたが提案した実装は、あなたが言う範囲内の値を生成し、閉じた範囲と半閉じた範囲では、エンドポイントを確率で生成し1/(RAND_MAX+1)ます。多くの、またはほとんどの実用的な目的には、これで十分です。

正確に表現できるRAND_MAX+2範囲内であれば、+1 と +2 をいじっても問題ありません。doubleこれは IEEE 倍精度と 32 ビットintに当てはまりますが、実際には C 標準では保証されていません。

(あなたの の使用はlong double少し混乱させるので無視しています。少なくとも と同じ大きさであることが保証されていますがdouble、 とまったく同じである一般的な実装があるdoubleため、longは不確実性以外には何も追加しません)。

Answer 2

問題が完全に特定されていないため、この質問には答える準備ができていません。特に、生成できる値のセットをどの程度細かく分散するかについての仕様は述べられていません。説明のために、[0, 1] の値を生成することを検討し、表現可能な値を持つ浮動小数点形式を検討してください。

0、1/16、2/16、3/16、4/16、6/16、8/16、12/16、1.

これらの値に対するいくつかの分布は、「一様」と見なされる場合があります。

• それぞれを同じ確率で選択します。これは離散値では均一ですが、値間の実際の距離では密度が均一ではありません。
• その付近の表現可能な値の密度に比例する確率でそれぞれを選択します。
• 間隔全体で同じ粒度を維持するには、0、4/16、8/16、12/16、および 1 を同じ確率で選択します。

これらの最初のものは意図されたものではないので、却下します。2 番目は Steve Jessop の提案に似ていますが、まだ完全には指定されていません。その中点から次の点までの間隔に比例した確率で 0 を選択する必要がありますか? (これにより、1/32 の確率が得られます。) または、-1/32 から 1/32 までの間隔を中心に関連付ける必要がありますか? (これは、1/17 の確率を与え、1 がそれ​​自体より 1/32 延長された間隔も割り当てられたと仮定します。)

これは閉区間なので、0 と 1 で停止する必要があると考えるかもしれません。しかし、あるアプリケーションで [0, 2] の分布を区間 [0, 1] と (1, 2]. 後者の 2 つの区間の分布の結合が、前者の区間の分布と等しくなるようにしたいので、分布はうまくかみ合う必要があります。

3 番目のケースにも同様の問題があります。おそらく、このような粒度を維持したい場合は、0 を 1/8 の確率で選択し、1/4、1/2、3/4 の 3 つの点をそれぞれ 1/4 の確率で選択し、1 を 1/8 の確率で選択する必要があります。 .

ジェネレーターの望ましいプロパティを指定するこれらの問題とは別に、質問者によって提案されたコードにはいくつかの問題があります。

• RAND_MAX+1 が 2 の累乗であると仮定すると (したがって、2 進浮動小数点演算ではそれで割るのは「適切」です)、RAND_MAX または RAND_MAX+2 で割ると、生成される値に不規則性が生じる場合があります。それらには奇妙な量子化があるかもしれません。

• 1/(RAND_MAX+1) ≤ 1/4 ULP(1) の場合、間隔が [0, 1) であるため、RAND_MAX/(RAND_MAX+1) は切り上げられて 1 を返します。(「ULP(1)」は、使用されている浮動小数点形式の値 1 の最小精度の単位を意味します。) (これは、RAND_MAX が仮数のビット内に収まる long double を使用したテストでは確認されていませんが、たとえば、RAND_MAX が 2147483647 で、浮動小数点型が float で仮数部が 24 ビットの場合に発生します。)

• 乗算(b-a)および加算によってa丸め誤差が発生し、その結果を評価する必要があります。b-aが小さい場合とaが大きい場合、aおよびがbゼロにまたがる場合 (したがって、より細かい結果が表現可能であっても、b 付近の粒度が失われる) など、さまざまなケースがあります。

• (0, 1) の結果の下限は、1/(RAND_MAX+2) に最も近い浮動小数点値です。この境界は、浮動小数点値の精度や目的の分布とは関係ありません。これは単に rand の実装の成果物です。(0, 1/(RAND_MAX+2)) 内の値は、問題の仕様に起因する理由なしに省略されています。同様のアーティファクトが上端に存在する可能性があります (特定の浮動小数点形式、rand の実装、および間隔のエンドポイント b によって異なります)。

質問者がこの「単純な」問題に対して不満足な回答に遭遇した理由は、それが単純な問題ではないということです。

Answer 3

まず、[a,b] で乱数を生成します。[a,b) で乱数を生成するには、[a,b] で乱数を生成し、それが b と等しいかどうかを確認し、そうであればもう一度試してください。他のすべての開区間バリアントについても同様です。

Answer 4

頭のてっぺんから、さまざまな浮動小数点型と整数型（テンプレート化されたC ++実装のボーナスポイント）のすべてのバリアントを提供し、rand()より良いものに置き換えます（drand48()頭に浮かぶ）

Answer 5

以下は、生成された数値の基本的なエラーを見つけるために使用している (非常に大雑把な) テストです。生成された数値が良好であることを示すことは意図されていませんが、悪くはないことを示しています。

#include<stdio.h>
#include<stdlib.h>
#include<time.h>

int main(int argc, char *argv[]) {

    long double x1,x2,x3,x4;
    if ( argc!=2 ) {
        printf("USAGE: %s [1,2,3,4]\n",argv[0]);
        exit(EXIT_SUCCESS);
    }

    srand((unsigned int)time(NULL));

    printf("This program simply generates random numbers in the chosen interval\n"
               "and looks for values on the boundary or outside it. When an\n"
               "allowable boundary is found, it reports it. Unexpected \"impossible\"\n"
               "values will be reported and the program will terminte. Under\n"
               "normal circumstances, the program should not terminate. Use ctrl-c.\n\n");

    switch ( atoi(argv[1]) ) {
        case 1:
            /* x1 will be an element of [0,1] */
            printf("NOTE: Testing [0,1].\n");
            while ( 1 ) {
                x1=((long double)rand()/RAND_MAX);
                if ( x1==0 ) {
                    printf("x1=0 ENCOUNTERED.\n");
                } else if ( x1==1 ) {
                    printf("x1=1 ENCOUNTERED.\n");
                } else if ( x1 < 0 ) {
                    printf("x1<0 ENCOUNTERED. Abnormal termination.\n");
                    exit(EXIT_FAILURE);
                } else if ( x1 > 1 ) {
                    printf("x1>0 ENCOUNTERED. Abnormal termination.\n");
                    exit(EXIT_FAILURE);
                }
            }
            break;
        case 2:
            /* x2 will be an element of [0,1) */
            printf("NOTE: Testing [0,1).\n");
            while ( 1 ) {
                x2=((long double)rand()/((long double)RAND_MAX+1));
                if ( x2==0 ) {
                    printf("x2=0 ENCOUNTERED.\n");
                } else if ( x2==1 ) {
                    printf("x2=1 ENCOUNTERED. Abnormal termination.\n");
                    exit(EXIT_FAILURE);
                } else if ( x2 < 0 ) {
                    printf("x2<0 ENCOUNTERED. Abnormal termination.\n");
                    exit(EXIT_FAILURE);
                } else if ( x2 > 1 ) {
                    printf("x2>0 ENCOUNTERED. Abnormal termination.\n");
                    exit(EXIT_FAILURE);
                }
            }
            break;
        case 3:
            /* x3 will be an element of (0,1] */
            printf("NOTE: Testing (0,1].\n");
            while ( 1 ) {
                x3=(((long double)rand()+1)/((long double)RAND_MAX+1));
                if ( x3==1 ) {
                    printf("x3=1 ENCOUNTERED.\n");
                } else if ( x3==0 ) {
                    printf("x3=0 ENCOUNTERED. Abnormal termination.\n");
                    exit(EXIT_FAILURE);
                } else if ( x3 < 0 ) {
                    printf("x3<0 ENCOUNTERED. Abnormal termination.\n");
                    exit(EXIT_FAILURE);
                } else if ( x3 > 1 ) {
                    printf("x3>0 ENCOUNTERED. Abnormal termination.\n");
                    exit(EXIT_FAILURE);
                }
            }
            break;
        case 4:
            /* x4 will be an element of (0,1) */
            printf("NOTE: Testing (0,1).\n");
            while ( 1 ) {
                x4=(((long double)rand()+1)/((long double)RAND_MAX+2));
                if ( x4==0 ) {
                    printf("x4=0 ENCOUNTERED. Abnormal termination.\n");
                    exit(EXIT_FAILURE);
                } else if ( x4==1 ) {
                    printf("x4=1 ENCOUNTERED. Abnormal termination.\n");
                    exit(EXIT_FAILURE);
                } else if ( x4 < 0 ) {
                    printf("x4<0 ENCOUNTERED. Abnormal termination.\n");
                    exit(EXIT_FAILURE);
                } else if ( x4 > 1 ) {
                    printf("x4>0 ENCOUNTERED. Abnormal termination.\n");
                    exit(EXIT_FAILURE);
                }
            }
            break;
        default:
            printf("ERROR: invalid argument. Enter 1, 2, 3, or 4 for [0,1], [0,1), (0,1], and (0,1), respectively.\n");
            exit(EXIT_FAILURE);
    }

    exit(EXIT_SUCCESS);
}