Всем привет!
Возникла потребность написать быстрое вычисление Sin и Cos. За основу для вычислений взял разложение по ряду Тейлора. Использую в 3D-системах (OpenGL и графическая библиотека своей разработки). К сожалению свести ряд «идеально» для Double не получается, но это компенсируется хорошим ускорением. Код написан на встроенном в Delphi XE6 ассемблере. Используется SSE2.
Для научных вычислений не подходит, а для использования в играх вполне.
Точности хватает, чтобы покрыть разрядность числа Single, которое используется
для умножения на матрицу.
В итоге:
- Достигнутая точность результата равна: 10.e-13
- Максимальное расхождение с CPU — 0.000000000000045.
- Скорость увеличена в сравнении с CPU в 4.75 раза.
- Скорость увеличена в сравнении с Math.Sin и Math.Cos в 2.6 раза.
Для теста использовал процессор Intel Core-i7 6950X Extreme 3.0 ГГц.
Исходный текст на Delphi встроен в комментарии к ассемблеру.
Исходный код:
Заголовок спойлера
var
gSinTab: array of Double;
gSinAddr: UInt64;
const
AbsMask64: UInt64 = ($7FFFFFFFFFFFFFFF);
FSinTab: array[0..7] of Double =
(1.0 / 6.0, //3!
1.0 / 120.0, //5!
1.0 / 5040.0, //7!
1.0 / 362880.0, //9!
1.0 / 39916800.0, //11!
1.0 / 6227020800.0, //13!
1.0 / 1307674368000.0, //15!
1.0 / 355687428096000.0); //17!
//Максимумы функции Sin для положительных значений
QSinMaxP: array[0..3] of Double = (0.0, 1.0, 0.0, -1.0);
//Максимумы функции Sin для отрицательных значений
QSinMaxM: array[0..3] of Double = (0.0, -1.0, 0.0, 1.0);
//Знаки квадрантов для положительных значений
QSinSignP: array[0..3] of Double = (1.0, 1.0, -1.0, -1.0);
//Знаки квадрантов для отрицательных значений
QSinSignM: array[0..3] of Double = (-1.0, -1.0, 1.0, 1.0);
function Sinf(const A: Double): Double;
asm
.NOFRAME
// На входе A в xmm0
// Четверть окружности
// X := IntFMod(Abs(A), PI90, D);
// Result := FNum - (Trunc(FNum / FDen) * FDen);
pxor xmm3, xmm3
comisd A, xmm3
jnz @ANZ
ret // Result := 0.0;
@ANZ:
//Флаг отрицательного угла
//Minus := (A < 0.0);
setc cl
movsd xmm1, [AbsMask64] //Abs(A)
pand A, xmm1
movsd xmm1, [PI90] //PI90 = FDen
movsd xmm2, A //Копия A = FNum
divsd xmm2, xmm1 //(FNum / FDen)
roundsd xmm2, xmm2, 11b //11b - Trunc
cvtsd2si rax, xmm2 //Целая часть (X в EAX)
mulsd xmm2, xmm1
subsd xmm0, xmm2
//-----------------------------------
//Нормализация квадранта
and rax, 3 // D := (D and 3);
//-----------------------------------
// Проверка на максимум функции
// if (X = 0.0) then
// begin
// if Minus then
// Exit(QSinMaxM[D]) else
// Exit(QSinMaxP[D]);
// end;
comisd xmm0, xmm3
jnz @XNZ
shl rax, 3 //умножить номер квадранта на размер Double (8 байт)
cmp cl, 1 //Если угол отрицательынй, то переход
je @MaxMinus
@MaxPlus:
lea rdx, qword ptr [QSinMaxP]
movsd xmm0, qword ptr [rdx + rax]
ret
@MaxMinus:
lea rdx, qword ptr [QSinMaxM]
movsd xmm0, qword ptr [rdx + rax]
ret
//-----------------------------------
@XNZ:
//При нечетном квадранте нужно вычесть градусы для симметрии
// if (D and 1) <> 0 then X := (PI90 - X);
mov edx, eax
and edx, 1
cmp edx, 0
je @DZ
subsd xmm1, xmm0
movsd xmm0, xmm1
//-----------------------------------
@DZ:
// Result остается в xmm0
// X в xmm0
// N := (X * X) в xmm2
// F := (N * X) в xmm1
// N
movsd xmm2, xmm0 // Копия X
mulsd xmm2, xmm2 // N := (X * X)
// F
movsd xmm1, xmm2 // Копия N
mulsd xmm1, xmm0 // F := (N * X)
//Загружаем таблицу факториалов для синуса
mov rdx, [gSinTab]
movaps xmm4, dqword ptr [rdx]
movaps xmm6, dqword ptr [rdx + 16]
movaps xmm8, dqword ptr [rdx + 32]
movaps xmm10, dqword ptr [rdx + 48]
//Извлекаем нечетную часть таблицы
movhlps xmm5, xmm4
movhlps xmm7, xmm6
movhlps xmm9, xmm8
movhlps xmm11, xmm10
// Result := X - F * PDouble(gSinAddr)^;
mulsd xmm4, xmm1 //FSinTab[0]
subsd xmm0, xmm4
// F := (F * N);
// Result := Result + F * PDouble(gSinAddr + 8)^;
mulsd xmm1, xmm2
mulsd xmm5, xmm1 //FSinTab[1]
addsd xmm0, xmm5
// F := (F * N);
// Result := Result - F * PDouble(gSinAddr + 16)^;
mulsd xmm1, xmm2
mulsd xmm6, xmm1 //FSinTab[2]
subsd xmm0, xmm6
// F := (F * N);
// Result := Result + F * PDouble(gSinAddr + 24)^;
mulsd xmm1, xmm2
mulsd xmm7, xmm1 //FSinTab[3]
addsd xmm0, xmm7
// F := (F * N);
// Result := Result - F * PDouble(gSinAddr + 32)^;
mulsd xmm1, xmm2
mulsd xmm8, xmm1 //FSinTab[4]
subsd xmm0, xmm8
// F := (F * N);
// Result := Result + F * PDouble(gSinAddr + 40)^;
mulsd xmm1, xmm2
mulsd xmm9, xmm1 //FSinTab[5]
addsd xmm0, xmm9
// F := (F * N);
// Result := Result - F * PDouble(gSinAddr + 48)^;
mulsd xmm1, xmm2
mulsd xmm10, xmm1 //FSinTab[6]
subsd xmm0, xmm10
// F := (F * N);
// Result := Result + F * PDouble(gSinAddr + 56)^;
mulsd xmm1, xmm2
mulsd xmm11, xmm1 //FSinTab[7]
addsd xmm0, xmm11
//-----------------------------------
// if Minus then
// Result := (Result * QSinSignM[D]) else
// Result := (Result * QSinSignP[D]);
shl rax, 3 //Умножаем на 8
cmp cl, 1
je @Minus
lea rdx, qword ptr [QSinSignP]
mulsd xmm0, qword ptr [rdx + rax]
ret
@Minus:
lea rdx, qword ptr [QSinSignM]
mulsd xmm0, qword ptr [rdx + rax]
end;
//------------------------------------
function Cosf(const A: Double): Double;
asm
.NOFRAME
// A в xmm0
// Четверть окружности
// X := IntFMod(AMod, PI90, D);
// Result := FNum - (Trunc(FNum / FDen) * FDen);
pxor xmm3, xmm3
comisd A, xmm3
jnz @ANZ
movsd xmm0, [DoubleOne]
ret // Result := 1.0;
//-----------------------------------
@ANZ:
//Флаг отрицательного угла
//Minus := (A < 0.0);
setc cl
movsd xmm1, [AbsMask64] //Abs(A)
pand A, xmm1
//-----------------------------------
movsd xmm1, [PI90] //PI90 = FDen
//-----------------------------------
// if Minus then
// AMod := Abs(A) - PI90 else
// AMod := Abs(A) + PI90;
cmp cl, 1
je @SubPI90
addsd A, xmm1
jmp @IntFMod
@SubPI90:
subsd A, xmm1
//-----------------------------------
@IntFMod:
movsd xmm2, A //Копия A = FNum
divsd xmm2, xmm1 //(FNum / FDen)
roundsd xmm2, xmm2, 11b //11b - Trunc
cvtsd2si rax, xmm2 //Целая часть (X в EAX)
mulsd xmm2, xmm1
subsd xmm0, xmm2
//-----------------------------------
//Нормализация квадранта
and rax, 3 // D := (D and 3);
//-----------------------------------
// Проверка на максимум функции
// if (X = 0.0) then
// begin
// if Minus then
// Exit(QSinMaxM[D]) else
// Exit(QSinMaxP[D]);
// end;
comisd xmm0, xmm3
jnz @XNZ
shl rax, 3 //умножить номер квадранта на размер Double (8 байт)
cmp cl, 1 //Если угол отрицательынй, то переход
je @MaxMinus
@MaxPlus:
lea rdx, qword ptr [QSinMaxP]
movsd xmm0, qword ptr [rdx + rax]
ret
@MaxMinus:
lea rdx, qword ptr [QSinMaxM]
movsd xmm0, qword ptr [rdx + rax]
ret
//-----------------------------------
@XNZ:
//При нечетном квадранте нужно вычесть градусы для симметрии
// if (D and 1) <> 0 then X := (PI90 - X);
mov edx, eax
and edx, 1
cmp edx, 0
je @DZ
subsd xmm1, xmm0
movsd xmm0, xmm1
@DZ:
// Result остается в xmm0
// X в xmm0
// N := (X * X) в xmm2
// F := (N * X) в xmm1
// N
movsd xmm2, xmm0 // Копия X
mulsd xmm2, xmm2 // N := (X * X)
// F
movsd xmm1, xmm2 // Копия N
mulsd xmm1, xmm0 // F := (N * X)
//Загружаем таблицу факториалов для синуса
mov rdx, [gSinTab]
movaps xmm4, dqword ptr [rdx]
movaps xmm6, dqword ptr [rdx + 16]
movaps xmm8, dqword ptr [rdx + 32]
movaps xmm10, dqword ptr [rdx + 48]
//Извлекаем нечетную часть таблицы
movhlps xmm5, xmm4
movhlps xmm7, xmm6
movhlps xmm9, xmm8
movhlps xmm11, xmm10
// Result := X - F * PDouble(gSinAddr)^;
mulsd xmm4, xmm1 //FSinTab[0]
subsd xmm0, xmm4
// F := (F * N);
// Result := Result + F * PDouble(gSinAddr + 8)^;
mulsd xmm1, xmm2
mulsd xmm5, xmm1 //FSinTab[1]
addsd xmm0, xmm5
// F := (F * N);
// Result := Result - F * PDouble(gSinAddr + 16)^;
mulsd xmm1, xmm2
mulsd xmm6, xmm1 //FSinTab[2]
subsd xmm0, xmm6
// F := (F * N);
// Result := Result + F * PDouble(gSinAddr + 24)^;
mulsd xmm1, xmm2
mulsd xmm7, xmm1 //FSinTab[3]
addsd xmm0, xmm7
// F := (F * N);
// Result := Result - F * PDouble(gSinAddr + 32)^;
mulsd xmm1, xmm2
mulsd xmm8, xmm1 //FSinTab[4]
subsd xmm0, xmm8
// F := (F * N);
// Result := Result + F * PDouble(gSinAddr + 40)^;
mulsd xmm1, xmm2
mulsd xmm9, xmm1 //FSinTab[5]
addsd xmm0, xmm9
// F := (F * N);
// Result := Result - F * PDouble(gSinAddr + 48)^;
mulsd xmm1, xmm2
mulsd xmm10, xmm1 //FSinTab[6]
subsd xmm0, xmm10
// F := (F * N);
// Result := Result + F * PDouble(gSinAddr + 56)^;
mulsd xmm1, xmm2
mulsd xmm11, xmm1 //FSinTab[7]
addsd xmm0, xmm11
//-----------------------------------
// if Minus then
// Result := (Result * QSinSignM[D]) else
// Result := (Result * QSinSignP[D]);
shl rax, 3 //Умножаем на 8
cmp cl, 1
je @Minus
lea rdx, qword ptr [QSinSignP]
mulsd xmm0, qword ptr [rdx + rax]
ret
@Minus:
lea rdx, qword ptr [QSinSignM]
mulsd xmm0, qword ptr [rdx + rax]
end;
Initialization
//Выровненный массив
SetLength(gSinTab, 8);
//Адрес таблицы
gSinAddr := UInt64(@FSinTab[0]);
//Копируем таблицу в выровненный массив
Move(FSinTab[0], gSinTab[0], SizeOf(Double) * 8);
Finalization
SetLength(gSinTab, 0);
Конструктивные предложения и замечания приветствуются.
Автор: тащит всю команду
Предлагаю быструю версию функции для Delphi, которая параллельно вычисляет синус и косинус в градусах с одинарной точностью. Алгоритм основан на разложении функции sin(x*2/pi)/x в ряд Чебышёва с последующим подбором коэффициентов. Функция обеспечивает погрешность, не превышающую примерно 2 ULP для |x|<=2^24. Для бОльших значений аргумента начинает снижаться точность вплоть до полной её потери. Функция написана под 32-битную архитектуру x86 с использованием команд SSE3, т.е. подходит практически для любых машин. Код приведён ниже.
// Параллельное вычисление синуса и косинуса в градусах с одинарной точностью
procedure SinCosD(x: Single; var s,c: Single);
const
ct: array[0..7] of Single = ( // Таблица констант, где b0, …, b4 – коэф-ты аппрокс. многочлена
5.55555569E-3, // Множитель 1/180 для приведения x
1.74532924E-2, // b0/90
-0.0, // 80000000h
90, // Константа для перехода от cos к sin
1.34955597E-11, // b2/90^5
3.90824210E-22, // b4/90^9
-8.86095734E-7, // b1/90^3
-9.77194052E-17); // b3/90^7
asm
mov ecx, offset ct // В ecx – адрес таблицы констант
movss xmm0, [x] // В xmm0 поместить аргумент x
movups xmm3, [ecx] // xmm3 = 90 # 80000000h : b0/90 # 1/180
mulss xmm3, xmm0 // xmm3 = x/180
cvtss2si ebp, xmm3 // ebp = k – округлённое до целых значение x/180
movhlps xmm2, xmm3 // xmm2 = 90 # 80000000h
imul ebp, 180 // ebp = 180*k; of=1, если переполнение
jno @@sc_cont // В случае слишком большого |x| результат, как при x=0
sub ebp, ebp // Для этого делаем ebp = 0
xorps xmm0, xmm0 // xmm0 = 0
@@sc_cont: // Продолжаем для корректного значения x
cvtsi2ss xmm1, ebp // xmm1 = 180*k
shl ebp, 29 // При нечётном k установить знаковый бит ebp
subss xmm0, xmm1 // xmm0 = x-k*180 = 90*t – число в диапазоне [-90; 90]
movd xmm1, ebp // В xmm1 – знаковая маска результата
orps xmm2, xmm0 // xmm2 = 90 # -|90*t|
movlhps xmm1, xmm1 // Знаковую маску распространить на старшую часть xmm1
haddps xmm2, xmm0 // xmm2 = 90*t : 90-|90*t| – приведённые значения аргументов
movddup xmm4, [ecx+16] // 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, [ecx+24] // xmm0 = b3/90^7 # b1/90^3 : b3/90^7 # b1/90^3
mulps xmm4, xmm2 // Умножить векторно коэффициенты 2 и 4 на аргументы многочлена
addps xmm0, xmm4 // Векторно прибавить к произведению коэффициенты 1 и 3
mulps xmm0, xmm2 // xmm0 = b1/90^3 * ts^2 + b2/90^5 * ts^4 : b1/90^3 * tc^2 + b2/90^5 * tc^4
shufps xmm3, xmm3, 153 // xmm3 = b0/90 : b0/90 – свободный член
mulps xmm2, xmm2 // xmm2 = ts^4 : 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 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 – суммы для синуса и для косинуса, кроме свободных членов
addps xmm0, xmm3 // xmm0 = сумма для синуса : сумма для косинуса
mulps xmm0, xmm1 // Получаем в xmm0 готовый результат (-1)^k*t*f(t)
movss [edx], xmm0 // Сохранить косинус в переменной c
movhlps xmm0, xmm0 // Переслать старшую часть xmm0 в младшую
movss [eax], xmm0 // Сохранить синус в переменной s
end;