查找大于双精度值的最小浮点数

Question

我在 cudnnBatchNormalizationForwardTraining 函数中使用 CUDNN_BN_MIN_EPSILON 值时遇到了问题（请参阅 the docs here），事实证明这是因为我传递了 float 值 1e-5f 而不是 double（我'正在使用 float 值以节省内存并加快计算速度），并且该值一旦转换为浮点数就略小于 1e-5，这是该常量的实际值。

经过反复试验，我找到了一个不错的近似值：

const float CUDNN_BN_MIN_EPSILON = 1e-5f + 5e-13f;

我确信有更好的方法来解决这样的问题，所以问题是：

给定一个正的double 值，找到最小可能的float 值的最佳（如“可靠”）方法是什么（单独和如果/当转换为double 时）严格更大比最初的double 值？

另一种表述这个问题的方法是，给定一个double 值d1 和一个float 值f1，d1 - (float)f1 应该是最小可能的 negative 值（如否则，这意味着f1 小于 d1，这不是我们想要的）。

我做了一些基本的试验和错误（使用1e-5 作为我的目标值）：

// Check the initial difference
> 1e-5 - 1e-5f
2,5262124918247909E-13 // We'd like a small negative value here

// Try to add the difference to the float value
> 1e-5 - (1e-5f + (float)(1e-5 - 1e-5f))
2,5262124918247909E-13 // Same, probably due to approximation

// Double the difference (as a test)
> 1e-5 - (1e-5f + (float)((1e-5 - 1e-5f) * 2))
-6,5687345259044915E-13 // OK

使用这个近似值，最终的float 值为1,00000007E-05，看起来不错。

但是，* 2 乘法对我来说完全是任意的，我不确定它是否可靠或在那里做的最佳可能的事情。

有没有更好的方法来实现这一点？

谢谢！

---

编辑：这是我现在使用的（坏）解决方案，很乐意用更好的解决方案替换它！

/// <summary>
/// Returns the minimum possible upper <see cref="float"/> approximation of the given <see cref="double"/> value
/// </summary>
/// <param name="value">The value to approximate</param>
public static float ToApproximatedFloat(this double value)
    => (float)value + (float)((value - (float)value) * 2);

---

解决方案：这是最终的正确实现（感谢 John Bollinger）：

public static unsafe float ToApproximatedFloat(this double value)
{
    // Obtain the bit representation of the double value
    ulong bits = *((ulong*)&value);

    // Extract and re-bias the exponent field
    ulong exponent = ((bits >> 52) & 0x7FF) - 1023 + 127;

    // Extract the significand bits and truncate the excess
    ulong significand = (bits >> 29) & 0x7FFFFF;

    // Assemble the result in 32-bit unsigned integer format, then add 1
    ulong converted = (((bits >> 32) & 0x80000000u)
                        | (exponent << 23)
                        | significand) + 1;

    // Reinterpret the bit pattern as a float
    return *((float*)&converted);
}

Answer 1

由于您似乎对表示级别的细节感兴趣，您将依赖于类型 float 和 double 的表示。然而，在实践中，这很可能归结为IEEE-754 的基本“binary32”和“binary64”格式。它们具有一个符号位、几位有偏指数和一堆有效位的一般形式，对于标准化值，包括一个隐式位有效位。

简单案例

给定一个IEEE-754 binary64格式的double，其值不小于+2^-126，你要做的是

• 以可直接检查和操作的形式获取原始double 值的位模式。例如，作为无符号 64 位整数。

  double d = 1e-5;
  uint64_t bits;
  memcpy(&bits, &d, 8);
  

• 提取并重新偏置指数场

  uint64_t exponent = ((bits >> 52) & 0x7FF) - 1023 + 127;
  

• 提取有效位并截断多余的位

  uint64_t significand = (bits >> 29) & 0x7fffff;
  

• 将结果组装成 32 位无符号整数格式

  uint32_t float_bits = ((bits >> 32) & 0x80000000u)
          | (exponent << 23)
          | significand;
  

• 添加一个。由于您想要一个严格大于原始double 的结果，因此无论所有截断的有效位是否为0，这都是正确的。如果加法溢出有效位，它将正确增加指数字段。但是，它可能会产生无穷大的位模式。

  float_bits += 1;
  

• 将位模式存储/复制/重新解释为float的位模式

  float f;
  
  memcpy(&f, &float_bits, 4);
  

负数

给定一个二进制64格式的负double，其幅度不小于2^-126，按照上述过程，除了从float_bits中减去1而不是加1。请注意，对于 -2^-126，这会产生一个次正规的 binary32（见下文），这是正确的结果。

零和非常小的数字，包括次正规数

IEEE 754 提供了非常小的非零数的降低精度表示。这样的表示被称为subnormal。在某些情况下，超过给定输入 binary64 的最小 binary32 是次正规的，包括对于一些不是 binary64 次正规的输入。

此外，IEEE 754 提供有符号零，-0 是一种特殊情况：严格大于 -0（任一格式）的最小 binary32 是最小的正次正规数。注意：不是 +0，因为根据 IEEE 754，+0 和 -0 通过普通比较运算符比较相等。最小正、非零、次正规 binary32 值的位模式为 0x00000001。

受这些考虑的 binary64 值具有偏差 binary64 指数字段，其值小于或等于 binary64 指数偏差和 binary32 指数偏差 (896) 之间的差值。这包括那些偏置指数正好为 0 的那些，它们表征了 binary64 零和次正规。检查简单案例程序中的重新偏置步骤应该会导致您正确地得出结论，该程序将针对此类输入产生错误的结果。

这些案例的代码留作练习。

无穷大和 NaN

具有偏置 binary64 指数字段集的所有位的输入表示正或负无穷大（当 binary64 有效位没有设置位时）或非数字 (NaN) 值。 Binary64 NaN 和正无穷应该转换为它们的 binary32 等价物。负无穷大可能应该转换为最大幅度的负 binary32 值。这些需要作为特殊情况处理。

这些案例的代码留作练习。

Answer 2

在 C 中：

#include <math.h>

float NextFloatGreaterThan(double x)
{
    float y = x;
    if (y <= x) y = nexttowardf(y, INFINITY);
    return y;
}

如果您不想使用库例程，则将上面的nexttowardf(y, INFINITY) 替换为-NextBefore(-y)，其中NextBefore 取自this answer 并进行了修改：

• 将double 更改为float 并将DBL_ 更改为FLT_。
• 将.625 更改为.625f。
• 将fmax(SmallestPositive, fabs(q)*Scale) 替换为SmallestPositive < fabs(q)*Scale ? fabs(q)*Scale : SmallestPositive。
• 将fabs(q) 替换为(q < 0 ? -q : q)。

（显然，例程可以从-NextBefore(-y) 转换为NextAfter(y)。留给读者作为练习。）