-1

我需要以下双精度快速 exp2 函数的版本，你能帮帮我吗？请不要回答说这是一个近似值，所以双重版本毫无意义，将结果转换为 (double) 就足够了..谢谢。我在某处找到并且对我有用的函数如下，它比 exp2f() 快得多，但不幸的是我找不到任何双精度版本：

inline float fastexp2(float p)
{
 if(p<-126.f) p= -126.f;
 int w=(int)p;
 float z=p-(float)w;
 if(p<0.f) z+= 1.f;
 union {uint32_t i; float f;} v={(uint32_t)((1<<23)*(p+121.2740575f+27.7280233f/(4.84252568f -z)-1.49012907f * z)) };
 return v.f;
}
4

1 回答 1

4

我的假设是问题中的现有代码假设 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 与 x86-64或gcc 与 ARM64。

在这一点上，很明显需要改变什么才能切换到double计算。更改 to 的所有实例，floattodouble的所有实例int32_t，int64_t更改文字常量和数学函数的类型后缀，并将 IEEE-754 的浮点格式特定参数更改为 IEEE-754binary32的参数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我熟悉的所有平台都慢，所以双精度端口似乎没有提供任何显着优势，除了如果调用代码使用双精度可能消除转换开销-精密计算。

于 2021-01-04T11:53:22.870 回答