Оптимизация sincos функции

Давненько я ничего не писал в свой блог и на это было несколько причин. Во-первых, я полностью погрузился в написание своего нового проекта - sdr-modem. Не сказать, что он простой, но он гармонично развивает идеи sdr-server и не содержит того, о чём хотелось бы написать. Во-вторых, ничего особенного не происходит.

Тем не менее хочется написать небольшую заметку о том, как я измерял производительность функции sincos, и что из этого получилось.

sincos

Вычисление синуса и косинуса - это достаточно дорогие операции. Именно поэтому их очень часто пытаются оптимизировать. Однако, иногда нужно посчитать одновременно и синус и косинус. Например, при частотной модуляции сигнала для каждого значения сигнала нужно увеличивать фазу и считать для неё синус и косинус. Невероятное совпадение, но это как раз то, над чем я сейчас работаю. В таких случаях есть очень странное ощущение, что вычислив синус угла, можно чуть быстрее посчитать косинус. И, действительно, в библиотеке C есть отдельная функция для этого - sincos. Она позволяет вычислить одновременно синус и косинус угла.

Однако, к тому моменту как я нашёл эту функцию, я знал ещё как минимум о двух других реализациях - volk и приближённое вычисление в gnuradio. При таком разнообразии подходов глаза начинают разбегаться и нужно срочно измерять! Что я и сделал, написав небольшое приложение для каждой из реализаций.

Наивная реализация

#include <math.h>

void calc_naive(float *input, size_t len, float complex *output, size_t output_len) {
    for (size_t i = 0; i < len; i++) {
        output[i] = cosf(input[i]) + I * sinf(input[i]);
    }
}

Ничего особенного, просто вычисление синуса и косинуса каждый раз. Здесь и далее все массивы заранее созданы и проинициализированы.

Стандартная функция sincos

#include <math.h>

void calc_sincos(float *input, size_t len, float complex *output, size_t output_len) {
    float real;
    float imag;
    for (size_t i = 0; i < len; i++) {
        __sincosf(input[i], &real, &imag);
        output[i] = real + I * imag;
    }
}

__sincosf берёт углы один за другим и последовательно рассчитывает значение. Результат кладётся в выходной массив комплексных чисел.

Когда я писал этот метод, то случайно навёл курсор на функцию почитать документацию и увидел интересное заявление:

Оказывается в MacOS есть специальная библиотека для работы с SIMD инструкциями и она поддерживает sincos!

sincos из библиотеки Accelerate

Accelerate - это Framework, если использовать терминологию Apple. Подключается он просто:

add_compile_options("-F/Library/Frameworks/")
link_libraries("-F/Library/Frameworks/")

find_library(LIB_ACCELERATE accelerate)
target_link_libraries(sdr_modemLib ${LIB_ACCELERATE})

После этого его можно использовать:

#include <Accelerate/Accelerate.h>

void calc_accelerate(float *sin, float *cos, float *input, size_t len, float complex *output, size_t output_len) {
    vvsincosf(sin, cos, input, (const int *) &len);
    for (size_t i = 0; i < len; i++) {
        output[i] = cos[i] + I * sin[i];
    }
}

Эта функция возвращает отдельно массив синусов и массив косинусов, так что всё равно придётся выполнять de-interleaving, чтобы получить массив комплексных чисел.

volk

В volk тоже есть функции для работы с синусом и косинусом. Однако, это две совершенно разные функции, которые работают независимо. Но это не беда:

#include <volk/volk.h>

void calc_volk(float *sin, float *cos, float *input, size_t len, float complex *output, size_t output_len) {
    volk_32f_sin_32f(sin, input, (unsigned int) len);
    volk_32f_cos_32f(cos, input, (unsigned int) len);
    for (size_t i = 0; i < len; i++) {
        output[i] = cos[i] + I * sin[i];
    }
}

Вызываем сначала расчёт синусов, затем косинусов, а в конце de-interleaving.

Приближённое вычисление и таблицы поиска

Этот метод активно используется в gnuradio. Насколько я понял, они сначала конвертируют угол в целочисленное значение, а потом хитрым способом ищут в таблице поиска. Код из gnuradio не так-то просто добавить в проект, поэтому я скопировал его как-есть.

void cals_gnuradio(float *input, size_t len, float complex *output, size_t output_len) {
    float oi, oq;
    for (size_t i = 0; i < len; i++) {
        float d_phase = fmod(input[i] + PI, 2.0f * PI) - PI;
        int32_t angle = float_to_fixed(d_phase);
        gnuradio_sincos(angle, &oq, &oi);
        output[i] = oi + I * oq;
    }
}

Результаты

Результат выполнения программы для каждого из метода:

Второй столбец показывает время выполнения в попугаях на случайно большом количестве входящих данных. Третий и четвёртый столбцы - это значение комплексного числа из результата по случайному индексу. Я добавил его просто, чтобы убедиться в одинаковых результатах.

Как видно из времени выполнения, самый быстрый способ - это volk. Он почти в 3.2 раза быстрее стандартной библиотеки. И даже быстрее Accelerate.

Ещё один удивительный факт - это крайне медленный способ приближённого расчёта. Я бы ожидал, что он на порядок быстрее всех остальных, но на практике он в 7 раз медленнее самого быстрого способа! Я выполнил программу несколько раз и каждый раз получил примерно одинаковые результаты. Так что, это вряд ли связано с загрузкой системы.

Функция sincos значительно проигрывает наивной реализации. Я не знаю почему.

Несмотря на то, что Accelerate не самый быстрый способ расчёта sincos (примерно в 2.8 раз быстрее стандартной функции), я был приятно удивлён. Буквально за пару строчек кода я получил стандартную библиотеку, которая неплохо оптимизирована для MacOS. А если пойти на официальный сайт, то можно увидеть достаточно интересный список поддерживаемых функций. Помимо DSP, эта библиотека поддерживает работу с нейронными сетями, векторами, матрицами, линейной алгеброй. Зная, насколько Apple заморачивается с поддержкой нейронных сетей в своих процессорах, я представляю насколько круто иметь железо и библиотеку работающих в тандеме. Не удивительно, что софт в MacOS по ощущениям просто летает.