Быстрая фиксированная точка pow, log, exp и sqrt
У меня есть класс с фиксированной запятой (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, поэтому любые указатели на этом фронте были бы оценены, даже если бы просто использовать методы, которые я описал выше.
Обратите внимание: это ДОЛЖНО быть кроссплатформенным и в чистом коде C/C++, поэтому я не могу использовать какие-либо оптимизации на ассемблере.
Ответы
Ответ 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, поэтому у вас не будет много итераций.
Вы бы использовали аналогичный алгоритм "сначала целое число" для 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)
простым простым двоичным поиском a в этой таблице. Улучшение с 7 до 3 шагов не так уж и велико, но это означает, что вам нужно делить только один раз на e^n
вместо n
раз на e
.
[edit 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
, сосчитать количество 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)
1/8 (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 напрямую.
Ответ 2
Ниже приведен пример реализации алгоритма C с логической базой 2 с фиксированной запятой Clay S. Turner [ 1 ]. Алгоритм не требует никакой таблицы поиска. Это может быть полезно в системах, где ограничения памяти ограничены, а в процессоре отсутствует FPU, как, например, в случае многих микроконтроллеров. База журналов e и база журналов 10 также поддерживаются с помощью свойства логарифмов, которое для любой базы n:
log (x)
y
log (x) = _______
n log (n)
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 г.
Ответ 3
Хорошей отправной точкой является Книга Джека Креншоу, "Математический инструментарий для программирования в реальном времени" . Он хорошо обсуждает алгоритмы и реализации для различных трансцендентных функций.
Ответ 4
Проверьте реализацию sqrt с фиксированной точкой, используя только целые операции.
Было интересно придумать. Довольно старый.
https://groups.google.com/forum/?hl=fr%05aacf5997b615c37&fromgroups#!topic/comp.lang.c/IpwKbw0MAxw/discussion
В противном случае установите CORDIC набор алгоритмов. Это способ реализовать все перечисленные вами функции и тригонометрические функции.
EDIT: Я опубликовал обзорный источник на GitHub здесь