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