самодельных пау() C++

Question

Я читал через How can I write a power function myself? и ответ дается dan04 привлекло мое внимание главным образом потому, что я не уверен, что ответ дается Fortran, но я взял, что и реализовано это:

#include <iostream> 
using namespace std; 
float pow(float base, float ex){ 
    // power of 0 
    if (ex == 0){ 
     return 1; 
    // negative exponenet 
    }else if(ex < 0){ 
     return 1/pow(base, -ex); 
    // even exponenet 
    }else if ((int)ex % 2 == 0){ 
     float half_pow = pow(base, ex/2); 
     return half_pow * half_pow; 
    //integer exponenet 
    }else{ 
     return base * pow(base, ex - 1); 
    } 
} 
int main(){ 
    for (int ii = 0; ii< 10; ii++){\ 
     cout << "pow(" << ii << ".5) = " << pow(ii, .5) << endl; 
     cout << "pow(" << ii << ",2) = " << pow(ii, 2) << endl; 
     cout << "pow(" << ii << ",3) = " << pow(ii, 3) << endl; 
    } 
} 

, хотя я не уверен, что я перевел это право, потому что все вызовы, дающие .5, как показатель экспоненты, возвращают 0. В ответе он утверждает, что ему может понадобиться log2 (x) на основе a^b = 2^(b * log2(a)), но я не уверен, потому что я не уверен, где это сказать, или если я даже думаю об этом праве.

ПРИМЕЧАНИЕ: Я знаю, что это может быть определено в математической библиотеке, но мне не нужны все дополнительные затраты на всю математическую библиотеку для нескольких функций.

EDIT: кто-нибудь знает реализацию с плавающей запятой для дробных показателей? (Я видел двойную реализацию, но это использовало трюк с регистрами, и мне нужна плавающая точка, и добавление библиотеки просто для того, чтобы сделать трюк, мне было бы лучше, включая математическую библиотеку)

Answer 1

Ниже приведены ссылки на реальные реализации powf.

Я ожидаю, что более простые решения не будут иметь точность в результате или не будут обрабатывать параметры InF и NaN.

http://opensource.apple.com/source/Libm/Libm-2026/Source/Intel/expf_logf_powf.c

http://opensource.apple.com/source/Libm/Libm-315/Source/ARM/powf.c

Answer 2

Я думаю алгоритм, который вы ищете, может быть «nth root». С начальным приближением 1 (при к == 0):

#include <iostream> 
using namespace std; 


float pow(float base, float ex); 

float nth_root(float A, int n) { 
    const int K = 6; 
    float x[K] = {1}; 
    for (int k = 0; k < K - 1; k++) 
     x[k + 1] = (1.0/n) * ((n - 1) * x[k] + A/pow(x[k], n - 1)); 
    return x[K-1]; 
} 

float pow(float base, float ex){ 
    if (base == 0) 
     return 0; 
    // power of 0 
    if (ex == 0){ 
     return 1; 
    // negative exponenet 
    }else if(ex < 0){ 
     return 1/pow(base, -ex); 
    // fractional exponent 
    }else if (ex > 0 && ex < 1){ 
     return nth_root(base, 1/ex); 
    }else if ((int)ex % 2 == 0){ 
     float half_pow = pow(base, ex/2); 
     return half_pow * half_pow; 
    //integer exponenet 
    }else{ 
     return base * pow(base, ex - 1); 
    } 
} 
int main_pow(int, char **){ 
    for (int ii = 0; ii< 10; ii++){\ 
     cout << "pow(" << ii << ", .5) = " << pow(ii, .5) << endl; 
     cout << "pow(" << ii << ", 2) = " << pow(ii, 2) << endl; 
     cout << "pow(" << ii << ", 3) = " << pow(ii, 3) << endl; 
    } 
    return 0; 
} 

тест:

pow(0, .5) = 0.03125 
pow(0, 2) = 0 
pow(0, 3) = 0 
pow(1, .5) = 1 
pow(1, 2) = 1 
pow(1, 3) = 1 
pow(2, .5) = 1.41421 
pow(2, 2) = 4 
pow(2, 3) = 8 
pow(3, .5) = 1.73205 
pow(3, 2) = 9 
pow(3, 3) = 27 
pow(4, .5) = 2 
pow(4, 2) = 16 
pow(4, 3) = 64 
pow(5, .5) = 2.23607 
pow(5, 2) = 25 
pow(5, 3) = 125 
pow(6, .5) = 2.44949 
pow(6, 2) = 36 
pow(6, 3) = 216 
pow(7, .5) = 2.64575 
pow(7, 2) = 49 
pow(7, 3) = 343 
pow(8, .5) = 2.82843 
pow(8, 2) = 64 
pow(8, 3) = 512 
pow(9, .5) = 3 
pow(9, 2) = 81 
pow(9, 3) = 729 

Answer 3

Я смотрел на этой бумаге here, которая описывает, как аппроксимировать экспоненциальную функцию для двойной точности. После небольшого исследования в Википедии о представлении с плавающей точкой с одинарной точностью я разработал эквивалентные алгоритмы. Они только реализовали функцию ехра, так что я нашел обратную функцию для журнала, а затем просто сделали

POW(a, b) = EXP(LOG(a) * b). 

составление этого gcc4.6.2 дает POW функцию почти в 4 раза быстрее, чем реализации стандартной библиотеки (компиляция с O2) ,

Примечание: код для EXP скопирован почти дословно из прочитанной бумаги и функция LOG скопирована с here.

Вот соответствующий код:

#define EXP_A 184 
    #define EXP_C 16249 

    float EXP(float y) 
    { 
     union 
     { 
     float d; 
     struct 
     { 
    #ifdef LITTLE_ENDIAN 
      short j, i; 
    #else 
      short i, j; 
    #endif 
     } n; 
     } eco; 
     eco.n.i = EXP_A*(y) + (EXP_C); 
     eco.n.j = 0; 
     return eco.d; 
    } 

    float LOG(float y) 
    { 
     int * nTemp = (int*)&y; 
     y = (*nTemp) >> 16; 
     return (y - EXP_C)/EXP_A; 
    } 

    float POW(float b, float p) 
    { 
     return EXP(LOG(b) * p); 
    } 

Существует еще некоторые оптимизации вы можете сделать здесь, или, возможно, это достаточно хорошо. Это грубое приближение, но если бы вы были удовлетворены ошибками, введенными с использованием двойного представления, я предполагаю, что это будет удовлетворительным.

Answer 4

Я думаю, что вы могли бы попытаться решить его, используя серию Тейлора, . http://en.wikipedia.org/wiki/Taylor_series

С серией Тейлора вы можете решить любой трудный для решения расчет, такой как 3^3.8, используя уже известные результаты, такие как 3^4. В этом случае у вас есть 3^4 = 81, поэтому

3^3.8 = 81 + 3.8 * 3 (3.8 - 4) + .. + .. и так далее зависит от того, насколько велика ваша n, вы получите более близкое решение вашей проблемы.

Answer 5

Я и мой друг столкнулись с подобной проблемой, пока мы работаем над проектом OpenGL, а в некоторых случаях недостаточно.У нашего преподавателя также была та же проблема, и он сказал нам, чтобы разделить мощность на целые и плавающие части. Например, если вы хотите рассчитать x^11.5, вы можете рассчитать sqrt (x^115, 10), что может привести к более точному результату.