Почему в Switch SDK имеются три различные функции sin?
To sine or not to sine?
Работая в компании 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(temp);
float floatTemp = static_cast(intTemp);
float r5 = floatTemp;
floatTemp = r5 * DUIE0_2;
float r4 = x - floatTemp;
double r7 = r4 * r4;
floatTemp = static_cast(r7);
int index = intTemp & 255;
float floatTemp2 = *(DUIE0_3 + index);
double r2 = floatTemp2;
double r3 = static_cast(intTemp >> 8);
double r0 = static_cast(r5);
double r1 = static_cast(r4);
r0 = r0 + r1;
r0 = r0 * r2;
r0 = r0 + r3;
float result = static_cast(r0);
return result;
} Если еще покопаться в интернетах на тему быстрого синуса, то можно найти например еще такую аппроксимацию https://en.wikipedia.org/wiki/Bhaskara_I's_sine_approximation_formula ошибка вычислений по этой формуле находится в пределах 0.001, что зачастую хватает для большинства задач.
Или обратить внимание на статью Ларса Ниланда (Lars Nyland, nVidia), одного из авторов CUDA https://www.researchgate.net/publication/2930018_Fast_Trigonometric_Functions_Using_Intel's_Sse2
Для синусов в движке Dagor используется своя реализация, она одинаково хорошо и быстро ведет себя на всех платформах, найти его можно тут (https://github.com/GaijinEntertainment/vecmath/blob/master/dag_vecMath_trig.h). Не сочтите за рекламу, я давно уже не работаю в компании, может кому пригодится. Самых добрых приветов моим бывшим коллегам.
Выводы:
Никакой практической пользы мы с этого не получили, но было интересно разобраться в особенностях работы сдк. Честь и хвала японским инженерам, подарившим миру эту замечательную консоль. У меня только один вопрос, нафига три?
UPD: в личку написал товарищ, что так Нинтендо отслеживала использование старых и паленых сдк, без активного ключа разработчика всегда вызывалась первая реализация, таким образом билд, который пытался пройти сертификацию помогал находить украденные аккаунты и ключи. Проверить эти данные я не могу, остается только поверить 🙂
Похожие статьи
Хидео Кодзима рассказал об ограничениях Death Stranding из-за небольшого бюджета и команды
Компания Inart приобрела права на создание коллекционных фигурок по франшизе League of Legends
Инсайдер назвал Marvel’s Wolverine второй по стоимости игрой Insomniac Games
Blizzard дала фанатам новый повод ждать выхода World of Warcraft на консолях
Даниэль Вавра встал на защиту DLSS из-за того, что в Kingdom Come: Deliverance II технология позволила реализовать его видение
Низкая оценка Marvel’s Wolverine может лишить Sony главной номинации на The Game Awards
Экшен X-Men Origins: Wolverine удалось разработать без переработок благодаря переносу фильма
Разработчик Total War Warhammer 40000 рассказал о влиянии настольной игры на стратегию