使用 rand() 在 (a,b)、[a,b)、(a,b] 和 [a,b] 上生成均匀分布的浮点数答案

【问题标题】：Use rand() to generate uniformly distributed floating point numbers on (a,b), [a,b), (a,b], and [a,b]使用 rand() 在 (a,b)、[a,b)、(a,b] 和 [a,b] 上生成均匀分布的浮点数
【发布时间】：2012-09-07 18:03:33
【问题描述】：

我想在一个地方收集在所有四种类型的间隔上生成随机数的“最佳”方式。我厌倦了谷歌搜索。搜索结果出现了很多废话。甚至相关的结果也是页面或博客，这些页面或博客通常是完全错误的，或者在讨论中自封的专家在某些技术上存在分歧，通常他们的“答案”似乎暴露了他们不了解不同的类型（关闭、开、半开）的区间。对于这样一个“简单”的问题，我厌倦了阅读有关在 C 中生成随机数的不良信息。

请告诉我如何生成均匀分布的浮点数。这是我在 (a,b)、[a,b)、(a,b] 和 [a,b] 上的典型方式（以“long double”为例）：

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() 的返回值进行强制转换是必要的，不仅因为我们想避免整数除法，而且因为它们是整数，否则加一（或二）可能会使它们溢出。

我认为“double”和“float”的版本完全相同，只是替换了类型。不同的浮点类型有什么微妙之处吗？

您发现上述实现有什么问题吗？如果是这样，您将如何解决？

编辑：上述实现通过了必要的测试以使其正确（至少在运行 64 位 Linux 的 64 位 Intel Core 2 Duo 机器上）：x1 可以生成 0 和 1，x2 可以生成 0 但没有'未见生成 1，x3 可生成 1 但未见生成 0，x4 未见生成 0 或 1。

【问题讨论】：

您描述的测试不足以测试正确性。它们仅测试边界违规，并且仅针对浮点格式和一兰特实现间隔的特定情况。适当的测试还需要解决分布是否满足所需的均匀性标准以及其他格式是否可能发生绑定违规。您的问题没有指定您想要什么分布：浮点格式更精细，结果应该更精细，还是应该在整个间隔内保持粒度？
我使用了“必要”这个词的技术含义。我同意你的观点，测试是不够的。
我建议您将以下内容添加到问题规范中：在为区间 [a 生成值时，结果应该恰好是 a 的概率是多少, b) 和 [a, b]？ b 也是如此？（我假设 a 和 b 始终可以以目标格式表示。如果不是，则必须解决其他问题。）在 ( a, b)，尤其是在浮点格式的精细度发生变化的子区间中。一旦问题得到充分说明，就可以设计解决方案。

标签： c random

【解决方案1】：

如果您希望该范围内的每个双精度值都是可能的，并且概率与它与其相邻双精度值之间的差异成正比，那么实际上真的很难。

考虑范围[0, 1000]。在该范围的非常小的第一部分中存在绝对的值桶负载：0 和 1000000*DBL_MIN 之间的一百万个值，DBL_MIN 大约是 2 * 10^-308。该范围内总共有多个 2^32 值，因此很明显，一次调用 rand() 不足以生成所有值。你需要做的是均匀地生成你的双精度数的尾数，并选择一个指数分布的指数，然后稍微捏造一些东西以确保结果在范围内。

如果您不要求范围内的每个双精度数都是可能的，那么开放范围和封闭范围之间的差异是相当不相关的，因为在“真正的”连续均匀随机分布中，概率any 的确切值无论如何都是 0。所以你不妨只在开放范围内生成一个数字。

所有这一切：是的，您提出的实现生成的值在您所说的范围内，并且对于封闭和半封闭范围，它们以概率1/(RAND_MAX+1) 左右生成端点。这对于许多或大多数实际用途来说已经足够了。

只要RAND_MAX+2 在double 可以准确表示的范围内，您就可以使用+1 和+2。这对于 IEEE 双精度和 32 位 int 是正确的，但 C 标准实际上并不能保证这一点。

（我忽略了你对long double 的使用，因为它有点混淆了。它保证至少和double 一样大，但是在一些常见的实现中它与double 完全相同，所以long 除了不确定性没有添加任何东西）。

【讨论】：

关于生成随机尾数和随机指数的优点。我不同意完全跳过封闭间隔。我认为有时问题的逻辑表明应该使用封闭或半封闭间隔（例如，如果您正在模拟来自可以在边界上返回值的测量设备的数据）。您关于添加 1 和 2 的最后评论是我自己的理解变得模糊的地方。我不确定这样做在系统之间是否符合犹太教规。
在实践中很好，在任何“正常”系统上，double 的精度位比 int 的位多。非正常系统可能包括不寻常的架构，如 ILP64 大型机，其中int 是 64 位。因此，如果您想非常谨慎，那么您将记录和/或断言对实施细节的确切要求。如果您正在模拟来自测量设备的数据并且您非常关心，那么我认为您希望它根据设备的粒度分布，而不是均匀分布，也不是根据(long double)(b - a)/ RAND_MAX 的粒度分布。跨度>
另一个很好的观点。我的问题的解决方案似乎比我所说的黑白方式更复杂。你的答案可能是最好的。如果在一两周内没有人回答得更好，我会很满意地接受你的回复。谢谢你的时间。
断言要求很容易，顺便说一句。你只需要((double)RAND_MAX < (double)RAND_MAX + 1) && ((double)RAND_MAX + 1 < (double)RAND_MAX + 2) 来确保调整确实排除了端点，而不是没有任何区别。然后，万一有人试图为不满足您要求的平台编译您的代码，他们可以自行修复或向您提交错误报告，具体取决于您对代码的处理方式。
在 0 和 1000000*DBL_MIN 之间有超过一百万个双精度值。大约有 4.6 个 quintillion 值。对于浮点数，在 0 到 1000000*FLT_MIN 之间大约有 1.75 亿。

【解决方案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，以在区间内保持相同的粒度。

我怀疑第一个是故意的，我会忽略它。第二个类似于 Steve Jessop 的建议，但仍未完全指定。是否应该以与从它到中点到下一个点的间隔成比例的概率选择 0？（这将给出 1/32 的概率。）或者它是否应该与以它为中心的区间相关联，从 -1/32 到 1/32？（这将给它 1/17 的概率，假设 1 也被分配了一个超出其自身 1/32 的间隔。）

您可能会认为这是一个闭区间，因此它应该在 0 和 1 处停止。但是假设我们在某些应用程序中将 [0, 2] 上的分布切分为区间 [0, 1] 和(1, 2]。我们希望后两个区间上的分布联合等于前一个区间上的分布。所以我们的分布应该很好地啮合。

第三种情况也有类似的问题。也许，如果我们希望保持这样的粒度，应该以 1/8 的概率选择 0，以 1/4 的概率选择 1/4、1/2 和 3/4 这三个点，以 1/8 的概率选择 1 .

除了指定生成器所需属性的这些问题之外，提问者提出的代码还有一些问题：

假设 RAND_MAX+1 是 2 的幂（因此除以它在二进制浮点算术中“很好”），除以 RAND_MAX 或 RAND_MAX+2 可能会导致生成的值出现一些不规则性。其中可能有奇怪的量化。
当 1/(RAND_MAX+1) ≤ 1/4 ULP(1) 时，RAND_MAX/(RAND_MAX+1) 将四舍五入并在不应该返回 1 时返回 1，因为区间为 [0, 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）。

我提交提问者对这个“简单”问题的回答不满意的原因是它不是一个简单的问题。

【讨论】：

您提出了一些新的有效观点，但我很困惑为什么您写道我“遇到了不令人满意的答案”，尤其是当我写道我对 Jessop 的回答“满意”时？我还注意到这个问题比最初提出的要复杂。
@Dr.PersonPersonII：我想，Eric 的意思是“我厌倦了谷歌搜索”。
@SteveJessop：我不明白。我是不是做错了什么？
@Dr.PersonPersonII：嗯，我不确定在 SO 上是否鼓励发誓，但这不是我的意思。我的意思是，考虑到您对体验的评价，当您在 Google 上搜索时，您“遇到不满意的答案”是合理的 :-)
@Dr.PersonPersonII：正如史蒂夫·杰索普所说，你的问题陈述表明你厌倦了谷歌搜索，并且搜索结果包含很多无价值的材料。我的评论表明，这些观察结果是由于问题并不简单。

【解决方案3】：

首先，在 [a,b] 上生成随机数。要在 [a,b) 上生成随机数，只需在 [a,b] 上生成一个随机数，检查它是否等于 b，如果是，请重试。对于所有其他开区间变体也是如此。

【讨论】：

这是一个很好的观点，我过去曾使用过它。但是上述方法行不通吗？它们是单行语句，因此比使用边界测试要优雅得多。

【解决方案4】：

我只想提供不同浮点和整数类型的所有变体（模板化 C++ 实现的奖励点），我会用更好的东西替换 rand() (drand48()想到了）

【讨论】：

【解决方案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);
}

【讨论】：