Поскольку двоичные числа с плавающей точкой IEEE-754 являются дробными, вычисление a b технически является алгебраической операцией. Однако общий подход к реализации powf(float a, float b) заключается в том, что e b * log a , то есть с использованием трансцендентных функций.
Однако есть несколько предостережений. log a не определено для a < 0, тогда как powf() допускает вычисления с некоторым отрицательным значением a. Экспонентация в форме expf() страдает от увеличения ошибки, как я объяснил в своем ответе на этот вопрос . Это требует, чтобы мы вычислили log a с точностью выше единичной для получения точного результата powf(). Существуют различные методы для достижения этой цели, простой способ состоит в том, чтобы использовать ограниченное количество двойных float вычислений, ссылки на которые я привел в своем ответе на этот вопрос . Суть double- float заключается в том, что каждый операнд с плавающей запятой представлен в виде пары значений float, называемых "head" и "tail", которые удовлетворяют соотношению | tail | ≤ ½ * ulp (| head |) при правильной нормализации.
Код c ниже показывает примерную реализацию этого подхода. Предполагается, что в IEEE 754-2008 доступна операция FMA (плавное умножение-сложение), которая представлена в C как стандартные математические функции fma(), fmaf(). Он не предусматривает обработку errno или исключений с плавающей запятой, но он обеспечивает правильную обработку всех 18 особых случаев, перечисленных в стандарте ISO C. Тесты были выполнены с включенной денормальной поддержкой; код может работать, а может и не работать должным образом в среде, не поддерживающей стандарт IEEE-754 (FTZ).
Часть возведения в степень использует простое сокращение аргумента, основанное непосредственно на полулогарифмическом кодировании чисел с плавающей точкой, затем применяет полиномиальное минимаксное приближение на интервале первичного приближения. Вычисление логарифма основано на log(x) = 2 atanh ((x-1) / (x+1)) в сочетании с выборочным использованием вычисления двойного float для достижения относительной точности 3,65e-9. Вычисление b * log a выполняется как операция с двойным float, точность окончательного возведения в степень повышается за счет линейной интерполяции, наблюдая, что e x является его собственной производной и, следовательно, e x + y ≅ e x + y ⋅ e x , когда | y | ≪ | x |.
Double- 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;
}