pow（浮动，浮动）算法

Question

我需要一个有效的算法在两个浮点数之间做数学::幂函数，你知道怎么做，（我需要算法不使用函数本身）

Answer 1

通用算法倾向于将浮动功率计算为整数幂与剩余根的组合。整数幂相当简单，可以使用Newton - Raphson method或Taylor series计算根。 IIRC numerical recipes in C就此有一些文字。还有其他（可能更好）的方法也可以这样做，但这将为实现一个令人惊讶的复杂问题提供一个合理的起点。另请注意，某些实现使用查找表和许多技巧来减少所需的计算。

Answer 2

由于IEEE-754二进制浮点数是分数，因此计算ab在技术上是一种代数运算。然而，实现powf(float a, float b)的常用方法是eb * log a，即使用超越函数。

但是，有一些警告。 log a未定义为a < 0，而powf()允许使用一些负a进行计算。 expf()形式的指数会出现误差放大，正如我在this question的回答中所解释的那样。这要求我们计算高于单精度的log a，以获得准确的powf()结果。有各种技术可以实现这一点，一种简单的方法是使用有限数量的双float计算，我在this question的答案中提供了参考。 double-float的本质是每个浮点操作数被表示为一对称为“head”和“tail”的float值，它们满足关系| tail |正确归一化时，≤½* ulp（| head |）。

下面的c代码显示了该方法的示例性实现。它假设IEEE 754-2008操作FMA（融合乘法 - 加法）可用，它在C中作为标准数学函数fma()，fmaf()公开。它不提供errno或浮点异常的处理，但它确实提供了对ISO C标准所列举的所有18种特殊情况的正确处理。已启用非正规支持执行测试;在非IEEE-754清零到零（FTZ）环境中，代码可能正常工作，也可能无法正常工作。

取幂部分采用直接基于浮点数的半对数编码的简单参数约简，然后在初级近似区间上应用多项式minimax approximation。对数计算基于log(x) = 2 atanh ((x-1) / (x+1))，结合选择性使用双float计算，实现3.65e-9的相对精度。 b * log a的计算是作为双float运算执行的，通过线性插值，通过观察ex是它自己的导数，因此ex +y≅ex+ y·ex，当| y时，最终取幂的精度得到改善。 | «| x |。

双float计算在溢出边界附近变得有点不对，在代码中有两个这样的实例。当exp的输入的头部刚刚导致结果溢出到无穷大时，尾部可能是负的，使得powf()的结果实际上是有限的。解决这个问题的一种方法是在这种情况下将“头部”的值减少一个ulp，另一种方法是通过在易于获得的圆形到零模式下添加来计算头部，因为这将确保类似的符号首尾。另一个警告是，如果取幂确实溢出，我们就不能插入结果，因为这样做会产生NaN。

应该注意的是，这里使用的对数计算的准确性不足以确保忠实地舍入powf()实现，但它提供了一个相当小的误差（我在广泛测试中发现的最大误差小于2 ulps）并且它为了证明相关的设计原则，允许代码保持相当简单。

#include <math.h> // for helper functions, e.g frexpf(), ldexpf(), isinff(), nextafterf()

/* Compute log(a) with extended precision, returned as a double-float value 
   loghi:loglo. Maximum relative error about 8.36545e-10.
*/
void my_logf_ext (float a, float *loghi, float *loglo)
{
    const float LOG2_HI =  0x1.62e430p-1f;  //  6.93147182e-1
    const float LOG2_LO = -0x1.05c610p-29f; // -1.90465421e-9

    float m, r, i, s, t, p, qhi, qlo;
    int e;

    /* Reduce argument to m in [181/256, 362/256] */
    m = frexpf (a, &e);
    if (m < 0.70703125f) { /* 181/256 */
        m = m + m;
        e = e - 1;
    }
    i = (float)e;

    /* Compute q = (m-1)/(m+1) as a double-float qhi:qlo */
    p = m + 1.0f;
    m = m - 1.0f;
    r = 1.0f / p;
    qhi = m * r;
    t = fmaf (qhi, -2.0f, m);
    s = fmaf (qhi, -m, t);
    qlo = r * s;

    /* Approximate atanh(q), q in [-75/437, 53/309] */ 
    s = qhi * qhi;
    r =             0x1.08c000p-3f;  // 0.1292724609
    r = fmaf (r, s, 0x1.22cde0p-3f); // 0.1419942379
    r = fmaf (r, s, 0x1.99a160p-3f); // 0.2000148296
    r = fmaf (r, s, 0x1.555550p-2f); // 0.3333332539
    t = fmaf (qhi, qlo + qlo, fmaf (qhi, qhi, -s)); // s:t = (qhi:qlo)**2
    p = s * qhi;
    t = fmaf (s, qlo, fmaf (t, qhi, fmaf (s, qhi, -p))); // p:t = (qhi:qlo)**3
    s = fmaf (r, p, fmaf (r, t, qlo));

    /* log(a) = 2 * atanh(q) + i * log(2) */
    t = fmaf ( LOG2_HI * 0.5f, i, qhi);
    p = fmaf (-LOG2_HI * 0.5f, i, t);
    s = (qhi - p) + s;                        
    s = fmaf ( LOG2_LO * 0.5f, i, s);
    r = t + t;
    *loghi = t = fmaf (2.0f, s, r);
    *loglo = fmaf (2.0f, s, r - t);
}

/* Compute exp(a). Maximum ulp error = 0.86565 */
float my_expf_unchecked (float a)
{
    float f, r;
    int i;

    // exp(a) = exp(i + f); i = rint (a / log(2))
    r = fmaf (0x1.715476p0f, a, 0x1.8p23f) - 0x1.8p23f; // 1.442695, 12582912.0
    f = fmaf (r, -0x1.62e400p-01f, a); // log_2_hi // -6.93145752e-1
    f = fmaf (r, -0x1.7f7d1cp-20f, f); // log_2_lo // -1.42860677e-6
    i = (int)r;
    // approximate r = exp(f) on interval [-log(2)/2,+log(2)/2]
    r =             0x1.694000p-10f; // 1.37805939e-3
    r = fmaf (r, f, 0x1.125edcp-7f); // 8.37312452e-3
    r = fmaf (r, f, 0x1.555b5ap-5f); // 4.16695364e-2
    r = fmaf (r, f, 0x1.555450p-3f); // 1.66664720e-1
    r = fmaf (r, f, 0x1.fffff6p-2f); // 4.99999851e-1
    r = fmaf (r, f, 0x1.000000p+0f); // 1.00000000e+0
    r = fmaf (r, f, 0x1.000000p+0f); // 1.00000000e+0
    // exp(a) = 2**i * exp(f);
    r = ldexpf (r, i);
    return r;
}

/* a**b = exp (b * log (a)), where a > 0, and log(a) is computed with extended 
   precision as a double-float.  The maximum ulp error found so far is 1.95083
   ulp for a = 0.698779583, b = 224.566528
*/
float my_powf_core (float a, float b)
{
    const float MAX_IEEE754_FLT = 0x1.fffffep127f;
    const float EXP_OVFL_BOUND = 0x1.62e430p+6f;
    const float EXP_OVFL_UNFL_F = 104.0f;
    const float MY_INF_F = 1.0f / 0.0f;
    float lhi, llo, thi, tlo, phi, plo, r;

    /* compute lhi:llo = log(a) */
    my_logf_ext (a, &lhi, &llo);
    /* compute phi:plo = b * log(a) */
    thi = lhi * b;
    if (fabsf (thi) > EXP_OVFL_UNFL_F) { // definitely overflow / underflow
        r = (thi < 0.0f) ? 0.0f : MY_INF_F;
    } else {
        tlo = fmaf (lhi, b, -thi);
        tlo = fmaf (llo, b, +tlo);
        /* normalize intermediate result thi:tlo, giving final result phi:plo */
#if FAST_FADD_RZ
        phi = __fadd_rz (thi, tlo);// avoid premature ovfl in exp() computation
#else // FAST_FADD_RZ
        phi = thi + tlo;
        if (phi == EXP_OVFL_BOUND){// avoid premature ovfl in exp() computation
            phi = nextafterf (phi, 0.0f);
        }
#endif // FAST_FADD_RZ
        plo = (thi - phi) + tlo;
        /* exp'(x) = exp(x); exp(x+y) = exp(x) + exp(x) * y, for |y| << |x| */
        r = my_expf_unchecked (phi);
        /* prevent generation of NaN during interpolation due to r = INF */
        if (fabsf (r) <= MAX_IEEE754_FLT) {
            r = fmaf (plo, r, r);
        }
    }
    return r;
}

float my_powf (float a, float b)
{
    const float MY_NAN_F = 0.0f / 0.0f;
    const float MY_INF_F = 1.0f / 0.0f;
    int expo_odd_int;
    float r;

    /* special case handling per ISO C specification */
    expo_odd_int = fmaf (-2.0f, floorf (0.5f * b), b) == 1.0f;
    if ((a == 1.0f) || (b == 0.0f)) {
        r = 1.0f;
    } else if (isnanf (a) || isnanf (b)) {
        r = a + b;  // convert SNaN to QNanN or trigger exception
    } else if (isinff (b)) {
        r = ((fabsf (a) < 1.0f) != (b < 0.0f)) ? 0.0f :  MY_INF_F;
        if (a == -1.0f) r = 1.0f;
    } else if (isinff (a)) {
        r = (b < 0.0f) ? 0.0f : MY_INF_F;
        if ((a < 0.0f) && expo_odd_int) r = -r;
    } else if (a == 0.0f) {
        r = (expo_odd_int) ? (a + a) : 0.0f;
        if (b < 0.0f) r = copysignf (MY_INF_F, r);
    } else if ((a < 0.0f) && (b != floorf (b))) {
        r = MY_NAN_F;
    } else {
        r = my_powf_core (fabsf (a), b);
        if ((a < 0.0f) && expo_odd_int) {
            r = -r;
        }
    }
    return r;
}