самостоятельно сделал pow () c ++

Question

Я читал Как я могу самостоятельно написать степенную функцию? , и ответ, данный dan04, привлек мое внимание главным образом потому, что я не уверен в ответе, данном Фортраном, но я принял его и выполнил это:

#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)), но я не уверен в том, чтобы поставить это, поскольку я не уверен, где это поставить, или если я даже думаю об этом праве.

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

РЕДАКТИРОВАТЬ: кто-нибудь знает реализацию с плавающей запятой для дробных показателей? (Я видел двойную реализацию, но это было использование трюка с регистрами, и мне нужна плавающая точка, и добавление библиотеки просто для выполнения трюка, мне было бы лучше просто включить библиотеку математики)

Answer 1

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

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

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

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

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

    #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 2

Я думаю, что вы можете попытаться решить эту проблему, используя серию Тейлора, проверьте это.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) + .. + .. и т. Д. В зависимости от того, насколько велик ваш пболее близкое решение вашей проблемы.

Answer 3

Я думаю, что алгоритм, который вы ищете, может быть ' nth root '.С начальным предположением 1 (для k == 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 4

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