Зачем в Switch SDK три разных sin?

Работая в компании 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(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(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(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: в личку написал товарищ, что так Нинтендо отслеживала использование старых и паленых сдк, без активного ключа разработчика всегда вызывалась первая реализация, таким образом билд, который пытался пройти сертификацию помогал находить украденные аккаунты и ключи. Проверить эти данные я не могу, остается только поверить 🙂