Наиболее эффективный способ реализации функции pow() в плавающей точке

Question

Я пытаюсь реализовать мою собственную версию функций pow() и sqrt(), так как моя пользовательская библиотека не имеет поддержки pow()/sqrt() с плавающей точкой.

Кто-нибудь может помочь?

Answer 1

Конечно - это легко, если у вас есть экспоненциальные и натуральные лог-функции.

поскольку y = x^nВы можете взять натуральный логарифм обеих сторон:

ln(y) = n*ln(x)

Тогда взятие экспоненты обеих сторон дает вам то, что вы хотите:

y = exp(n*ln(x))

Если вы хотите чего-то лучшего, лучшее место, которое я знаю, это Абрамовиц и Стегун.

Answer 2

Да, Солнце может (Oracle сейчас, я думаю):

fdlibm, "свободно распространяемая математическая библиотека", имеет sqrt и pow вместе со многими другими математическими функциями.

Это довольно высокотехнологичные реализации, и, конечно, ничто не является "самой эффективной" реализацией чего-то подобного. Вы после исходного кода, чтобы сделать это, или вы действительно не так много ищете pow а также sqrt, но на самом деле ищете образование в программировании алгоритмов с плавающей точкой?

Answer 3

Обратите внимание, что если в вашем наборе инструкций есть инструкция для квадратного корня или мощности, вам будет намного лучше использовать это. Например, в инструкциях с плавающей запятой x87 есть инструкция fsqrtи дополнения SSE2 включают другую инструкцию sqrtsd, что, вероятно, будет намного быстрее, чем большинство решений, написанных на C. На самом деле, по крайней мере, gcc использует две инструкции, когда компиляция происходит на компьютере с архитектурой x86.

Однако для власти все становится несколько мутно. В наборе команд с плавающей запятой x87 есть инструкция, которую можно использовать для вычисления n*log2(n), а именно fyl2x, Еще одна инструкция, fldl2e, хранит log2(e) в стеке с плавающей запятой. Возможно, вы захотите взглянуть на них.

Вы также можете взглянуть на то, как отдельные библиотеки C делают это. dietlibc, например, просто использует fsqrt:

sqrt:
    fldl 4(%esp)
    fsqrt
    ret

glibc использует реализацию Sun для машин, где аппаратная инструкция квадратного корня недоступна (в разделе sysdeps/ieee754/flt-32/e-sqrtf.c) и использует fsqrt в наборе команд x86 (хотя gcc можно указать вместо этого использовать sqrtsd инструкция).

Answer 4

Квадратный корень правильно реализован с помощью итерационного метода Ньютона.

Answer 5

double ipow(int base, int exp)
{

bool flag=0;
if(exp<0) {flag=1;exp*=-1;}
int result = 1;
while (exp)
{
    if (exp & 1)
        result *= base;
    exp >>= 1;
    base *= base;
}
if(flag==0)
return result;
else
return (1.0/result);
}
//most suitable way to implement power function for integer to power integer

Answer 6

Для вычисления квадратного корня с плавающей точкой в ​​C я бы рекомендовал использовать fsqrt если вы нацелены на x86. Вы можете использовать такую ​​инструкцию ASM с:

asm("fsqrt" : "+t"(myfloat));

Для GCC или

asm {

fstp myfloat

fsqrt

fldp myfloat

}

Или что-то подобное для Visual Studio.

Для реализации pow следует использовать большой оператор switch, такой как на http://www.upitasoft.com/link/powLUT.txt. Это может вызвать некоторые проблемы с кешем, но если вы сохраните его таким, чтобы это не было проблемой, просто ограничьте диапазон (обратите внимание, вы все равно можете оптимизировать код, который я предоставил)

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

float result = exp(number * log(power));

Но обычно это медленно и / или неточно.

Надеюсь, я помог.

Answer 7

Самый быстрый способ, которым я могу придумать выполнение pow(), будет выглядеть следующим образом (обратите внимание, это довольно сложно):

//raise x^y
double pow(double x, int y) {
    int power;
    map<int, double> powers;
    for (power = 1; power < y; power *= 2, x *= x)
        powers.insert(power, x);

    while (power > y) {
        //figure out how to get there
        map<int, double>::iterator p = powers.lower_bound(power - y);
        //p is an iterator that points to the biggest power we have that doesn't go over power - y
        power -= p->first;
        x /= p->second;
    }

    return x;
}

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

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

Редактировать 2: Хе-хе, мой плохой. power *= 2 вместо power++