Проблема реализации exp из math.h

Question

У меня есть функция-копия exp из math.h. Она работает корректно, но на больших значениях (25<) начинает теряться точность по сравнению с оригиналом. Подскажите, пожалуйста, есть ли какая-нибудь формула или свойства чтобы работать в данном случае с большими значениями?

long double my_exp(double x) {
    if(x < 0){
        return 1.0/my_exp(-x);
    }
    long double sum = 1.0;
    long double term = 1.0;
    for(int i = 1; i < 100; i++) {
        term *= x / i;
        sum += term;
    }
    return sum;
}

Пример при my_exp(25), моя функция возвращает 72004899337.385918 а оригинальный exp(25) 72004899337.385880

my_exp(30) 10686474581524.459685 exp(30) 10686474581524.462891

Answer 1

Как вариант. Причем широко используемый, по крайней мере, когда машины были большими :)

exp(a+b) == exp(a) * exp(b)

exp(n*x) == exp(x) ^ n

Так что считайте точно в каком-то приемлемом диапазоне, а затем приводите исходное значение x к этому диапазону. Например, вы умеете точно считать для x от 0 до 1, а надо для 30? Значит, возводите (быстрым возведением в степень) в 30-ю степень функцию от 1.

Update

Ладно, вот конкретный пример. Я взял аппроксимацию exp(x) именно из тех времен "больших машин"; вы напишите свою, более точную (учтите, что при возведении в степень N относительная точность упадет примерно в N раз...).
Привожу для простоты к диапазону 0-0.5.

//// Быстрое возведение в степень
double qpow(double x, unsigned int e)
{
    double res = 1.;
    for(;e;e>>=1)
    {
        if (e&1) res *= x;
        x *= x;
    }
    return res;
}

double apprexp(double x)  // exp(x), 0 << x << ln2 ~ 0.69
{
    return 1./(((((((-0.0001413161*x+0.0013298820)*x
           -0.0083013598)*x+0.0416573475)*x
           -0.1666653019)*x+0.4999999206)*x
           -0.9999999995)*x+1);
}

double qexp(double x)
{
    bool inv = x < 0;
    if (x < 0) x = -x;

    unsigned int ex = floor(x*2);
    x -= ex/2.;

    x = qpow(apprexp(0.5),ex)*apprexp(x);
    if (inv) x = 1/x;

    return x;
}

int main(int argc, char * argv[])
{
    for(double x = 0.001; x < 50; x *= 2)
    {
        cout << setprecision(10) << x << "\t\t" << exp(x) << "\t\t" << qexp(x) << endl;
    }
}

Вот, примерно так... Что сложного?

БЛОГ НА HUSL

Проблема реализации exp из math.h

Ответы (1 шт):