Comments 26
x87 fpu со стеком это legacy, которое не рекомендуется к использованию. Современные компиляторы генерируют sse-инструкции. В них регистры не организованы в стек, а доступны непосредственно, по номерам. Современные процессоры с их миллиардами транзисторов реализуют out-of-order superscalar pipeline, на который стек ложится плохо.
Код на x87 скалярный, 80-bit регистры загружаются по невыровненному адресу за несколько тактов. Единственное возможное преимущество - повышенная точность, но ее можно достичь другими методами (например см метод Кэхена). Если оставаться в рамках стандартных 64-bit, то можно задействовать SSE, AVX, AVX-512. Это легко перекроет выгоду от x87 даже с учетом более медленного алгоритма.
А ещё доводилось реализовывать алгоритм Блюстейна для дискретного преобразования Фурье произвольного размера. И внезапно оказалось, что он крайне критичен к точности и на размере порядка 50000 семплов не дотягивает даже до half precision, не говоря уже о single. Простое переписывание на FPU с сохранением промежуточных вычислений с 80-битной точностью позволило вернуть близкую к double точность.
Ну таких опций компилятора наверное еще не придумали. Я лично этот метод руками писал для нахождения суммы разностей квадратов на fp32 (вот пример на SSE см SquaredDifferenceKahanSum32f ).
Зависит от задачи. Но если суммировать порядка 10^6 чисел, то легко получить погрешность во 2-3 знаке. Причем результат будет различасться в зависимости от того, что используешь: Scalar/SSE/AVX/AVX-512 (см подробнее здесь). Впрочем чем больше длина вектора, тем погрешность меньше.
x87 fpu со стеком это legacy
FPU станет legacy только тогда, когда появится нативная поддержка 128-битной точности чисел с плавающей точкой — тогда использовать 80 бит действительно не станет никакого смысла, особенно в случае выравнивания данных по 16 байтам. Не зря же FPU в 64-битном режиме никуда не делся — хотя очевидно, что с поддержкой старых программ это никак не связано, потому как они в 64-битном режиме в любом случае работать не будут.
И проблема стековой организации FPU вовсе не в том, что вычисления не распараллеливаются и оттого выполняются чуточку медленнее — а в том, что на него нельзя навесить интрисики, как это сделано для MMX/SSE/AVX. Поэтому если в ЯП нет его нативной поддержки — то приходится программировать на ассемблере, что далеко не всем нравится.
Алгоритм Кэхена тут тоже особо не поможет, потому что FPU помимо 20 знаков точности обеспечивает ещё и больший диапазон значений (10^4932 против 10^308 у double). Кроме того, алгоритм Кэхена определён только для сложения — умножать им не получится, а погрешности после умножения нивелируют всю накопленную Кэхеном точность — проверял лично на методе наименьших квадратов.
Intel optimization manual рекомендует не использовать x87, если у вас нет требований к точности, которых можно добиться только с 80-битной арифметикой. Только это я хотел сказать, когда писал, что это legacy. Понятно, что если вам для вашей задачи нужна повышенная точность, то придётся использовать x87. Я просто призываю не использовать эти инструкции бездумно, когда без них можно обойтись.
Кроме того, вычисления всё чаще выгружают на gpu, где 80-битная арифметика не поддерживается в принципе. Если в этом случае нужна повышенная точность, приходится что-то изобретать. Может быть, эти подходы и на CPU применимы?
Часть про стек и out-of-order superscalar относилась к окончанию статьи, в котором автор упоминает про несколько стеков и переключение между ними.
С повышенной точностью есть один нюанс — её необходимость далеко не всегда очевидна, и возникает обычно просто при масштабировании задач. Да даже если просто синус через ряд считать — не получится точность double получить, и посчитать его на FPU будет проще, чем лукап-таблицы городить. Синус, понятно, тут для примера — функция элементарная и готовых реализаций её достаточно, о чём в руководстве Intel и говорится. А вот если вам потребуется что-то менее распространённое типа функции ошибок — тут без FPU не обойтись чтобы получить просто обычную, а не расширенную точность (без использования лукап-таблиц и прочих хитростей). Например, вот в этой реализации используется switch на 100 ответвлений — такое себе удовольствие, и не каждый компилятор такое автоматически оптимизирует.
Красивый код.
static double erfcx_y100(double y100)
{
switch ((int) y100) {
case 0: {
double t = 2*y100 - 1;
return 0.70878032454106438663e-3 + (0.71234091047026302958e-3 + (0.35779077297597742384e-5 + (0.17403143962587937815e-7 + (0.81710660047307788845e-10 + (0.36885022360434957634e-12 + 0.15917038551111111111e-14 * t) * t) * t) * t) * t) * t;
}
...
}
// we only get here if y = 1, i.e. |x| < 4*eps, in which case
// erfcx is within 1e-15 of 1..
return 1.0;
}LLVM 10 (-O3 -march=skylake) превращает switch со 100 ветками в переход по таблице:
_ZL10erfcx_y100d: # @_ZL10erfcx_y100d
.cfi_startproc
# %bb.0:
vcvttsd2si eax, xmm0
cmp eax, 99
ja .LBB12_1
# %bb.2:
jmp qword ptr [8*rax + .LJTI12_0]
.LBB12_3:
vaddsd xmm0, xmm0, xmm0
vaddsd xmm0, xmm0, qword ptr [rip + .LCPI12_793]
vmulsd xmm1, xmm0, qword ptr [rip + .LCPI12_794]
vaddsd xmm1, xmm1, qword ptr [rip + .LCPI12_795]
vmulsd xmm1, xmm0, xmm1
vaddsd xmm1, xmm1, qword ptr [rip + .LCPI12_796]
vmulsd xmm1, xmm0, xmm1
vaddsd xmm1, xmm1, qword ptr [rip + .LCPI12_797]
vmulsd xmm1, xmm0, xmm1
vaddsd xmm1, xmm1, qword ptr [rip + .LCPI12_798]
vmulsd xmm1, xmm0, xmm1
vaddsd xmm1, xmm1, qword ptr [rip + .LCPI12_799]
vmulsd xmm0, xmm0, xmm1
vaddsd xmm0, xmm0, qword ptr [rip + .LCPI12_800]
ret
...
.section .rodata,"a",@progbits
.p2align 3
.LJTI12_0:
.quad .LBB12_3
.quad .LBB12_4
...Кстати, если уж мы говорим о точности, мне удивительно, что в этом коде не используются инструкции FMA.
Если уж смотреть на эту конкретную реализацию, то не будет ли эффективнее вместо switch на 100 веток написать таблицу коэффициентов? Как я вижу, во всех ветках выполняется одно и то же действие, только числа разные. Ну и использовать инструкции fma --- ошибка будет медленнее накапливаться.
Сейчас ещё проверил: если подать -ffast-math, компилятор начинает генерировать FMA, но ну его к чёрту. Нельзя включать -ffast-math, если есть требования к точности:
.LBB12_3:
vmovsd xmm1, qword ptr [rip + .LCPI12_1] # xmm1 = mem[0],zero
vfmadd213sd xmm1, xmm0, qword ptr [rip + .LCPI12_794] # xmm1 = (xmm0 * xmm1) + mem
vmovsd xmm0, qword ptr [rip + .LCPI12_795] # xmm0 = mem[0],zero
vfmadd213sd xmm0, xmm1, qword ptr [rip + .LCPI12_796] # xmm0 = (xmm1 * xmm0) + mem
vfmadd213sd xmm0, xmm1, qword ptr [rip + .LCPI12_797] # xmm0 = (xmm1 * xmm0) + mem
vfmadd213sd xmm0, xmm1, qword ptr [rip + .LCPI12_798] # xmm0 = (xmm1 * xmm0) + mem
vfmadd213sd xmm0, xmm1, qword ptr [rip + .LCPI12_799] # xmm0 = (xmm1 * xmm0) + mem
vfmadd213sd xmm0, xmm1, qword ptr [rip + .LCPI12_800] # xmm0 = (xmm1 * xmm0) + mem
vfmadd213sd xmm0, xmm1, qword ptr [rip + .LCPI12_801] # xmm0 = (xmm1 * xmm0) + mem
ret
а триногометрические функции Вы тоже будете командами SSE считать? И уж сразу полиномами Чебышева?
Вообще, да. Intel 64 and IA-32 Architectures Optimization Reference Manual именно это и рекомендует, если не нужна 80-битная точность. Intel заявляет, что программная реализация на SSE будет быстрее.
User/Source Coding Rule 15. (M impact, ML generality) Usually, math libraries take advantage of the transcendental instructions (for example, FSIN) when evaluating elementary functions. If there is no critical need to evaluate the transcendental functions using the extended precision of 80 bits, applications should consider an alternate, software-based approach, such as a look-up-table-based algorithm using interpolation techniques. It is possible to improve transcendental performance with these techniques by choosing the desired numeric precision and the size of the look-up table, and by taking advantage of the parallelism of the SSE and the SSE2 instructions.
Шикарный совет не пользоваться FPU! И точность 80 бит здесь ни при чем. Попробуйте составить таблицы хотя бы для точности 32 бита. Как говорится, счастливого пути!
Конечно, если в игрушках какие-нибудь текстуры нужно быстро чуть повернуть, можно упрощенный расчет применять. Но не все сводится к текстурам...
Собственно, это уже реализовано в стандартной библиотеке языка C. В частности, реализация синуса в glibc: https://github.com/lattera/glibc/blob/master/sysdeps/ieee754/dbl-64/s_sin.c
А почему вы не используете LLVM чтобы не воевать с бинарным кодом самостоятельно?
Подобные вопросы один персонаж еще более 200 лет назад задавал: к чему учить географию, если извозчики сами куда надо довезут?
Например, как через LLVM задать начальное значение стека в программе? То самое, которое в PE-заголовке называется size_of_stack_reserve?
Я полагаю, эту задачу решает линкер, который не является частью LLVM в общем случае. Не знаю, как оно на Windows, на Linux есть 3 опции: GNU ld, GNU gold и LLD.
Например, в PL/1 глубина стека задается в заголовке главной программы. Т.е. это принадлежность самого языка, а не редактора связей. Но это ладно. В представлении многих LLVM - это такая волшебная палочка в руках волшебников, которая все превращает в идеальный код. Но я не случайно спросил о стеке, потому что, похоже, в LLVM нет такого понятия вообще. Т.е. виртуальный процессор LLVM не очень похож на настоящие, как и у байт-кода, к примеру.
Интересно, отчего нельзя было проводить большие вычисления на том компьютере, где уже был сопроцессор?
Как увеличить стек FPU