Комментарии 30
Дизассемблировав несколько программ, использующих косинус, я не смог найти ни одной, задействующей команду fcos.
Вспоминается анекдот "А Вы куда звонили, сэр?"
Встроенные в язык sin, cos или, например, арктангенс, для процессоров x86 естественно реализуются встроенными функциями FPU. Именно для этого они в FPU и вводились. Другое дело, что для игрушек часто достаточно очень приближенных значений.
В FPU есть также замечательная функция sincos, которая сразу считает и sin и cos, быстрее, чем отдельными расчетами. При этом sin и cos часто идут такими парами, например, в комплексных числах.
В используемом нами PL/1-KT из-за этого даже добавлена встроенная функция sincos и она активно используется, например, при расчете орбитальных данных.
А автору статьи (который сам признается, что не математик) остается пожелать познакомиться с полиномами Чебышёва, о которых он, похоже, не слышал.
Современные процессоры поддерживают инструкции FPU 8087 для совместимости, но в новых программах они не используются. Отчасти это связано с повсеместным более-менее буквальным внедрением стандарта IEEE 754, не предусматривающего 80-битных чисел, как в 8087, зато налагающего очень строгие требования на совпадение результатов вычислений с математически предсказываемыми.
Насколько мне известно, синус и косинус в современных компиляторах для x86_64 вычисляются библиотечными процедурами на основе наборов команд SSE*.
Что касается содержания статьи, то исходно ясно, что для небольшой точности достаточно таблицы.
Команды SSE* предназначены для параллельного выполнения одной операции с разными наборами операндов и непосредственно отношения к алгоритмам вычисления триногометрических функций не имеют. Насколько мне известно, в современных FPU триногометрические функции вычисляются полиномами Чебышёва, а квадратный корень - по методу Ньютона.
Я глубоко не вникал, но думаю, что вычисление полинома Чебышёва векторизуется и в таком качестве обрабатывается. Ну и SSE* просто включает набор операций вещественной арифметики, ничто на худой конец не мешает их выполнять в том числе и над векторами из единственного элемента.
Вот, к примеру, что происходит в gfortran:
real (8) :: x, y
read *, x
y = sin (x)
print *, y
call __gfortran_st_read
movq %rsp, %rsi
movl $8, %edx
movq %rbx, %rdi
call __gfortran_transfer_real
movq %rbx, %rdi
call __gfortran_st_read_done
vmovsd (%rsp), %xmm0
call _sin
movq lC2(%rip), %rax
movq %rbx, %rdi
movq %rbp, 24(%rsp)
vmovsd %xmm0, 8(%rsp)
movq %rax, 16(%rsp)
movl $7, 32(%rsp)
call __gfortran_st_writeВызывается некая подпрограмма _sin, обменивающаяся значениями через SSEшный регистр XMM0.
А у меня в системной библиотеке вот так
PUBLIC COSD: ;ВЫЗЫВАЕТСЯ ИЗ PL/1
;---- ПЕРЕВОД ГРАДУСОВ В РАДИАНЫ ----
MOV RAX,[RBX]
FILD16 ГРАДУС
FLDPI
JC @
;---- ДЛЯ FLOAT-32 ----
FMUL32 [RAX]
FDIVR
JMPS CS1
;---- ДЛЯ FLOAT-64 ----
@: FMUL64 [RAX]
FDIVR
JMPS CS2
PUBLIC COS: ;ВЫЗЫВАЕТСЯ ИЗ PL/1
MOV RAX,[RBX] ;АДРЕС АРГУМЕНТА
JC @
;---- ДЛЯ FLOAT-32 ----
FLD32 [RAX] ;ЗАГРУЗИЛИ АРГУМЕНТ
CS1: FCOS
POP RAX
LEA RSP,[RSP]-8
FST32P [RSP]
JMP RAX
;---- ДЛЯ FLOAT-64 ----
@: FLD64 [RAX] ;ЗАГРУЗИЛИ АРГУМЕНТ
CS2: FCOS
POP RAX
LEA RSP,[RSP]-8
FST64P [RSP]
JMP RAX
DSEG
ГРАДУС DW 180
думаю, что вычисление полинома Чебышёва векторизуется и в таком качестве обрабатываетсяНет, полиномы Чебышёва используются в аналитическом виде для расчёта коэффициентов многочлена. Подробнее об этом можно почитать тут, там же есть сравнение производительности.
SSE это ещё и скалярные операции. И прочитайте наконец интеловское руководство по оптимизации. Там чёрным по белому написано: не используйте инструкции x87 FPU, они медленные, и сохранены только для совместимости со старым софтом.
Если Вы имеете ввиду "Intel® 64 and IA-32 Architectures
Optimization Reference Manual", то там не говорится, что "сохранены только для совместимости". Там, например утверждается, что
"Использование скалярного кода SSE/SSE2 вместо кода x87 может обеспечить значительное повышение производительности при использовании микроархитектуры Intel NetBurst и скромное преимущество при использовании процессоров Intel Core Solo и Intel Core Duo."
А кроме этого, изначально речь шла о триногометрических функциях, которых нет в SSE. Единственная команда FCOS во время выполнения не обращается к памяти, что сейчас очень ускоряет работу. Кстати, во время ее выполнения можно заниматься, например, целыми вычислениями.
Использование 80 бит в FPU вместо 64 в SSE2 замедляет работу, но зато намного корректней результат округляется до 64 бит.
И если в игрушках этого не требуется, то во многих задачах требования к точности сосвсем другие. Например, один градус отклонения от надира на высоте орбиты МКС - это 7 км на поверхности Земли. В игрушке наплевать, что монстр чуть-чуть не там, где был бы "в натуре". А в реальных задачах может и не наплевать.
Поэтому вычисления FPU более точные, а для триногометрических (и трансцендентных) функций и более быстрые.
Почему-то у меня в документе написано другое: "Scalar floating-point SIMD instructions have lower latencies than equivalent x87 instructions. Scalar SIMD floating-point multiply instruction may be pipelined, while x87 multiply instruction is not. Although x87 supports transcendental instructions, software library implementation of transcendental function can be faster in many cases."
Если только один синус — то конечно хуже, надо же данные сначала отдать в GPU, а затем забрать, это всё дорогостоящие накладные расходы. Если миллион синусов — то конечно лучше. Ну и есть нюанс, что не все GPU поддерживают double, хотя накостылять повышенную точность через 2 single при желании можно.
CUDA предназначена для векторно-конвейерных вычислений с большим объёмом одинаковых операций. На одиночной операции инициализация займёт гораздо дольше времени, чем сама операция.
Вспоминается анекдот "А Вы куда звонили, сэр?"
Не гуглится :(
Возможно вот этот: https://www.anekdot.ru/id/-321919003/
Ваша спираль эпилепсии, - как раз отличный случай применить CORDIC, но не в версии "цифра за цифрой", а в версии "посчитать синусы и косинусы последовательных близких углов".
Честно посчитав (co)sin(x0=0) используйте формулу
sin (x+a) = sin(x) + a*cos(x)
cos (x+a) = cos(x) - a*sin(x)Здесь **a**-постоянное приращение угла, то есть, малая константа, в идеале - степень двойки, собственно, именно тогда умножение превращается в элегантный сдвиг.
То есть, для каждого следующего угла сохраняйте и используйте в расчётах синус и косинус предыдущего.
Трэш начинается уже в ряду Тейлора. Уж на что я - не математик и то - знаю о симметрии и периодичности тригонометрических функций. Дальше читать не стал.
Фактически синус и косинус в стандартной библиотеке Linux в настоящее время считаются так.
Как видно, область определения функций делится на множество интервалов, и в каждом из них используется своё приближение.
Ну уж коль была попытка найти таблицы, то это хорошо забытое старое... Таблицы Брадиса.
И да, вычисления можно свести к целым числам.
А вот интересно, есть ли достаточно эффективный (но не табличный) способ вычисления синуса и косинуса для чисел с фиксированной точкой (т.е. обычный целочисленный тип, где INT_MAX принят за 1)? Для простоты пусть будет даже обычный int8.
Интересно в том числе на аппаратном уровне: наверное можно составить таблицу истинности, по ней провести оптимизацию через карты Карно и попытаться вывести битовую формулу, из которой уже можно составить схему некоего синусатора на логических элементах, преобразующего int8 значение угла в int8 значение синуса.
На Zx Spectrum кто-то делал, но не факт, что без таблиц. На современных процессорах вряд ли можно преимущество по скорости получить, единственный смысл — когда точность больше 15 цифр требуется.
Мне в этом вопросе интересна скорее сама математическая суть данного процесса. Как будет выглядеть логическая схема для синуса, будет ли она простой и элегантной или бессистемным нагромождением и мешаниной элементов, будут ли в ней какие-то закономерности, позволяющие ее масштабировать на большие разрядности и т.д.
Ряд Тейлора из статьи чем плох? пара умножений и деление на число кратное 2. Деление на 2! это сдвиг побитовый числа вправо например.
Как вариант промежуточный вариант с упрощенным рядом Тейлора из одного члена и таблицей из нескольких значений. Или линейной / квадратичной апроксимацией между значениями таблицы. Можно подобрать баланс между размером таблицы и сложностью вычислений.
У меня для микроконтроллеров была идея сделать 2-3 таблицы. Первая, например 10 абсолютных значений функции. Потом 10 значений первой производной. Потом 10 значений второй производной. И результат соответственно
Y = K1 + dX*K2 + dХ^2*K3
где K значения из таблиц, dX расстояние между узлом записанным в таблицу и значением на входе. Если они совпадают, то dx=0 и значение берется сразу из таблицы первой. Плюс вручную можно оптимизировать таблицу, например K3 обнулить в начале синусоиды, нам небольшой изгиб и хватает линейной интерполяции, синус почти как прямая линия. Тогда и лишнее умножение не потребуется.
На Спектруме обычно делают табличку из 256 байт для четверти периода.
Если автор использовал RISC-V, то возможно стоило использовать функции из библиотеки NMSIS DSP: float32_t riscv_sin_f32(float32_t x), q15_t riscv_sin_q15(q15_t x), q31_t riscv_sin_q31(q31_t x) ?
Интересно подсчитать синус и косинус сразу для 4-х членов регистра xmm, команды SINPS, COSPS. Тут надо потоковое вычисление делать, а значит никаких if'ов. Я это частично сделал, но пока расчёт ведётся в пределах нормализированного угла -Pi*2 .. Pi*2, а это высокая погрешность на углах ближе к Pi*2 при малом количестве итераций, и высокая общая погрешность и время выполнения для большом количестве итераций.
В общем, надо найти формулу которая нормализирует угол до -Pi .. Pi либо даже до -Pi/2 .. Pi/2. Тогда итераций можно сделать немного и будет меньше погрешностей при вычислении float. Да, float сам наводит погрешность при большом количестве итераций. И всё это используя только SSE, SSE2, SSE3.
Реализуем с нуля функцию косинуса на языке C