Random.nextGaussian() 的真正最大值(和最小值)是多少?随机实现答案仔细看看下一个高斯倒车下一个双倍将一切整合在一起

2022-09-03 07:04:04

从理论上讲,的边界意味着正无穷大和负无穷大。但是,由于 用于计算高斯随机数的 ()不会无限接近 0 和 1,因此 存在一个实际极限。而且也不是一个完全均匀的分布。nextGaussianRandom.nextDoublenextGaussianRandom.next

从理论上讲,最大值应约为2.2042 * 10 ^ 17,并且与(参考)的53位移位有关,但这可能只是一个上限。nextDouble

答案可能取决于 和 的分布和确切实现。我也找不到太多关于它的信息。Random.nextStrictMath.sqrtStrictMath.log

是的,我知道外部值极不可能,但它可能是相关的,例如在游戏中的RNG操作的背景下。


答案 1

随机实现

对于这个答案,您必须知道的最重要的事情是实现:Random.nextGaussian

synchronized public double nextGaussian() {
    // See Knuth, ACP, Section 3.4.1 Algorithm C.
    if (haveNextNextGaussian) {
        haveNextNextGaussian = false;
        return nextNextGaussian;
    } else {
        double v1, v2, s;
        do {
            v1 = 2 * nextDouble() - 1; // between -1 and 1
            v2 = 2 * nextDouble() - 1; // between -1 and 1
            s = v1 * v1 + v2 * v2;
        } while (s >= 1 || s == 0);
        double multiplier = StrictMath.sqrt(-2 * StrictMath.log(s)/s);
        nextNextGaussian = v2 * multiplier;
        haveNextNextGaussian = true;
        return v1 * multiplier;
    }
}

并实现:Random.nextDouble

public double nextDouble() {
    return (double) (((long)(next(26)) << 27) + next(27)) / (1L << 53);
}

首先,我想提请您注意一次生成 2 个值的事实,并且根据您是否知道自上次设置种子以来已通过的调用数,您可以对奇数与偶数调用使用略低的最大值。从现在开始,我将调用两个最大值v1_max和v2_max,指的是值是由 还是 生成的。nextGaussiannextGaussianv1 * multiplierv2 * multiplier

答案

有了这个,让我们直接切入追逐,稍后解释:

|      |Value             |Seed*          |
|------|------------------|---------------|
|v1_max|7.995084298635286 |97128757896197 |
|v2_max|7.973782613935931 |10818416657590 |
|v1_min|-7.799011049744149|119153396299238|
|v2_min|-7.844680087923773|10300138714312 |
* Seeds for v2 need to have nextGaussian called twice before you see the value listed.

仔细看看下一个高斯

@KaptainWutax和@Marco13的答案已经详细介绍了同样的事情,但我认为在图表上看到事情会让事情变得更清晰。让我们专注于v1_max,其他三个值具有非常相似的逻辑。我将在 x 轴、y 轴和 z 轴上绘制。v1v2v1 * multiplier

Graph

我们的眼睛立即跳到最大点 = 0, = 0, = 无穷大。但是,如果您在 do-while 循环中注意到,它明确不允许这种情况。因此,从图中可以清楚地看出,实际v1_max的值必须略高,但不会高得多。同样值得注意的是,对于任何> 0 的值,最大值为 = 0。v1v2v1 * multiplierv1v1v1 * multiplierv2

我们查找v1_max的方法是从零开始计数(或者更具体地说,从0.5开始计算生成它,按照 的实现,以2^-53的步长递增)。但是,只要知道,我们如何得到其他变量,以及v1nextDoublenextDoublev1v1 * multiplierv1

倒车下一个双倍

事实证明,知道调用的输出就足以确定当时生成它的对象的种子。直观地说,这是因为从实现来看,它“看起来”应该有2 ^ 54个可能的输出 - 但的种子只有48位。此外,有可能在比蛮力快得多的时间内恢复这个种子。nextDoubleRandomnextDoubleRandom

我最初尝试了一种天真的方法,方法是直接使用来获取种子的位,然后暴力破解剩余的21位,但这被证明太慢而没有用处。然后 SicksonFSJoe 给了我一种更快的方法,从单个调用中提取种子。请注意,要了解此方法的详细信息,您必须知道 的实现和一些模块化算术。next(27)nextDoubleRandom.next

private static long getSeed(double val) {
    long lval = (long) (val * (1L << 53));
    // let t = first seed (generating the high bits of this double)
    // let u = second seed (generating the low bits of this double)
    long a = lval >> 27; // a is the high 26 bits of t
    long b = lval & ((1 << 27) - 1); // b is the high 27 bits of u

    // ((a << 22) + c) * 0x5deece66d + 0xb = (b << 21) + d (mod 2**48)
    // after rearranging this gives
    // (b << 21) - 11 - (a << 22) * 0x5deece66d = c * 0x5deece66d - d (mod 2**48)
    // and because modular arithmetic
    // (b << 21) - 11 - (a << 22) * 0x5deece66d + (k << 48) = c * 0x5deece66d - d
    long lhs = ((b << 21) - 0xb - (a << 22) * 0x5deece66dL) & 0xffffffffffffL;

    // c * 0x5deece66d is 56 bits max, which gives a max k of 375
    // also check k = 65535 because the rhs can be negative
    for (long k = 65535; k != 376; k = k == 65535 ? 0 : k + 1) {
        // calculate the value of d
        long rem = (0x5deece66dL - (lhs + (k << 48))) % 0x5deece66dL;
        long d = (rem + 0x5deece66dL) % 0x5deece66dL; // force positive
        if (d < (1 << 21)) {
            // rearrange the formula to get c
            long c = lhs + d;
            c *= 0xdfe05bcb1365L; // = 0x5deece66d**-1 (mod 2**48)
            c &= 0xffffffffffffL;
            if (c < (1 << 22)) {
                long seed = (a << 22) + c;
                seed = ((seed - 0xb) * 0xdfe05bcb1365L) & 0xffffffffffffL; // run the LCG forwards one step
                return seed;
            }
        }
    }

    return Long.MAX_VALUE; // no seed
}

现在我们可以从 中获取种子,这样我们就可以迭代值而不是种子。nextDoublev1

将一切整合在一起

该算法的概述如下:

  1. 初始化(代表 1)为 0.5nd1nextDouble
  2. 当上限和我们当前v1_max尚未跨越时,请重复步骤 3-7
  3. 递增 2^-53nd1
  4. 计算自(如果存在),并生成 、 和seednd1nd2v1v2s
  5. 检查s
  6. 生成高斯,与v1_max进行比较
  7. 通过假设 = 0 设置新的上限v2

这是一个Java实现。如果需要,您可以自己验证我上面给出的值。

public static void main(String[] args) {
    double upperBound;
    double nd1 = 0.5, nd2;
    double maxGaussian = Double.MIN_VALUE;
    long maxSeed = 0;
    Random rand = new Random();
    long seed;
    int i = 0;
    do {
        nd1 += 0x1.0p-53;
        seed = getSeed(nd1);

        double v1, v2, s;
        v1 = 2 * nd1 - 1;

        if (seed != Long.MAX_VALUE) { // not no seed
            rand.setSeed(seed ^ 0x5deece66dL);
            rand.nextDouble(); // nd1
            nd2 = rand.nextDouble();

            v2 = 2 * nd2 - 1;
            s = v1 * v1 + v2 * v2;
            if (s < 1 && s != 0) { // if not, another seed will catch it
                double gaussian = v1 * StrictMath.sqrt(-2 * StrictMath.log(s) / s);
                if (gaussian > maxGaussian) {
                    maxGaussian = gaussian;
                    maxSeed = seed;
                }
            }
        }

        upperBound = v1 * StrictMath.sqrt(-2 * StrictMath.log(v1 * v1) / (v1 * v1));
        if (i++ % 100000 == 0)
            System.out.println(maxGaussian + " " + upperBound);
    } while (upperBound > maxGaussian);
    System.out.println(maxGaussian + " " + maxSeed);
}

最后要注意的一个问题是,此算法将为您提供 .要在 中使用它,您必须使用 的乘数(在上表中已经为您完成)RandomsetSeedRandom0x5deece66dL


答案 2

因此,我在这里要说的一切都纯粹是理论上的,我仍在研究一个GPU程序来扫描整个种子库。

nextGaussian() 方法就是这样实现的。

private double nextNextGaussian;
private boolean haveNextNextGaussian = false;

 public double nextGaussian() {

   if (haveNextNextGaussian) {

     haveNextNextGaussian = false;
     return nextNextGaussian;

   } else {

     double v1, v2, s;

     do {
       v1 = 2 * nextDouble() - 1;   // between -1.0 and 1.0
       v2 = 2 * nextDouble() - 1;   // between -1.0 and 1.0
       s = v1 * v1 + v2 * v2;
     } while (s >= 1 || s == 0);

     double multiplier = StrictMath.sqrt(-2 * StrictMath.log(s)/s);
     nextNextGaussian = v2 * multiplier;
     haveNextNextGaussian = true;
     return v1 * multiplier;

   }

 }

最有趣的部分必须在最后,[返回v1 *乘数]。由于 v1 不能大于 1.0D,因此我们需要找到一种方法来增加乘数的大小,如下所示。

double multiplier = StrictMath.sqrt(-2 * StrictMath.log(s)/s);

唯一的变量是“s”,可以安全地确定较低的“s”是乘数越大。都很好吗?让我们继续前进。

 do {
   v1 = 2 * nextDouble() - 1;   // between -1.0 and 1.0
   v2 = 2 * nextDouble() - 1;   // between -1.0 and 1.0
   s = v1 * v1 + v2 * v2;
 } while (s >= 1 || s == 0);

这告诉我们“s”必须属于]0,1[数字集,并且我们正在寻找的最低值仅略大于零。“S”用“v1”和“v2”的平方和来声明。要获得最小的理论值,v2 需要为零,v1 需要尽可能小。为什么是“理论”?因为它们是从 nextDouble() 调用生成的。不能保证种子基数包含这 2 个连续数字。

现在让我们玩得开心!

“v1”可以保持的最低值是双精度值的 epsilon,即 2^(-1022)。回到过去,要得到这样的数字,nextDouble需要生成(2^(-1022)+ 1)/ 2。

那是。。。非常非常令人不安。我不是专家,但我非常确定许多位会丢失,浮点错误是意料之中的。

下一个Double可能(绝对)不可能生成这样的值,但目标是找到一个尽可能接近该数字的值。

只是为了好玩,让我们做完整的数学运算来找到答案。StrictMath.log() 被实现为自然日志。我还没有研究它的精度,但让我们假设在这个层面上没有限制。最高的下一个高斯将计算为...

= (-2 * ln(v1 * v1) / (v1 * v1)) * v1 
= (-2 * ln(EPSILON^2) / (EPSILON^2)) * EPSILON

where EPSILON is equal to 2^(-1022).

信不信由你,我几乎找不到任何接受如此小数字的计算器,但我最终选择了这个高精度计算器

通过插入这个等式,

(-2 * ln((2^(-1022))^2) / ((2^(-1022))^2)) * (2^(-1022))

我得到了,

1.273479378356503041913108844696651886724617446559145569961266215283953862086306158E+311

相当大吧?井。。。它绝对不会那么大...但考虑到这一点很好。希望我的推理是有道理的,不要害羞地指出我犯的任何错误。

正如我在开始时所说,我正在开发一个程序,以暴力破解所有种子并找到实际的最低值。我会随时向你通报最新情况。

编辑:

很抱歉回复晚了。在大约10个小时内暴力破解了2 ^ 48个种子后,我发现了与地球计算机完全相同的答案。


推荐