Работая в компании Gaijin несколько лет назад, мне довелось поучаствовать в портировании пары игр компании на консоль Nintendo Switch, тогда вовсю завоевывающую новые рынки. Для меня это стало первым крупным проектом на этой платформе. А с учетом, что ни команда, ни разработчик движка с платформой, системой сборки и вообще экосистемой Нинтендо знакомы не были, то все грабли приходилось искать и бережно на них наступать. Чтобы опробовать возможности новой платформы, параллельно с портированием игры, был написан внутренний middleware (связка dagor engine + nxsdk + jam) и код обрастал всевозможными тестами, build matrix, бенчмарками, прогоном стабильности и другими внутренними проверками. Надо отметить что на момент 2018 года, в самом switch sdk не было реализовано часть posix функций вроде poll и send/receive, и большая часть функций для работы с файлами, posix прослойку нужно было писать самим. Дошли тогда руки и до написания различных бенчмарков для функций стандартной библиотеки, и были замечены некоторые аномалии в поведении части тригонометрических функций в различных режимах сборки. Для справки, sdk использует урезанный вариант musl libc (https://www.musl-libc.org/), все статически линкуется в один большой бинарник clang'ом от Нинтендо 9 версии (2018 год), который потом запускается на консоли. Доступа к исходникам самой libc в исполнении Нинтендо у нас не было, но всегда можно посмотреть дизасм и боле менее представить что происходит.
Свои действия опишу в настоящем времени, просто помните что на дворе был 2018. Код бенчмарка абсолютно детский, гоняем синус (в сдк вся тригонометрия лежала в tr1) по кругу на консоли
static void nxsdk_sin(benchmark::State& state) { // Code inside this loop is measured repeatedly float i = 0; for( auto _ : state) { float p = tr1::sin(p); i += 0.1f; benchmark::DoNotOptimize(p); } } BENCHMARK(nxsdk_sin);
Получаем такие результаты (меньшее время - лучше)
-----------(меньше лучше)----------------------------- Benchmark Time CPU Iterations ------------------------------------------------------ nxsdk_sin 8.68 ns 8.58 ns 74666667 <<< https://git.musl-libc.org/cgit/musl/tree/src/math/sin.c nxsdk_sin_vec 8.35 ns 8.30 ns 75825003 <<< sin(vec4f) dagor_sin 5.48 ns 5.47 ns 102504001 glibc_sin 8.28 ns 8.08 ns 77623654 <<< https://github.com/lattera/glibc/blob/master/sysdeps/ieee754/dbl-64/s_sin.c
Пока все ожидаемо, формулы для вычисления sin в musl примерно те же, что и в glibc, и времена показывают примерно одинаковые. Но если мы включим -03 + -math:precise, синус от Нинтендо стал почти на треть быстрее, тут либо компилятор слишком умный, либо бенчмарк врет и луна не в той фазе, либо что-то еще. Собственно, когда мой коллега показал мне эти результаты, я ему так и ответил - мол, у тебя бенч сломался, иди перепроверь, и забыл про этот момент. На следующий день мы смотрели этот тест уже вместе, потому что перфу взяться там действительно неоткуда (ну это я сначала так думал). Ожидаемо видно небольшое снижение времени glibc_sin за счет агрессивных оптимизаций clang'a.
-----------(меньше лучше)----------------------------- Benchmark Time CPU Iterations ------------------------------------------------------ nxsdk_sin 6.03 ns 6.00 ns 87625024 nxsdk_sin_vec 5.82 ns 5.78 ns 90022043 dagor_sin 5.38 ns 5.27 ns 104602002 glibc_sin 8.08 ns 8.05 ns 78214356
Еще более "другой" результат, получился когда тест был запущен с флагами -O3 + -math:fast, тут я действительно удивился потому, что быстрее dagor_sin на тот момент не было при сопоставимой точности, а результаты тестов на правильность все синусы проходили одинаково хорошо, величина расхождения между реализациями с std::sin не превышала 1e-6
-----------(меньше лучше)----------------------------- Benchmark Time CPU Iterations ------------------------------------------------------ nxsdk_sin 4.03 ns 3.99 ns 132307692 nxsdk_sin_vec 4.09 ns 4.08 ns 131403251 dagor_sin 5.13 ns 5.10 ns 110400021 glibc_sin 7.43 ns 7.16 ns 79544722
И вот тут пришлось лезть в дизасм, чтобы понять чего такого гении из Нинтендо придумали чтобы добиться таких впечатляющих цифр в бенчмарке.
Начнем с режима -01 (speed), я выложу только важные участки, чтобы было понимание что происходит, и постараюсь не выкладывать весь дизасм, мало ли там какое IP есть. Функция sin из поставки nx sdk, является алиасом для __sin_nx_502, (2018 год, актуальный sdk 4.0 и 5.0) видимо имя функции генерится на основе этих данных
__sin_nx_502
__sin_nx_502(double): sub sp, sp, #32 stp x29, x30, [sp, #16] add x29, sp, #16 fmov x8, d0 mov w9, #8699 ubfx x8, x8, #32, #31 movk w9, #16361, lsl #16 cmp w8, w9 b.hi .ERIT1_3 // .ERITX_X -> .BBX_X arm lsr w9, w8, #20 cmp w9, #996 b.hi .ERIT1_5 mov x9, #4066750463515557888 mov x10, #5147614374084476928 fmov d1, x9 fmov d2, x10 fmul d1, d0, d1 fadd d2, d0, d2 cmp w8, #256, lsl #12 fcsel d1, d1, d2, lo str d1, [sp] ldp x29, x30, [sp, #16] add sp, sp, #32 ret .ERIT1_3: lsr w8, w8, #20 cmp w8, #2047 b.lo .ERIT1_6 fsub d0, d0, d0 ldp x29, x30, [sp, #16] add sp, sp, #32 ret .ERIT1_5: adrp x8, .ERIT1_0 ... неинтересный код add sp, sp, #32 ret .ERIT1_6: mov x0, sp bl __rem_pio2_nx_502(double, double*) // вызов __rem_pio2 and w8, w0, #0x3 ... неинтересный код ldp x29, x30, [sp, #16] add sp, sp, #32 ret .ERIT1_10: ldp d1, d0, [sp] ... fadd d2, d2, d3 fmov d5, #0.50000000 ldr d3, [x8, :lo12:.ERIT1_5] ... ldp x29, x30, [sp, #16] add sp, sp, #32 ret .ERIT1_11: ldp d0, d1, [sp] bl __cos_nx_502(double, double) // вызов __сos ldp x29, x30, [sp, #16] add sp, sp, #32 ret .ERIT1_12: ldp d0, d1, [sp] bl __cos_nx_502(double, double) // вызов __сos fneg d0, d0 ldp x29, x30, [sp, #16] add sp, sp, #32 ret
быстро пробежимся по структуре кода, и онa почти полностью повторяет стандартную функцию из musl - видно сигнатуру вызова __rem_pio2/__cos
double sin(double x) { double y[2]; uint32_t ix; unsigned n; /* High word of x. */ GET_HIGH_WORD(ix, x); ix &= 0x7fffffff; /* |x| ~< pi/4 */ if (ix <= 0x3fe921fb) { if (ix < 0x3e500000) { /* |x| < 2**-26 */ /* raise inexact if x != 0 and underflow if subnormal*/ FORCE_EVAL(ix < 0x00100000 ? x/0x1p120f : x+0x1p120f); return x; } return __sin(x, 0.0, 0); } /* sin(Inf or NaN) is NaN */ if (ix >= 0x7ff00000) return x - x; /* argument reduction needed */ n = __rem_pio2(x, y); switch (n&3) { case 0: return __sin(y[0], y[1], 1); case 1: return __cos(y[0], y[1]); case 2: return -__sin(y[0], y[1], 1); default: return -__cos(y[0], y[1]); } }
Но каким бы ни был умным компилятор, он не даст нам здесь больше 10-15% прироста, какие режимы не включай. Копаем глубже, включаем -03 + -math:precise и смотрим что нагенерило в этот раз. На этот раз, получили алиас на другую функцию, которая называется __sin_dorf_nx_502, код очень линейный, мало ифов, сложения и умножения
__sin_dorf_nx_502
__sin_dorf_nx_502(float): push {r4, r5, r6, r7, r8, r9, r10, r11, lr} add r11, sp, #28 sub sp, sp, #4 ldr r1, .EVRT0_0 mov r4, r0 bl __aeabi_fmul and r1, r4, #-2147483648 orr r1, r1, #1056964608 bl __aeabi_fadd bl __aeabi_f2uiz mov r9, r0 bl __aeabi_ui2f ldr r8, .EVRT0_1 mov r5, r0 mov r1, r5 ldr r0, [r8, #24] bl __aeabi_fmul mov r1, r0 mov r0, r4 bl __aeabi_fsub mov r4, r0 ldr r0, [r8, #28] mov r1, r5 bl __aeabi_fmul mov r1, r0 mov r0, r4 bl __aeabi_fsub mov r1, r0 mov r7, r0 bl __aeabi_fmul mov r5, r0 mov r0, r7 mov r1, r5 bl __aeabi_fmul mov r6, r0 ldr r0, [r8, #8] ldm r8, {r4, r10} mov r1, r5 str r0, [sp] ldr r0, [r8, #12] bl __aeabi_fmul ldr r1, [r8, #16] bl __aeabi_fadd mov r1, r0 mov r0, r5 bl __aeabi_fmul ldr r1, [r8, #20] bl __aeabi_fadd mov r1, r0 mov r0, r6 bl __aeabi_fmul mov r1, r0 mov r0, r7 bl __aeabi_fadd mov r7, r0 mov r0, r4 mov r1, r5 bl __aeabi_fmul mov r1, r0 mov r0, r10 bl __aeabi_fadd mov r1, r0 mov r0, r5 bl __aeabi_fmul mov r1, r0 ldr r0, [sp] bl __aeabi_fadd mov r1, r0 mov r0, r5 bl __aeabi_fmul mov r1, #1065353216 bl __aeabi_fadd tst r9, #1 moveq r0, r7 tst r9, #2 eorne r0, r0, #-2147483648 sub sp, r11, #28 pop {r4, r5, r6, r7, r8, r9, r10, r11, lr} bx lr .EVRT0_0: .long 1059256694 @ 0x3f22f976 .EVRT0_1: .long .KI_MergedGlobals ... .KI_MergedGlobals: // вот тут блок параметров для полинома, он располагался гдето в секции данных .long 3132246419 @ float -0.001360224 .long 1026203702 @ float 0.04165669 .long 3204448223 @ float -0.4999990 .long 3108801646 @ float -1.950727E-4 .long 1007190850 @ float 0.008332075 .long 3190467233 @ float -0.1666665 .long 1070141402 @ float 1.570796 /// вот это Pi/2 .long 866263401 @ float 7.549790E-8
Очень похоже на аппроксимацию синуса каким-то полиномом, если попробовать восстановить в псевдо коде примерную логику работы, то будет похоже на код ниже. И это уже похоже на то, что показывает бенчмарк.
float __sin_dorf_nx_502(float x) { const float EVRT0_0 = 0.636618972f; const float EVRT0_1 = 1.0f; const float KI_MergedGlobals[8] = { -0.001360224f, 0.04165669f, -0.4999990f, -1.950727E-4f, 0.008332075f, -0.1666665f, 1.570796f, // PI / 2 7.549790E-8f }; float result; float temp = EVRT0_0 * x; temp = temp + 1056964608.0f; int intTemp = (int)temp; float floatTemp = (float)intTemp; float r5 = floatTemp; floatTemp = r5 * KI_MergedGlobals[0]; float r4 = x - floatTemp; floatTemp = r5 * KI_MergedGlobals[1]; r4 = r4 - floatTemp; float r7 = r4 * r4; floatTemp = r7 * KI_MergedGlobals[2]; floatTemp = floatTemp + KI_MergedGlobals[3]; floatTemp = floatTemp * r7; floatTemp = floatTemp + KI_MergedGlobals[4]; floatTemp = floatTemp * r7; floatTemp = floatTemp + KI_MergedGlobals[5]; floatTemp = floatTemp * r7; floatTemp = floatTemp + EVRT0_1; if (intTemp & 1) { floatTemp = floatTemp * r4; } else { floatTemp = -floatTemp * r4; } result = floatTemp; return result; }
Думаю вы уже представляете, что будет когда соберем сэмпл с -03 + -math:fast? Правильно, алиас еще на одну функцию. Финальный вариант, который используется в релизной сборке, с которым игры отправлялись на сертификацию
__sin_corf_nx_502
__sin_corf_nx_502(float): push {r4, r5, r6, r7, r8, r9, r10, r11, lr} add r11, sp, #28 sub sp, sp, #4 bl __aeabi_f2d ldr r2, .DUIE0_0 ldr r3, .DUIE0_1 mov r5, r0 mov r6, r1 bl __aeabi_dmul bl __aeabi_d2iz mov r10, r0 bl __aeabi_i2d ldr r3, .DUIE0_2 mov r2, #1610612736 bl __aeabi_dmul mov r2, r0 mov r3, r1 mov r0, r5 mov r1, r6 bl __aeabi_dadd bl __aeabi_d2f ldr r1, .DUIE0_3 mov r7, r0 and r0, r10, #255 ldr r8, [r1] ldr r0, [r8, r0, lsl #2] bl __aeabi_f2d mov r3, #266338304 mov r2, #0 str r0, [sp] mov r9, r1 orr r3, r3, #-1342177280 bl __aeabi_dmul mov r5, r0 mov r0, r7 mov r6, r1 bl __aeabi_f2d mov r7, r0 mov r4, r1 mov r0, r5 mov r1, r6 mov r2, r7 mov r3, r4 bl __aeabi_dmul mov r5, r0 add r0, r10, #64 mov r6, r1 and r0, r0, #255 ldr r0, [r8, r0, lsl #2] bl __aeabi_f2d mov r2, r5 mov r3, r6 bl __aeabi_dadd mov r2, r7 mov r3, r4 bl __aeabi_dmul ldr r2, [sp] mov r3, r9 bl __aeabi_dadd bl __aeabi_d2f sub sp, r11, #28 pop {r4, r5, r6, r7, r8, r9, r10, r11, lr} bx lr .DUIE0_0: .long 1682373044 @ 0x6446f9b4 .DUIE0_1: .long 1078222640 @ 0x40445f30 .DUIE0_2: .long 3214483963 @ 0xbf9921fb .DUIE0_3: .long __sin_corf_duie_nx_502
Псевдокодом для этой ассемблерной части будет такой, и судя по логике работы это вычисление синуса с помощью таблицы. Быстрее всех, точность не сильно страдает.
float __sin_corf_nx_502(float x) { const double DUIE0_0 = 1.4636916370976118; const double DUIE0_1 = 3.141592653589793; // PI const double DUIE0_2 = -0.0009424784718423055; const float* DUIE0_3 = __sin_corf_duie_nx_502; double temp = x * DUIE0_0; temp = temp + DUIE0_1; int intTemp = static_cast<int>(temp); float floatTemp = static_cast<float>(intTemp); float r5 = floatTemp; floatTemp = r5 * DUIE0_2; float r4 = x - floatTemp; double r7 = r4 * r4; floatTemp = static_cast<float>(r7); int index = intTemp & 255; float floatTemp2 = *(DUIE0_3 + index); double r2 = floatTemp2; double r3 = static_cast<double>(intTemp >> 8); double r0 = static_cast<double>(r5); double r1 = static_cast<double>(r4); r0 = r0 + r1; r0 = r0 * r2; r0 = r0 + r3; float result = static_cast<float>(r0); return result; }
Если еще покопаться в интернетах на тему быстрого синуса, то можно найти например еще такую аппроксимацию https://en.wikipedia.org/wiki/Bhaskara_I's_sine_approximation_formula ошибка вычислений по этой формуле находится в пределах 0.001, что зачастую хватает для большинства задач.
Или обратить внимание на статью Ларса Ниланда (Lars Nyland, nVidia), одного из авторов CUDA.
Для синусов в движке Dagor используется своя реализация, она одинаково хорошо и быстро ведет себя на всех платформах, найти его можно тут. Не сочтите за рекламу, я давно уже не работаю в компании, может кому пригодится. Самых добрых приветов моим бывшим коллегам.
Выводы
Никакой практической пользы мы с этого не получили, но было интересно разобраться в особенностях работы сдк. Честь и хвала японским инженерам, подарившим миру эту замечательную консоль. У меня только один вопрос, нафига три?
UPD: в личку написал товарищ, что так Нинтендо отслеживала использование старых и паленых сдк, без активного ключа разработчика всегда вызывалась первая реализация, таким образом билд, который пытался пройти сертификацию помогал находить украденные аккаунты и ключи. Проверить эти данные я не могу, остается только поверить :)
