Быстрая фиксированная точка pow, log, exp и sqrt

Question

У меня есть класс с фиксированной точкой (10.22), и мне нужна функция pow, sqrt, exp и log.

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

Я предполагаю, что, как только у меня будет функция exp, будет довольно легко реализовать pow и sqrt, так как онипросто станьте.

pow( x, y ) => exp( y * log( x ) )
sqrt( x )   => pow( x, 0.5 )

Это просто те функции exp и log, которые я нахожу трудными (как будто я помню некоторые из моих правил ведения журнала, я больше не могу вспомнить о них).

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

Обратите внимание: это HASбыть кроссплатформенным и в чистом коде C / C ++, поэтому я не могу использовать оптимизацию ассемблера.

Answer 1

Очень простое решение - использовать приличное табличное приближение. На самом деле вам не нужно много данных, если вы правильно уменьшите свои данные. exp(a)==exp(a/2)*exp(a/2), что означает, что вам действительно нужно рассчитать exp(x) для 1 < x < 2. В этом диапазоне приближение Рунга-Кутты дало бы разумные результаты с ~ 16 записями IIRC.

Аналогично, sqrt(a) == 2 * sqrt(a/4) == sqrt(4*a) / 2, что означает, что вам нужны только записи в таблице для 1 < a < 4. Log (a) немного сложнее: log(a) == 1 + log(a/e). Это довольно медленная итерация, но log (1024) составляет всего 6,9, поэтому у вас не будет много итераций.

Вы бы использовали аналогичный алгоритм "integer-first" для pow: pow(x,y)==pow(x, floor(y)) * pow(x, frac(y)). Это работает, потому что pow(double, int) тривиально (разделяй и властвуй).

[править] Для интегрального компонента log(a) может быть полезно сохранить таблицу 1, e, e^2, e^3, e^4, e^5, e^6, e^7, чтобы можно было уменьшить log(a) == n + log(a/e^n) с помощью простого жестко заданного двоичного поиска в этой таблице. Улучшение с 7 до 3 шагов не так уж и велико, но это означает, что вам нужно делить только один раз на e^n вместо n раз на e.

[править 2] И для этого последнего log(a/e^n) термина вы можете использовать log(a/e^n) = log((a/e^n)^8)/8 - каждая итерация выдает еще 3 бита при поиске в таблице . Это сохраняет ваш код и размер таблицы маленьким. Обычно это код для встроенных систем, и у них нет больших кешей.

[править 3] Это все еще не умный на моей стороне. log(a) = log(2) + log(a/2). Вы можете просто сохранить значение с фиксированной точкой log2=0.30102999566, сосчитать число ведущих нулей, сместить a в диапазон, используемый для вашей таблицы поиска, и умножить это смещение (целое число) на постоянную с фиксированной точкой log2. Может быть всего 3 инструкции.

Использование e для шага сокращения дает просто «хорошую» константу log(e)=1.0, но это ложная оптимизация. 0,30102999566 является такой же хорошей константой, как 1,0; оба являются 32-битными константами в фиксированной точке 10,22. Использование 2 в качестве константы для уменьшения диапазона позволяет использовать битовый сдвиг для деления.

Вы все еще получаете трюк от правки 2, log(a/2^n) = log((a/2^n)^8)/8. По сути, это дает вам результат (a + b/8 + c/64 + d/512) * 0.30102999566 - с b, c, d в диапазоне [0,7]. a.bcd действительно восьмеричное число. Не удивительно, так как мы использовали 8 в качестве силы. (Трюк одинаково хорошо работает с мощностью 2, 4 или 16.)

[править 4] Все еще был открытый конец. pow(x, frac(y) это просто pow(sqrt(x), 2 * frac(y)), и у нас есть приличное 1/sqrt(x). Это дает нам гораздо более эффективный подход. Скажите frac(y)=0.101 двоичный, то есть 1/2 плюс 1/8. Тогда это означает, что x^0.101 это (x^1/2 * x^1/8). Но x^1/2 это просто sqrt(x) и x^1/8 это (sqrt(sqrt(sqrt(x))). Сохраняя еще одну операцию, Ньютон-Рафсон NR(x) дает нам 1/sqrt(x), поэтому мы вычисляем 1.0/(NR(x)*NR((NR(NR(x))). Мы только инвертируем конечный результат, напрямую не используем функцию sqrt.

Answer 2

Ниже приведен пример реализации C алгоритма Clay S. Turner с логарифмической базой 2 с фиксированной запятой [ 1 ].Алгоритм не требует никакой таблицы поиска.Это может быть полезно в системах, где ограничения памяти ограничены, а в процессоре отсутствует FPU, как, например, в случае многих микроконтроллеров.База журналов e и база журналов 10 также поддерживаются с помощью свойства логарифмов, которое для любой базы n :

          log (x)
             y
log (x) = _______
   n      log (n)
             y

где, для этого алгоритма, y равно 2.

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

При использовании этих функций ожидается, что x будет значением с фиксированной запятой, масштабированным в соответствии с указанным precision.Например, если precision равно 16, то x должно быть масштабировано до 2 ^ 16 (65536).Результатом является значение с фиксированной запятой с тем же масштабным коэффициентом, что и для ввода.Возвращаемое значение INT32_MIN представляет отрицательную бесконечность.Возвращаемое значение INT32_MAX указывает на ошибку, а errno будет установлено на EINVAL, указывая, что точность ввода была недопустимой.

#include <errno.h>
#include <stddef.h>

#include "log2fix.h"

#define INV_LOG2_E_Q1DOT31  UINT64_C(0x58b90bfc) // Inverse log base 2 of e
#define INV_LOG2_10_Q1DOT31 UINT64_C(0x268826a1) // Inverse log base 2 of 10

int32_t log2fix (uint32_t x, size_t precision)
{
    int32_t b = 1U << (precision - 1);
    int32_t y = 0;

    if (precision < 1 || precision > 31) {
        errno = EINVAL;
        return INT32_MAX; // indicates an error
    }

    if (x == 0) {
        return INT32_MIN; // represents negative infinity
    }

    while (x < 1U << precision) {
        x <<= 1;
        y -= 1U << precision;
    }

    while (x >= 2U << precision) {
        x >>= 1;
        y += 1U << precision;
    }

    uint64_t z = x;

    for (size_t i = 0; i < precision; i++) {
        z = z * z >> precision;
        if (z >= 2U << precision) {
            z >>= 1;
            y += b;
        }
        b >>= 1;
    }

    return y;
}

int32_t logfix (uint32_t x, size_t precision)
{
    uint64_t t;

    t = log2fix(x, precision) * INV_LOG2_E_Q1DOT31;

    return t >> 31;
}

int32_t log10fix (uint32_t x, size_t precision)
{
    uint64_t t;

    t = log2fix(x, precision) * INV_LOG2_10_Q1DOT31;

    return t >> 31;
}

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

[ 1 ] CS Turner, «Быстрый алгоритм двоичного логарифма» , Обработка сигналов IEEE Mag. , с. 124, 140, сентябрь 2010.

Answer 3

Хорошей отправной точкой является Книга Джека Креншоу, "Математический инструментарий для программирования в реальном времени" . Здесь хорошо обсуждаются алгоритмы и реализации для различных трансцендентных функций.

Answer 4

Проверьте мою реализацию sqrt с фиксированной точкой, используя только целочисленные операции.Было весело изобретать.Довольно старый сейчас.

https://groups.google.com/forum/?hl=fr%05aacf5997b615c37&fromgroups#!topic/comp.lang.c/IpwKbw0MAxw/discussion

В противном случае проверьте набор алгоритмов CORDIC .Это способ реализации всех перечисленных вами функций и тригонометрических функций.

РЕДАКТИРОВАТЬ: Я опубликовал проверенный источник на GitHub здесь