Мне нужен эффективный алгоритм для выполнения функции math:: power между двумя поплавками,
вы знаете, как это сделать (мне нужно, чтобы алгоритм не использовал эту функцию)
Ответ 2
Поскольку двоичные числа с плавающей точкой IEEE-754 являются дробными, вычисление a b является технически алгебраической операцией. Однако общий подход к реализации powf(float a, float b) такой же, как e b * log a то есть с использованием трансцендентных функций.
Однако есть несколько предостережений. log a не определено для a < 0, тогда как powf() позволяет вычислять с некоторым отрицательным значением a. Экспонентация в форме expf() страдает от увеличения ошибки, как я объяснил в своем ответе на этот вопрос. Это требует от нас вычисления log a с точностью, превышающей единичную, для получения точного результата powf(). Существуют различные методы для достижения этой цели, простой способ заключается в использовании ограниченного количества вычислений double- с 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)) сочетании с выборочным использованием вычисления double- с float для достижения относительной точности 3,65e-9. Вычисление b * log a выполняется как операция с float double-, точность окончательного возведения в степень повышается за счет линейной интерполяции, наблюдая, что e x является собственной производной и, следовательно, e x + y ≅ e x + y ⋅ e x когда | y | X | x |.
double- вычисление с float становится немного сомнительным вблизи границ переполнения, в коде есть два примера этого. Когда головная часть ввода в exp просто вызывает переполнение результата до бесконечности, хвостовая часть может быть отрицательной, так что результат powf() на самом деле конечен. Одним из способов решения этой проблемы является уменьшение значения "головы" на один ulp в таком случае, и альтернативой является вычисление головы путем сложения в режиме округления до нуля, где это легко доступно, поскольку это обеспечит аналогичные знаки для голова и хвост. Другое предостережение состоит в том, что если возведение в степень действительно переполняется, мы не можем интерполировать результат, так как это создаст NaN.
Следует отметить, что точность вычислений логарифма, использованная здесь, недостаточна для обеспечения точно округленной реализации powf(), но она обеспечивает достаточно небольшую ошибку (максимальная ошибка, которую я обнаружил при экстенсивном тестировании, составляет менее 2 ульт). и это позволяет сохранять код достаточно простым для демонстрации соответствующих принципов проектирования.
#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;
}
