我的假设是问题中的现有代码假设 IEEE-754 二进制浮点计算，特别是将 C 的 float 类型映射到 IEEE-754 的 binary32 格式。
现有代码表明，只有在正常范围内的浮点结果才有意义：通过从下方钳位输入来避免次正常结果，并忽略溢出。所以对于float 计算，有效输入在区间 [-126, 128) 中。通过详尽的测试，我发现问题中函数的最大相对误差为7.16e-5，均方根（RMS）误差为1.36e-5。
我的假设是，对double 计算的期望更改应该将允许输入的范围扩大到 [-1022, 1024)，并且应该保持相同的相对精度。代码是以相当神秘的方式编写的。因此，作为第一步，我将其重新排列为更具可读性的版本。在第二步中，我调整了核心近似的系数，以最小化最大相对误差。这会产生以下 ISO-C99 代码：
/* compute 2**p, for p in [-126, 128). Maximum relative error: 5.04e-5; RMS error: 1.03e-5 */
float fastexp2 (float p)
{
const int FP32_MIN_EXPO = -126; // exponent of minimum binary32 normal
const int FP32_MANT_BITS = 23; // number of stored mantissa (significand) bits
const int FP32_EXPO_BIAS = 127; // binary32 exponent bias
float res;
p = (p < FP32_MIN_EXPO) ? FP32_MIN_EXPO : p; // clamp below
/* 2**p = 2**(w+z), with w an integer and z in [0, 1) */
float w = floorf (p); // integral part
float z = p - w; // fractional part
/* approximate 2**z-1 for z in [0, 1) */
float approx = -0x1.6e7592p+2f + 0x1.bba764p+4f / (0x1.35ed00p+2f - z) - 0x1.f5e546p-2f * z;
/* assemble the exponent and mantissa components into final result */
int32_t resi = ((1 << FP32_MANT_BITS) * (w + FP32_EXPO_BIAS + approx));
memcpy (&res, &resi, sizeof res);
return res;
}
系数的重构和重新调整导致精度略有提高，最大相对误差为 5.04e-5，RMS 误差为 1.03e-5。需要注意的是，浮点运算通常不具有关联性，因此任何浮点运算的重新关联，无论是通过手动代码更改还是编译器转换，都可能影响所述精度，需要仔细重新测试。
我预计不需要修改代码，因为它可以编译为通用架构的高效机器代码，这可以从 Compiler Explorer 的试用编译中看出，例如gcc with x86-64 或 gcc with ARM64。
此时，很明显需要更改哪些内容才能切换到double 计算。将float 的所有实例更改为double，将int32_t 的所有实例更改为int64_t，更改文字常量和数学函数的类型后缀，并将IEEE-754 的浮点格式特定参数binary32 更改为那些用于 IEEE-754 binary64。核心逼近需要重新调整，以便在核心逼近中充分利用双精度系数。
/* compute 2**p, for p in [-1022, 1024). Maximum relative error: 4.93e-5. RMS error: 9.91e-6 */
double fastexp2 (double p)
{
const int FP64_MIN_EXPO = -1022; // exponent of minimum binary64 normal
const int FP64_MANT_BITS = 52; // number of stored mantissa (significand) bits
const int FP64_EXPO_BIAS = 1023; // binary64 exponent bias
double res;
p = (p < FP64_MIN_EXPO) ? FP64_MIN_EXPO : p; // clamp below
/* 2**p = 2**(w+z), with w an integer and z in [0, 1) */
double w = floor (p); // integral part
double z = p - w; // fractional part
/* approximate 2**z-1 for z in [0, 1) */
double approx = -0x1.6e75d58p+2 + 0x1.bba7414p+4 / (0x1.35eccbap+2 - z) - 0x1.f5e53c2p-2 * z;
/* assemble the exponent and mantissa components into final result */
int64_t resi = ((1LL << FP64_MANT_BITS) * (w + FP64_EXPO_BIAS + approx));
memcpy (&res, &resi, sizeof res);
return res;
}
最大相对误差和均方根误差均略微下降至 4.93e-5 和 9.91e-6。这是预期的，因为对于大致精确到 15 位的近似值，中间计算是以 24 位精度还是 53 位精度执行的无关紧要。计算使用除法，在我熟悉的所有平台上，double 的速度往往比float 慢，因此双精度端口似乎没有提供任何显着优势，除了可能消除转换开销如果调用代码使用双精度计算。