C中的快速双exp2函数

如何解决C中的快速双exp2函数

我需要以下双精度快速 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;
}

解决方法

我的假设是问题中的现有代码假设 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),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,sizeof res);
    return res;
}

最大相对误差和均方根误差分别略微降低至 4.93e-5 和 9.91e-6。这正如预期的那样，因为对于大致精确到 15 位的近似值，中间计算是使用 24 位精度还是 53 位精度执行无关紧要。计算使用除法，在我熟悉的所有平台上，double 的速度往往比 float 慢，因此双精度端口似乎没有提供任何显着优势，除了如果调用代码使用双精度计算，可能会消除转换开销。

C中的快速双exp2函数

如何解决C中的快速双exp2函数

解决方法

相关推荐