- samrrr
- Постоялец
Существует 2 основных подхода для эффективной реализации трансцендентных функций:
1) Интерполяция в таблице.
Достоинства метода:
- наиболее высокая точность,
- очень высокая или даже самая высокая производительность.
Недостатки:
- занимает большой объём памяти.
2) Аппроксимация по методу Чебышёва или Паде-Чебышёва.
Достоинства метода:
- высокая точность (лишь ненамного уступающая табличному методу),
- очень высокая производительность (приблизительно соответствующая или превосходящая табличный метод);
- сравнительно компактный код.
Недостатки:
- для построения качественной аппроксимации требуются специальные знания и навыки.
Наряду с методами Чебышёва и Паде-Чебышёва можно ещё упомянуть метод равномерных приближений, но он лишь ненамного превосходит по точности первые два метода, и в нём сложнее определяются коэффициенты.
Метод Чебышёва (аппроксимацию многочленом), как правило, целесообразно применять для реализаций функций с одинарной или более низкой точностью, а Паде-Чебышёва (аппроксимацию дробно-рациональной функцией) - для реализаций функций с двойной или более высокой точностью, а также при медленной сходимости ряда Чебышёва (как, например, в случае арктангенса). Метод Паде-Чебышёва имеет дополнительное преимущество в производительности за счёт возможности параллельного вычисления многочленов в числителе и знаменателе с помощью векторных команд. Наконец, соотношение между точностью и производительностью регулируется за счёт выбора числа слагаемых в частичной сумме ряда Чебышёва.
Приведу два примера реализации тригонометрических функций с одинарной точностью на SSE-инструкциях по методу Чебышёва. Приведённый код допускает дальнейшую оптимизацию - главным образом, за счёт применения FMA-инструкций.
1. Функция, параллельно вычисляющая синус и косинус угла в градусах. Здесь b0, ..., b4 - коэффициенты многочлена, полученного в результате приведения подобных в частичной сумме ряда Чебышёва для функции f(t)=sin(pi/2*t)/t, где |t|<=1. Нужно заметить, что если для вычисления многочленов в этой функции применить схему Горнера, то немного возрастёт точность.
// Параллельное вычисление синуса и косинуса в градусах
void _declspec(naked) _vectorcall SinCosD(float x, float &s, float &c)
{
static const float ct[8] = // Таблица констант
{
1/180.0f, // Множитель для приведения x
1.74532924E-2f, // b0/90
-0.0f, // 80000000h
90.0f, // Константа для перехода от cos к sin
1.34955584E-11f, // b2/90^5
3.91499991E-22f, // b4/90^9
-8.86095693E-7f, // b1/90^3
-9.77248790E-17f // b3/90^7
};
_asm
{
movups xmm3, oword ptr [offset ct] // xmm3 = 90 # 80000000h : b0/90 # 1/180
mulss xmm3, xmm0 // xmm3 = x/180
cvtss2si eax, xmm3 // eax = k - округлённое до целых значение x/180
movhlps xmm2, xmm3 // xmm2 = 90 # 80000000h
imul eax, 180 // eax = 180*k; of=1, если переполнение
jno sc_cont // В случае слишком большого |x| считать, как при x=0
sub eax, eax // Для этого делаем eax = 0
xorps xmm0, xmm0 // и xmm0 = 0
sc_cont: // Продолжаем для корректного значения x
cvtsi2ss xmm1, eax // xmm1 = 180*k
shl eax, 29 // При нечётном k установить знаковый бит eax
subss xmm0, xmm1 // xmm0 = x-k*180 = 90*t - число в диапазоне [-90; 90]
movd xmm1, eax // В xmm1 - знаковая маска результата
orps xmm2, xmm0 // xmm2 = 90 # -|90*t|
movlhps xmm1, xmm1 // Знаковую маску распространить на старшую часть xmm1
haddps xmm2, xmm0 // xmm2 = 90*t : 90-|90*t| - приведённые значения аргументов
movddup xmm4, qword ptr [16+offset ct] // xmm4 = b4/90^9 # b2/90^5 : b4/90^9 # b2/90^5
xorps xmm1, xmm2 // В xmm1 - приведённые значения аргументов с учётом знака
movsldup xmm2, xmm2 // xmm2 = 90*t # 90*t : 90-|90*t| # 90-|90*t| = ts # ts : tc # tc
mulps xmm2, xmm2 // xmm2 = ts^2 # ts^2 : tc^2 # tc^2
movddup xmm0, qword ptr [24+offset ct] // xmm0 = b3/90^7 # b1/90^3 : b3/90^7 # b1/90^3
mulps xmm4, xmm2 // Умножить векторно коэффициенты 2 и 4 на аргументы многочлена
addps xmm0, xmm4 // Векторно прибавить к произведению коэффициенты 1 и 3
shufps xmm3, xmm3, 153 // xmm3 = b0/90 : b0/90 - свободный член
mulps xmm0, xmm2 // xmm0 = b1/90^3 * ts^2 + b2/90^5 * ts^4 : b1/90^3 * tc^2 + b2/90^5 * tc^4
movshdup xmm4, xmm0 // xmm4 = b3/90^7 * ts^2 + b4/90^9 * ts^4 : b3/90^7 * tc^2 + b4/90^9 * tc^4
mulps xmm2, xmm2 // xmm2 = ts^4 : tc^4
addps xmm0, xmm3 // xmm0 = сумма (0,1,2) для синуса : сумма (0,1,2) для косинуса
mulps xmm4, xmm2 // xmm4 = b3/90^7 * ts^6 + b4/90^9 * ts^8 : b3/90^7 * tc^6 + b4/90^9 * tc^8
addps xmm0, xmm4 // xmm0 = сумма для синуса : сумма для косинуса
mulps xmm0, xmm1 // Получаем в xmm0 готовый результат (-1)^k*t*f(t)
movss [edx], xmm0 // Сохранить косинус в переменной c
movhlps xmm0, xmm0 // Переслать старшую часть xmm0 в младшую
movss [ecx], xmm0 // Сохранить синус в переменной s
ret // Вернуться
}
}2. Функция, вычисляющая тангенс угла в градусах. Здесь b0, ..., b4 - коэффициенты многочлена, полученного в результате приведения подобных в частичной сумме ряда Чебышёва для функции f(t)=t/tg(pi/4*t), где |t|<=1.
_declspec(naked) float _vectorcall TgD(float x) // Тангенс в градусах
{
static const float ct[6] = // Таблица констант
{
-0.0111111111f, // -1/90 - множитель для приведения x
57.2957793f, // b0*45
-5.81775793E-3f, // b1/45
-1.18170581E-7f, // b2/45^3
-3.39436425E-12f, // b3/45^5
-1.22505855E-16f // b4/45^7
};
_asm
{
mov edx, offset ct // В edx - адрес таблицы констант
movss xmm1, [edx] // xmm1 = -1/90
mulss xmm1, xmm0 // xmm1 = -x/90
cvtss2si ecx, xmm1 // ecx = -k, где k - округлённое до целых значение x/90
imul ecx, 90 // ecx = -90*k
jno tg_cont // В случае слишком большого |x| считать, как при x=0
sub ecx, ecx // Для этого обнуляем ecx,
xorps xmm0, xmm0 // и обнуляем xmm0
tg_cont: // Продолжаем для корректного значения x
cvtsi2ss xmm2, ecx // xmm2 = -90*k
shl ecx, 30 // При нечётном k установить знаковый бит ecx
addss xmm2, xmm0 // xmm2 = x-90*k = 45*t - в диапазоне [-45; 45]
movd xmm1, ecx // В xmm1 - знаковая маска результата
movss xmm0, [edx+20] // xmm0 = b4/45^7 - инициализировать сумму
xorps xmm1, xmm2 // xmm1=(-1)^k*45*t - выставить правильный знак будущего результата
mulss xmm2, xmm2 // xmm2=(45*t)^2 - подготовить аргумент для вычисления многочлена
mov cl, 4 // ecx = 4 или 80000004h - инициализировать счётчик цикла
tg_loop: // Цикл вычисления многочлена по схеме Горнера
mulss xmm0, xmm2 // Умножить текущую сумму на аргумент многочлена
addss xmm0, [edx+ecx*4] // Прибавить текущий коэффициент
_emit 0x67 // Префикс изменения разрядности для получения команды loopw
loop tg_loop // Учесть все 5 коэффициентов. В xmm0 формируется 45*f(t)
jnz tg_div // Перейти для вычисления -f(t)/t, если k нечётное, т.е. |tg x|>1
movaps xmm2, xmm0 // Для чётного k обменять местами регистры xmm0 и xmm1
movaps xmm0, xmm1 // В этом случае |tg x|<=1, и нужно считать t/f(t)
movaps xmm1, xmm2 // Теперь xmm1 = 45*f(t), xmm0 = 45*t
tg_div: // Получение окончательного результата
divss xmm0, xmm1 // Найти -f(t)/t, если k нечётное, или t/f(t), если k чётное
ret // Вернуться
}
}Aenigma
> Aenigma
А arcsin arccos arcctg2 есть?
Есть готовый код для arctg и arctg2. Я правильно понял, вас эти функции интересуют? Арксинус и арккосинус вычисляются через арктангенс.
Выкладываю код функции арктангенс с одинарной точностью на SSE-инструкциях. Аппроксимация получена методом Паде-Чебышёва. Функции arcsin и arccos напишу позже, сейчас готового кода нет.
_declspec(naked) float _vectorcall Arctg(float x) // Быстрый арктангенс с одинарной точностью
{
// Для |x|<=1 используется формула:
//
// a0 + a0*a1*x^2 + a0*a2*x^4 + a0*a3*x^6
// arctg x = x * ------------------------------------------------------- = x*P/Q
// a0 + a0*(a1+b0)*x^2 + a0*(a2+b1)*x^4 + a0*(a3+b2)*x^6
//
// Константа a0 сокращается, но она нужна для уменьшения погрешности округления остальных констант
// Для |x|>1 используется формула:
//
// P
// arctg x = pi/2*sgn(x)-arctg(1/x) = pi/2*sgn(x) - -----
// x*Q
static const float ct[8] = // Таблица констант
{
1.0f, // 1
1.59940588f, // a0
0.472401142f, // a0*a2
0.939805329f, // a0*(a2+b1)
1.82871687f, // a0*a1
2.36185217f, // a0*(a1+b0)
0.0151078226f, // a0*a3
0.0844737813f // a0*(a3+b2)
};
_asm
{
mov eax, offset ct // В eax - адрес таблицы констант
movaps xmm2, xmm0 // xmm2 = x
movups xmm4, [eax] // xmm4 = a0*(a2+b1) # a0*a2 : a0 # 1
mulss xmm2, xmm2 // xmm2 = x^2
movups xmm3, [eax+16] // xmm3 = a0*(a3+b2) # a0*a3 : a0*(a1+b0) # a0*a1
comiss xmm2, xmm4 // Сравнить |x| с 1, и выставить флаги по результату сравнения
jbe arctg_cont // Перейти, если |x|<=1
divss xmm4, xmm2 // Для x>1 вычислить xmm4 = 1/x^2
movaps xmm2, xmm4 // xmm2 = 1/x^2
arctg_cont: // Продолжить со значением y=x^2 или y=1/x^2 в xmm2
shufps xmm2, xmm2, 0 // xmm2 = y # y : y # y
shufps xmm4, xmm4, 229 // xmm4/a0 = a2+b1 # a2 : 1 # 1
mulps xmm3, xmm2 // xmm3/a0 = (a3+b2)*y # a3*y : (a1+b0)*y # a1*y
addps xmm3, xmm4 // xmm3/a0 = (a3+b2)*y+(a2+b1) # a3*y+a2 : (a1+b0)*y+1 # a1*y+1
mulps xmm2, xmm2 // xmm2 = y^2 # y^2 : y^2 # y^2
movhlps xmm1, xmm3 // xmm1/a0 = (a3+b2)*y+(a2+b1) # a3*y+a2
mulps xmm1, xmm2 // xmm1/a0 = (a3+b2)*y^3+(a2+b1)*y^2 # a3*y^3+a2*y^2
addps xmm1, xmm3 // xmm1/a0 = (a3+b2)*y^3+(a2+b1)*y^2+(a1+b0)*y+1 # a3*y^3+a2*y^2+a1*y+1
movshdup xmm3, xmm1 // xmm1 = P; xmm3 = Q
ja arctg_big // Прейти, если нужно обработать случай |x|>1
divss xmm1, xmm3 // xmm1 = P/Q
mulss xmm0, xmm1 // Для |x|<1 в xmm0 получить x*P/Q
ret // Вернуться
arctg_big: // При |x|>1 находим pi/2*sgn(x)-arctg(1/x)
comiss xmm0, xmm2 // cf=1, если x<0
mov eax, 0x7F921FB6 // В eax - сдвинутое влево двоичное представление числа pi/2
rcr eax, 1 // eax = pi/2*sgn(x)
mulss xmm3, xmm0 // Вычислить делитель xmm3 = Q*x
movd xmm0, eax // xmm0 = pi/2*sgn(x)
divss xmm1, xmm3 // xmm1 = P/(x*Q)
subss xmm0, xmm1 // xmm0 = pi/2*sgn(x)-P/(x*Q)
ret // Вернуться
}
}Aenigma
> Аппроксимация по методу Чебышёва или Паде-Чебышёва.
Так-то максимум точности полинома - это всё-таки алгоритм Ремеза, вроде же?
Но, насколько я понимаю, разложение по базису полиномов Чебышёва - довольно близко по точности.
Вычисление синусов и т. д. - это не только polynomial approximation, но и range reduction, который у тебя, насколько вижу, примитивный.
Оно не всегда важно, но люди заморачиваются:
https://www.csee.umbc.edu/~phatak/645/supl/Ng-ArgReduction.pdf
Для подобных развлечений человечество придумало специальную тулзу:
https://www.sollya.org/
Кроме Горнера есть https://en.wikipedia.org/wiki/Estrin%27s_scheme .
В догонку, мои темы:
http://gamedev.ru/code/forum/?id=204567&page=7&m=4879202#m97
http://gamedev.ru/flame/forum/?id=220141&page=3&m=4360776#m35
Я выше писал, что наилучшей аппроксимацией в определённом семействе функций считается равномерное приближение. Каким методом оно получено, не столь важно. Важно, что вычисление коэффициентов равномерного приближения сложнее, а результат по точности близок к методам Чебышёва или Паде-Чебышёва (в зависимости от того, на каком семействе функций строится аппроксимация).
Вычисление синусов и т. д. - это не только polynomial approximation, но и range reduction, который у тебя, насколько вижу, примитивный.
А разве для аргумента в градусах требуется более сложное редуцирование?
Для радианов можно позаморачиваться, но большого смысла в этом нет, потому что сам формат исходных данных float всё равно сведёт на нет все попытки устранить потерю точности при редуцировании, если |x|>>1.
Метод Эстрина вычисления многочленов у меня используется в коде функций SinCosD и Arctg. Он немного эффективнее параллелится с помощью векторных команд, но точность у него получается хуже, чем у Горнера.
Aenigma
А, cорри, ты вроде алгоритм Ремеза упоминаешь под названием
> метод равномерных приближений
Не, Ремеза я не упоминаю, просто говорю о равномерных приближениях.
Aenigma
> А разве для аргумента в градусах требуется более сложное редуцирование?
Чисто теоретически - да. Для аргумента 2^35=34359738368° твой код, подозреваю, даст ответ, отличный от точного 248°.
На практике, вероятно, пофиг.
В стандартной библиотеке есть https://en.cppreference.com/w/cpp/numeric/math/fmod , который считает точный остаток (не особо быстро).
FordPerfect, я прекрасно понимаю, что такое приведение аргумента тригонометрических функций, и как оно влияет на точность. Особенность моего алгоритма в том, что на вход поступает число float, т.е. число с 7-8 верными значащими десятичными цифрами. А у тебя пример с числом, имеющим 11 значащих цифр. Поэтому применение более продвинутого редуцирования точность особо не повысит, разве что создаст такую иллюзию. А вот производительность снизится.
P.S. Посмотрел твои темы. Вижу, ты тоже глубоко копаешь. :)
Продублирую из соседней темы:
http://marc-b-reynolds.github.io/math/2020/03/11/SinCosPi.html
Присмотрелся, в C23 вроде, собираются добавить вкусностей про плавучку, в частности sinpi.
http://www.open-std.org/jtc1/sc22/wg14/www/docs/n2596.pdf (Annex F)
Aenigma
> поступает число float, т.е. число с 7-8 верными значащими десятичными цифрами. А у тебя пример с числом, имеющим 11 значащих цифр.
Так-то это число точно представимо во float. Понятно, что погрешность в младшем знаке для него превышает период редукции (единица младшего разряда для него равна 4096°). Т. е. с т. з. backward error analysis точность не выглядит оправданной. Но некоторые - заморачиваются и в этих случаях. Кэхэн неоднократно в своих статьях говорил, что бывает осмысленно делать точнее. Насколько его аргументы убедительны - дело личное. Примеры аргументов можно прочитать здесь (они у него чуть по разному озвучиваются из статьи в статью):
https://people.eecs.berkeley.edu/~wkahan/Qdrtcs.pdf
Конкретно здесь - я не призываю никоим образом делать более дотошную редукцию для SinDeg. Оно действительно выглядит не стоящим того.
Aenigma
> Есть готовый код для arctg и arctg2. Я правильно понял, вас эти функции интересуют?
Да, то что нужно.
Aenigma
> Арксинус и арккосинус вычисляются через арктангенс.
arctg(sqrt(1-a*a)/a) что-то неочень...
FordPerfect
Так-то это число точно представимо во float
А откуда мы знаем, что это точно это число? Оно могло случайно, причём с небольшой вероятностью, совпасть с точным представлением в формате float, но в общем случае такие числа появляются в результате округления некой другой мантиссы до 23-х бит. Трудно представить, что кто-то свои исходные данные для вычисления тригонометрических функций хранит с точностью, например, до нескольких десятков или тысяч градусов и ждёт при этом точного результата. Явно абсурдная ситуация.
Поэтому, что касается аргумента в градусах, то в моём случае дополнительных мер по точному приведению аргумента не требуется. Но я согласен, что при расчётах в радианах немного повысить точность по сравнению с самым примитивным редуцированием, нужно. Потому что в этом случае погрешности округления при редуцировании играют более существенную роль.
samrrr,
на досуге подумаю, как лучше написать.
насколько это будет производительнее решения через треугольник?
Aenigma
> очень высокая или даже самая высокая производительность.
Мы живем не в прошлом тысячелетии, поэтому доступ к некешированной памяти — одна из самых дорогих операций на наших компьютерах. Соответственно, часто табличные методы проигрывают аппроксимациям полиномами. На практике имеет смысл совмещать и использовать сплайны невысокого порядка. Размер таблицы сильно меньше и формулы достаточно простые.
> Важно, что вычисление коэффициентов равномерного приближения сложнее, а результат по точности близок к методам Чебышёва или Паде-Чебышёва
Важно, что нахождение коэффициентов аппроксимации — задача однократная, поэтому можно хоть неделю молотить тупым перебором, лишь бы в рантайме бесплатно получить увеличение точности.
Тема в архиве.