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

Question

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

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

Answer 1

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

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

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

Answer 2

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

Начиная с y = x^n, вы можете взять натуральный журнал с обеих сторон:

ln(y) = n*ln(x)

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

y = exp(n*ln(x))

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

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, такой как upitasoft.com / link / powLUT.h . Это может вызвать некоторые проблемы с кешем, но если вы сохраните его таким, чтобы это не было проблемой, просто ограничьте диапазон (обратите внимание, вы все равно можете оптимизировать код, который я предоставил).

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

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;
}

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

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

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