Я исправил код в Wikibooks и добавил несколько дополнительных комментариев (ответ Шута хороший), так что теперь он собирается и работает правильно (протестировано с GDB, одноступенчатое с layout ret / tui reg float). Это разница между ревизиями .Ревизия, которая представила ошибку fmul st1,st1 недопустимая инструкция , находится здесь , но даже до этого ей не удалось очистить стек x87, когда это было сделано.
Просто ради интереса, яхотел написать более эффективную версию, которая загружает a и b только один раз.
и которая обеспечивает больший параллелизм на уровне команд, выполняя все , а не , включая результат cosпервый.то есть подготовьте 2*a*b перед тем, как умножить это на cos(ang), чтобы оба этих вычисления могли выполняться параллельно.Предполагая, что fcos является критическим путем, моя версия имеет только одну fmul и одну fsubp задержку от fcos результата до fsqrt ввода.
default rel ; in case we assemble this in 64-bit mode, use RIP-relative addressing
... declare stuff, omitted.
fld qword [a] ;load a into st0
fld st0 ; st1 = a because we'll need it again later.
fmul st0, st0 ;st0 = a * a = a^2
fld qword [b] ;load b into st0 (pushing the a^2 result up to st1)
fmul st2, st0 ; st2 = a*b
fmul st0, st0 ;st0 = b^2, st1 = a^2, st2 = a*b
faddp ;st0 = a^2 + b^2 st1 = a*b; st2 empty
fxch st1 ;st0 = a*b st1 = a^2 + b^2 ; could avoid this, but only by using cos(ang) earlier, worse for critical path latency
fadd st0,st0 ;st0 = 2*a*b st1 = a^2 + b^2
fld qword [ang]
fcos ;st0 = cos(ang) st1 = 2*a*b st2 = a^2+b^2
fmulp ;st0=cos(ang)*2*a*b st1 = a^2+b^2
fsubp st1, st0 ;st0 = (a^2 + b^2) - (2 * a * b * cos(ang))
fsqrt ;take square root of st0 = c
fstp qword [c] ;store st0 in c and pop, leaving the x87 stack empty again ‒ and we're done!
OfКонечно, x87 в значительной степени устарел.На современном x86 обычно вы используете скаляр SSE2 (или упакованный!) Для чего-либо с плавающей запятой.
x87 имеет две вещи для современного x86: аппаратная точность 80 бит (против64-битный double), и это хорошо для небольшого размера кода (байты машинного кода, а не количество инструкций или размер исходного кода).Хорошие кэши команд обычно означают, что размер кода не является достаточно важным фактором, чтобы сделать x87 оправданным для производительности кода FP, потому что он, как правило, медленнее, чем SSE2, из-за дополнительных инструкций, связанных с неуклюжим стеком x87.
И для новичков, или из-за размера кода, x87 имеет трансцендентные функции, такие как fcos и fsin, и log / exp встроены в виде отдельных инструкций.Они микрокодируются многими мопами, и, вероятно, не быстрее, чем функция скалярной библиотеки, но на некоторых процессорах вы можете согласиться с компромиссом между скоростью и точностью, который они делают, и абсолютной скоростью.По крайней мере, если вы используете x87 в первую очередь, в противном случае вы должны отразить результаты в / из регистров XMM с сохранением / перезагрузкой.
Уменьшение диапазона для sin / cos не делает никакого расширенного-точный материал, чтобы избежать огромных относительных ошибок, очень близких к кратному Пи, просто используя внутреннее 80-битное (64-битное значение и) значение Пи.(Реализация библиотеки может или не может сделать это, в зависимости от желаемого соотношения скорости и точности.) См. Intel недооценивает границы ошибок на 1,3 квинтиллиона .
(И, конечно, x87 в 32-битовый код обеспечивает совместимость с Pentium III и другими процессорами, у которых не было SSE2 для двойного, только SSE1 для плавающих или вообще не было регистров XMM. x86-64 имеет SSE2 в качестве базового уровня, поэтому это преимущество не существует в x86-64.)
Для новичков огромным недостатком x87 является отслеживание регистров стека x87, а не накопление содержимого.Вы можете легко получить код, который работает один раз, но затем выдает NaN, когда вы помещаете его в цикл, потому что вы не уравновешивали свои операции стека x87.
extern cos
global cosine_law_sse2_scalar
cosine_law_sse2_scalar:
movsd xmm0, [ang]
call cos ; xmm0 = cos(ang). Avoid using this right away so OoO exec can do the rest of the work in parallel
movsd xmm1, [a]
movsd xmm2, [b]
movaps xmm3, xmm1 ; copying registers should always copy the full reg, not movsd merging into the old value.
mulsd xmm3, xmm2 ; xmm3 = a*b
mulsd xmm1, xmm1 ; a^2
mulsd xmm2, xmm2 ; b^2
addsd xmm3, xmm3 ; 2*a*b
addsd xmm1, xmm2 ; a^2 + b^2
mulsd xmm3, xmm0 ; 2*a*b*cos(ang)
subsd xmm1, xmm3 ; (a^2 + b^2) - 2*a*b*cos(ang)
sqrtsd xmm0, xmm3 ; sqrt(that), in xmm0 as a return value
ret
;; This has the work interleaved more than necessary for most CPUs to find the parallelism
В этой версии только 11 мопов послеcall cos возвращается.(https://agner.org/optimize/). Это довольно компактно и довольно просто. Не отслеживает стек x87. И имеет те же цепочки зависимостей, что и x87, не используя результат cos, пока у нас уже не будет 2*a*b.
Мы можем даже поиграть с загрузкой a и b вместе, как одним 128-битным вектором, но затем распаковать его, чтобы сделать разные вещи с двумя половинками, или получить b^2 из верхнего элемента какскалярный, неуклюжий. Если бы SSE3 haddpd был только 1 моп, это было бы здорово (и давайте сделаем a*b + a*b и a^2 + b^2 с одной инструкцией, учитывая правильные входные данные), но на всех процессорах, у которых это есть, это 3микрооперации.
(PS против PD имеет значение только для фактических математических инструкций, таких как MULSS / SD. Для перестановок FP и регистровых копий, просто используйте любую инструкцию FP, которая доставит ваши данные туда, куда вы хотите, с предпочтением для PS / SS, потому что они корочекодировки в машинном коде. Вот почему я использовал movaps; movapd - это всегда пропущенная оптимизация, тратящая впустую 1 байт, если только вы не выполняете инструкции с целью выравнивания.)
;; I didn't actually end up using SSE3 for movddup or haddpd, it turned out I couldn't save uops that way.
global cosine_law_sse3_less_shuffle
cosine_law_sse3_less_shuffle:
;; 10 uops after the call cos, if both extract_high_half operations use pshufd or let movhlps have a false dependency
;; or if we had AVX for vunpckhpd xmm3, xmm1,xmm1
;; and those 10 are a mix of shuffle and MUL/ADD.
movsd xmm0, [ang]
call cos ; xmm0 = cos(ang). Avoid using this right away so OoO exec can do the rest of the work in parallel
movups xmm1, [a] ; {a, b} (they were in contiguous memory in this order. low element = a)
movaps xmm3, xmm1
; xorps xmm3, xmm3 ; break false dependency by zeroing. (xorps+movhlps is maybe better than movaps + unpckhpd, at least on SnB but maybe not Bulldozer / Ryzen)
; movhlps xmm3, xmm1 ; xmm3 = b
; pshufd xmm3, xmm1, 0b01001110 ; xmm3 = {b, a} ; bypass delay on Nehalem, but fine on most others
mulsd xmm3, [b] ; xmm3 = a*b ; reloading b is maybe cheaper than shufling it out of the high half of xmm1
addsd xmm3, xmm3 ; 2*b*a
mulsd xmm3, xmm0 ; 2*b*a*cos(ang)
mulpd xmm1, xmm1 ; {a^2, b^2}
;xorps xmm2, xmm2 ; we don't want to just use xmm0 here; that would couple this dependency chain to the slow cos(ang) critical path sooner.
movhlps xmm2, xmm1
addsd xmm1, xmm2 ; a^2 + b^2
subsd xmm1, xmm3 ; (a^2 + b^2) - 2*a*b*cos(ang)
sqrtsd xmm0, xmm1 ; sqrt(that), in xmm0 as a return value
ret
И мы можемсделайте еще лучше с AVX, сохранив копию регистра MOVAPS, потому что 3-операндные неразрушающие версии инструкций VEX позволяют нам поместить результат в новый регистр, не разрушая ни один из входов.Это действительно здорово для FP shuffle, потому что в SSE * нет опций копирования и перемешивания для операндов FP, только pshufd, что может вызвать дополнительную задержку обхода на некоторых процессорах.Таким образом, он сохраняет MOVAPS и (закомментированные) XORPS, которые нарушают зависимость от того, что произвело старое значение XMM2 для MOVHLPS.(MOVHLPS заменяет младшие 64 бита места назначения старшими 64 битами src, поэтому он имеет входную зависимость от обоих регистров).
global cosine_law_avx
cosine_law_avx:
;; 9 uops after the call cos. Reloading [b] is good here instead of shuffling it, saving total uops / instructions
vmovsd xmm0, [ang]
call cos ; xmm0 = cos(ang). Avoid using this right away so OoO exec can do the rest of the work in parallel
vmovups xmm1, [a] ; {a, b} (they were in contiguous memory in this order. low element = a)
vmulsd xmm3, xmm1, [b] ; xmm3 = a*b
vaddsd xmm3, xmm3 ; 2*b*a. (really vaddsd xmm3,xmm3,xmm3 but NASM lets us shorten when dst=src1)
vmulsd xmm3, xmm0 ; 2*b*a*cos(ang)
vmulpd xmm1, xmm1 ; {a^2, b^2}
vunpckhpd xmm2, xmm1,xmm1 ; xmm2 = { b^2, b^2 }
vaddsd xmm1, xmm2 ; a^2 + b^2
vsubsd xmm1, xmm3 ; (a^2 + b^2) - 2*a*b*cos(ang)
vsqrtsd xmm0, xmm1,xmm1 ; sqrt(that), in xmm0 as a return value. (Avoiding an output dependency on xmm0, even though it was an ancestor in the dep chain. Maybe lets the CPU free that physical reg sooner)
ret
Я тестировал только первую версию x87, поэтому возможновозможно, пропустил шаг в одном из других.