Эффективно оценивайте большие полиномы с помощью SIMD

У меня есть довольно большие (20-40 градусов) медленно сходящиеся (иногда) многочлены с плавающей запятой. Я хотел бы оптимизировать их оценку с помощью SIMD (SSE2, AVX1, AVX-512). Мне нужны решения как с плавающей запятой-32, так и с двойным-64.

Значения коэффициентов являются константами, заданными заранее, а значение X для оценки поли задается в качестве аргумента функции.

Важное примечание. У меня есть только один вход X для моей функции. Поэтому я не могу выполнять вертикальную оптимизацию, вычисляя полигон для 8-16 Xs одновременно. Это означает, что мне нужна горизонтальная оптимизация в рамках оценки для одного X.

Я создал связанный вопрос, который помогает мне вычислить степени X (например, X^1, X^2, ..., X^8), необходимые для оценки SIMD.

Очевидно, что SIMD следует использовать только после некоторого порога полиномиальной степени, для довольно малых полигонов можно использовать метод Хорнера (или Эстрина), основанный на методе Хорнера (или Эстрина) как здесь. Также ширина SIMD (128, 256 или 512) должна выбираться в зависимости от степени полиномии.

Ниже я реализовал вариант AVX-256-Float32, используя своего рода модифицированный метод Хорнера, подходящий для SIMD (умножение на x^8 вместо x^1). Спасибо @PeterCordes за быструю горизонтальную сумму tutorial. Нажмите на ссылку «попробовать онлайн», она содержит более крупный код, который также имеет справочную простую оценку для сравнения и измерения времени:

Попробуйте онлайн!

template <size_t S, size_t I, typename MT = __m256, size_t Cnt>
inline MT EvalPoly8xF32Helper(MT xhi,
        std::array<float, Cnt> const & A, MT r = _mm256_undefined_ps()) {
    size_t constexpr K = 8;
    if constexpr(I + K >= S)
        r = _mm256_load_ps(&A[I]);
    else {
        #ifdef __FMA__
            r = _mm256_fmadd_ps(r, xhi, _mm256_load_ps(&A[I]));
        #else
            r = _mm256_add_ps(_mm256_mul_ps(r, xhi),
                _mm256_load_ps(&A[I]));
        #endif
    }
    if constexpr(I < K)
        return r;
    else
        return EvalPoly8xF32Helper<S, I - K>(xhi, A, r);
}

inline float _mm_fast_hsum_ps(__m128 v) {
    __m128 shuf = _mm_movehdup_ps(v);
    __m128 sums = _mm_add_ps(v, shuf);
    shuf        = _mm_movehl_ps(shuf, sums);
    sums        = _mm_add_ss(sums, shuf);
    return        _mm_cvtss_f32(sums);
}

template <size_t S, size_t Cnt>
inline float EvalPoly8xF32(
        float x, std::array<float, Cnt> const & A) {
    auto constexpr K = 8;
    auto const x2 = x * x, x4 = x2 * x2, x8 = x4 * x4, x3 = x2 * x;
    auto const powsx = _mm256_setr_ps(
        1, x, x2, x3, x4, x4 * x, x4 * x2, x4 * x3);
    auto r0 = EvalPoly8xF32Helper<S, (S - 1) / K * K>(
        _mm256_set1_ps(x8), A);
    r0 = _mm256_mul_ps(r0, powsx);
    return _mm_fast_hsum_ps(_mm_add_ps(
        _mm256_castps256_ps128(r0), _mm256_extractf128_ps(r0, 1)));
}

Как видно, SIMD-версия дает достаточно большое ускорение по сравнению с эталонной простой реализацией. Для AVX1-256-float32 и степени 32 это дает ускорение примерно в 4.5x раз (для степени 16 дает 1.8x ускорение, что тоже хорошо)! Очевидно, что даже простое использование FMA инструкций внутри эталонных реализаций уже заметно улучшит скорость вычислений.

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

Чаще всего будет использоваться целевой ЦП: Intel Xeon Gold 6230 с AVX-512, поэтому мне нужно оптимизировать код для него.