完美随机浮点数

3174 字

16 分钟

完美随机浮点数

2025-05-12

未分类

无标签

翻译自 https://specbranch.com/posts/fp-rand/

完美随机浮点数#

许多编程语言和库中常用的一种流程并不是真正的浮点算法：

根据格式的精度生成随机整数。
转换为浮点数。
除以某个值以生成一个介于 0 和 1 之间的输出。

在代码中，这看起来像：

func (r *Rand) Float64() float64 {
    randInt := r.Int63n(1 << 53)
    return float64(randInt) / (1 << 53)
}

此函数应该生成在区间 $[0, 1)$ 上均匀分布的浮点数。0 是可能的输出，但 1 不是，并且分布是均匀的。上面算法中的 “53” 是根据浮点感知选择的：双精度浮点数具有 53 位精度，因此该算法仅创建与数字系统精度相等的位。看起来它符合要求。

然而，该算法存在一些基本缺陷。首先，它无法访问大多数介于 0 和 1 之间的浮点数。该范围内有 $2^{62}$ 个 float64 数字，而此算法仅能访问 $2^{53}$ 个：仅占空间的 $\frac{1}{512}$ 。即使我们使用更大的整数范围（如 $2^{64}$ ），许多随机整数也会映射到相同的浮点结果，并且一些整数在除法过程中会四舍五入到 1。

其次，来自此生成器的最低有效位是有偏的。这实际上是一个定点随机数生成器，随后进行浮点转换。

过去曾有一些尝试来修复这些缺陷，但在避免巨大性能损失的同时都未能成功。我最近完成了一篇关于生成完美随机浮点数的新方法的论文，并在 GitHub 上发布了预印本和 Go 语言参考实现。多亏了分支预测，该完美算法的速度并不比上述定点算法慢多少。

0 到 1 之间的浮点数#

浮点数是三个组件的元组，打包到一个 32 位或 64 位的字中：

字段	含义	位数 (`float32`)	位数 (`float64`)
符号 ( $s$ )	正号或负号	1	1
指数 ( $e$ )	最高有效位的量级	8	11
尾数 ( $m$ )	除最高有效位后的所有有效位	23	52

每个数字的存储都是高效的：最高有效位不存储，因为它可以根据指数推断为 0（当 $e = 0$ 时，表示 0）或 1（其他情况）。浮点数还有无穷大和非数字值（NaN 表示 “不是数字”），它们用全为 1 的指数编码。否则，指数是一个无符号数，通过减去偏置可以得到负指数。

一个正常的浮点数的值由以下公式给出：

F = (-1)^s \, 2^{e - eBias} \,(1.m)

通常， $eBias$ 是指数字段可表示范围的一半。32 位浮点数具有 8 位指数，使用 $eBias = 127$ ，因此可表示的指数范围为 -126 到 127（ $e = 0$ 和 $e = 255$ 有特殊含义）。

一半的浮点数位于 -1 和 1 之间：所有指数在 0 到 $eBias$ 之间的数字都小于 1。这似乎违反直觉，但这是浮点数机制的产物：它们在不同量级之间具有相同的精度。数字 $0.000001$ 和 $10^{20}$ 都具有 53 位精度（如果是 32 位浮点则为 24 位）。

这类似于“有效数字”的概念：每个浮点数无论大小都有 53 位有效数字。浮点数按幂律分布，极端量级的数字比看起来更多。

浮点舍入#

当操作产生的位数过多时，我们需要进行舍入以得到浮点结果。有四种浮点舍入模式，但其中三种相对专业：

最近舍入-偶数舍入（Round-to-nearest-ties-to-even）：舍入到最接近的浮点数，遇到平局时舍入到偶数。与学校学到的舍入类似，但学校通常将平局舍入到远离零的方向。偶数舍入对数值稳定性更好。
向零舍入（Round toward zero，截断）：始终舍入到更接近零的浮点数，等价于截断。
向下舍入（Round down）：始终舍入到较小的数字。
向上舍入（Round up）：始终舍入到较大的数字。

最近-偶数舍入是默认的算术舍入，因为它最准确。对于原生浮点随机数生成，我们可能更倾向于遵循所选的舍入模式，而不是严格地在 $[0,1)$ 上操作。下面的数轴示例展示了在假设浮点格式 ( $p=2$ ) 下的向下舍入和最近舍入模式。

浮点的整数算术#

浮点数被打包到机器字中，可以使用按位逻辑操作和整数算术进行操作，这使得我们可以做一些稍显危险的事情：

我们可以通过按位或（OR）将整数与浮点数的不同部分组合来构造浮点数（例如将尾数附加到指数）。
向浮点数添加小整数是按单位最后一位（ULP）进行微扰，当跨越指数边界时可能会出现奇异的行为。
添加或减去到达指数字段的大整数是快速乘除浮点数的一种方式，只要正确处理边界条件。

这些技巧允许你构建自己的快速操作。更疯狂的技巧，如 Quake III 的快速逆平方根算法，正是基于这些基本构建块。

浮点数的均匀随机性#

基于此，我们可以定义完美浮点随机数生成的概念：

对于无偏随机浮点数，按概率生成浮点数，其概率由从实数范围抽取一个实数然后舍入到浮点数决定。

整数和定点随机数生成都遵循此规则，但它们将空间分割为均匀的片段，因此每个数字的概率相等。浮点数没有此优势。相反，概率取决于浮点数字之间的距离。

通常，位于 $0.5$ 到 $1$ 之间的数字出现概率是位于 $0.25$ 到 $0.5$ 之间数字的两倍，后者又是位于 $0.125$ 到 $0.25$ 之间数字的两倍，以此类推。每个指数范围边界的概率由舍入模式决定：向下舍入（或向零舍入）时，范围底部与范围内其他数字概率相同；向上舍入时，范围顶部与范围内其他数字概率相同；最近舍入时，我们取两者的平均。

下面给出 $p$ 位精度（64 位浮点 $p=53$ ，32 位浮点 $p=24$ ）下前两个指数范围在 $(0,1)$ 上均匀生成的概率表：

输出范围	舍入模式	最低数字	中间数字	最高数字
$begin:math:display$ 0.5, 1 $end:math:display$	向下	$2^{-p}$	$2^{-p}$	0
$begin:math:display$ 0.5, 1 $end:math:display$	向上	$2^{-p-1}$	$2^{-p}$	$2^{-p}$
$begin:math:display$ 0.5, 1 $end:math:display$	最近	$2^{-p-1} + 2^{-p-2}$	$2^{-p}$	$2^{-p-1}$
$begin:math:display$ 0.25, 0.5 $end:math:display$	向下	$2^{-p-1}$	$2^{-p-1}$	$2^{-p}$
$begin:math:display$ 0.25, 0.5 $end:math:display$	向上	$2^{-p-2}$	$2^{-p-1}$	$2^{-p-1}$
$begin:math:display$ 0.25, 0.5 $end:math:display$	最近	$2^{-p-2} + 2^{-p-3}$	$2^{-p-1}$	$2^{-p-2} + 2^{-p-1}$

表格继续向指数为零方向展开。对于所有三种模式，范围内的中间数字概率相同：它们上下的实数区域相等，因此舍入不影响它们的概率。但范围的边界会变化。当向下舍入时，二的幂的概率与上方数字相同；向上舍入时，二的幂的概率与下方数字相同；最近舍入时，我们取上下概率的平均，使它们比下一级范围内的数字概率高 1.5 倍。

典型的随机浮点生成更贴合“向下舍入”模式，此模式不会生成 1。其他舍入模式会给 1 一定概率。另一个有趣的事实是，这些数字的最低有效位应当均匀分布，尽管最低有效位的幅度不同。

因为我们把它视为实数范围，无论区间开闭，四舍五入的概率质量相同。

浮点随机算法#

我们可以通过一个两步算法来匹配上述概率：

用 $p$ 位精度的定点随机数生成器生成一个固定点结果，如果生成 0，则特殊处理，用于确定实数范围。
通过回填额外的精度位来完成最终的随机数抽取。

定点阶段确定浮点数字范围，完成阶段从该范围中选取数字。在定点阶段，只有最低范围（从 0 到 $2^{-p}$ ）包含多个指数。要处理此范围，我们递归运行算法，但缩放因子为 $2^{-p}$ 。

当我们达到次正规数时，进入基底情况，此时定点阶段的所有可能值映射到同一个指数。这允许我们通过重复定点抽取来确定输出的指数。

定点结果具有已知数目的尾部零，基于结果的量级。完成阶段填充这些位。

在完成阶段，我们通过生成一个位数等于从 1 的指数到定点数指数差值的整数来回填精度位。这差值就是将定点数转换为浮点数时要填补的精度位，回填它们即可。此阶段根据舍入模式略有不同。

以下是 Go 语言中 float64 的完整算法实现。我们使用略有不同的定点随机数生成：不生成大整数然后除法，而是生成随机尾数，将其放入 [1,2) 范围，然后减 1。其输出分布与除法方法完全相同，但在某些架构上可能较慢。

// Generate float64 with a given rounding mode
func (r *FPRand) Float64Rounding(mode RoundingMode) float64 {
    const PRECISION = 53
    const FP64_MANTISSA_SIZE = PRECISION - 1
    const FP64_MANTISSA_MASK = (1 << FP64_MANTISSA_SIZE) - 1

    // PHASE 1: 缩放定点阶段

    // 用于缩小指数范围的变量
    exp_range := uint64(FP64_MANTISSA_SIZE) << FP64_MANTISSA_SIZE
    var one = math.Float64bits(1)

    // 跟踪下溢时的尾部位数
    var tail_bits uint64 = 0

    // 生成尾数，直到非零（预期只运行一次）
    var mantissa = r.src.Int63n(1 << FP64_MANTISSA_SIZE)
    for mantissa == 0 {
        // 缩小范围，每次循环都在上一次的 0 和 1 之间
        one -= exp_range

        if one < exp_range {
            // 低概率下溢基底情况，约 19 次循环后
            const UF_TAIL_SIZE = 35
            const UF_SHIFT = FP64_MANTISSA_SIZE - UF_TAIL_SIZE

            mantissa = r.src.Int63n(1 << UF_TAIL_SIZE) << UF_SHIFT
            tail_bits = UF_SHIFT
            break
        }

        mantissa = r.src.Int63n(1 << FP64_MANTISSA_SIZE)
    }

    // 生成一个位于 [1,2) 之间的浮点数并减 1，得到定点结果
    num := math.Float64frombits(one | mantissa) - math.Float64frombits(one)
    num_as_int := math.Float64bits(num)

    // PHASE 2: 完成阶段

    // 通过减去指数确定需要回填的尾部位数
    tail_bits += (one >> FP64_MANTISSA_SIZE) - (num_as_int >> FP64_MANTISSA_SIZE)

    // 处理指数为 0 的情况（极低概率）
    if tail_bits > FP64_MANTISSA_SIZE {
        tail_bits = FP64_MANTISSA_SIZE
    }

    // 根据舍入模式生成尾部
    var tail uint64 = 0
    if mode == RoundDown {
        tail = r.src.Int63n(1 << tail_bits)
    } else if mode == RoundUp {
        tail = r.src.Int63n(1 << tail_bits) + 1
    } else {
        tail = (r.src.Int63n(1 << tail_bits) + 1) >> 1
    }

    return math.Float64frombits(num_as_int + tail)
}

// 功能等价于 Float64 - 在 [0,1) 上生成数字
func (r *FPRand) Float64() float64 {
    return Float64Rounding(RoundDown)
}