Самый быстрый и точный способ преобразования вектора целых чисел в числа с плавающей точкой от 0 до 1 - PullRequest
Новая
       41

Самый быстрый и точный способ преобразования вектора целых чисел в числа с плавающей точкой от 0 до 1

0 голосов
/ 25 февраля 2019

Рассмотрим случайно сгенерированный __m256i вектор.Есть ли более точный способ преобразовать их в __m256 вектор чисел с плавающей запятой между 0 (включительно) и 1 (исключительно), чем деление на float(1ull<<32)?

Вот что я пробовал до сих поргде iRand - это ввод, а ans - это вывод: 

const __m256 fRand = _mm256_cvtepi32_ps(iRand);
const __m256 normalized = _mm256_div_ps(fRand, _mm256_set1_ps(float(1ull<<32)));
const __m256 ans = _mm256_add_ps(normalized, _mm256_set1_ps(0.5f));

Ответы [ 2 ]

0 голосов
/ 26 февраля 2019

Во-первых, без деления, замените его умножением.Хотя @Soonts может быть достаточно хорошим для вас, я могу заметить только из-за использования отображения на интервал [1 ... 2), он выдает равномерные двоичные рациональные числа в форме k / 2 −23 ,что составляет половину того, что может быть сгенерировано.Я предпочитаю метод из S.Vigna (внизу), причем все двоичные рациональные числа вида k / 2 −24 одинаково вероятны.

Code, VC++ 2019, x64, Win10, Intel i7 Skylake 

#include <random>

#include "immintrin.h"

auto p256_dec_u32(__m256i in) -> void {
    alignas(alignof(__m256i)) uint32_t v[8];
    _mm256_store_si256((__m256i*)v, in);
    printf("v8_u32: %u %u %u %u %u %u %u %u\n", v[0], v[1], v[2], v[3], v[4], v[5], v[6], v[7]);
}

auto p256_dec_f32(__m256 in) -> void {
    alignas(alignof(__m256)) float v[8];
    _mm256_store_ps(v, in);
    printf("v8_float: %e %e %e %e %e %e %e %e\n", v[0], v[1], v[2], v[3], v[4], v[5], v[6], v[7]);
}

auto main() -> int {
    const float c = 0x1.0p-24f; // or (1.0f / (uint32_t(1) << 24));

    const int N = 1000000;

    std::mt19937 rng{ 987654321ULL };

    __m256 sum = _mm256_set1_ps(0.0f);

    for (int k = 0; k != N; ++k) {
        alignas(alignof(__m256i)) uint32_t rnd[8] = { rng(), rng(), rng(), rng(), rng(), rng(), rng(), rng() };

        __m256i r = _mm256_load_si256((__m256i*)rnd);
        __m256  q = _mm256_mul_ps(_mm256_cvtepi32_ps(_mm256_srli_epi32(r, 8)), _mm256_set1_ps(c));

        sum = _mm256_add_ps(sum, q);
    }

    sum = _mm256_div_ps(sum, _mm256_set1_ps((float)N)); // computing average

    p256_dec_f32(sum);

    return 0;
}

с выводом 

5.002970e-01 4.997833e-01 4.996118e-01 5.004955e-01 5.002163e-01 4.997193e-01 4.996586e-01 5.001499e-01
0 голосов
/ 25 февраля 2019

Версия ниже должна быть быстрее, по сравнению с вашей первоначальной версией, которая использует _mm256_div_ps

vdivps довольно медленно, например, на моем Haswell Xeon это задержка 18-21 циклов, пропускная способность 14 циклов.Более новые процессоры работают лучше, кстати, 11/5 на Skylake, 10/6 на Ryzen.

Как сказано в комментариях, производительность можно исправить, заменив деление на умножение, а затем улучшив FMA.Проблема с подходом - качество распространения.Если вы попытаетесь получить эти числа в выходном интервале с помощью режима округления или отсечения, вы введете пики в распределении вероятностей выходных чисел.

Моя реализация тоже не идеальна, она не выводитвсе возможные значения в выходном интервале пропускают многие представимые числа с плавающей точкой, особенно около 0. Но, по крайней мере, распределение очень равномерное. 

__m256 __vectorcall randomFloats( __m256i randomBits )
{
    // Convert to random float bits
    __m256 result = _mm256_castsi256_ps( randomBits );

    // Zero out exponent bits, leave random bits in mantissa.
    // BTW since the mask value is constexpr, we don't actually need AVX2 instructions for this, it's just easier to code with set1_epi32.
    const __m256 mantissaMask = _mm256_castsi256_ps( _mm256_set1_epi32( 0x007FFFFF ) );
    result = _mm256_and_ps( result, mantissaMask );

    // Set sign + exponent bits to that of 1.0, which is sign=0, exponent=2^0.
    const __m256 one = _mm256_set1_ps( 1.0f );
    result = _mm256_or_ps( result, one );

    // Subtract 1.0. The above algorithm generates floats in range [1..2).
    // Can't use bit tricks to generate floats in [0..1) because it would cause them to be distributed very unevenly.
    return _mm256_sub_ps( result, one );
}

Обновление: Если вы хотите повысить точность, используйтеследующая версия.Но он больше не «самый быстрый». 

__m256 __vectorcall randomFloats_32( __m256i randomBits )
{
    // Convert to random float bits
    __m256 result = _mm256_castsi256_ps( randomBits );
    // Zero out exponent bits, leave random bits in mantissa.
    const __m256 mantissaMask = _mm256_castsi256_ps( _mm256_set1_epi32( 0x007FFFFF ) );
    result = _mm256_and_ps( result, mantissaMask );
    // Set sign + exponent bits to that of 1.0, which is sign=0, exponent = 2^0.
    const __m256 one = _mm256_set1_ps( 1.0f );
    result = _mm256_or_ps( result, one );
    // Subtract 1.0. The above algorithm generates floats in range [1..2).
    result = _mm256_sub_ps( result, one );

    // Use 9 unused random bits to add extra randomness to the lower bits of the values.
    // This increases precision to 2^-32, however most floats in the range can't store that many bits, fmadd will only add them for small enough values.

    // If you want uniformly distributed floats with 2^-24 precision, replace the second argument in the following line with _mm256_set1_epi32( 0x80000000 ).
    // In this case you don't need to set rounding mode bits in MXCSR.
    __m256i extraBits = _mm256_and_si256( randomBits, _mm256_castps_si256( mantissaMask ) );
    extraBits = _mm256_srli_epi32( extraBits, 9 );
    __m256 extra = _mm256_castsi256_ps( extraBits );
    extra = _mm256_or_ps( extra, one );
    extra = _mm256_sub_ps( extra, one );
    _MM_SET_ROUNDING_MODE( _MM_ROUND_DOWN );
    constexpr float mul = 0x1p-23f; // The initial part of the algorithm has generated uniform distribution with the step 2^-23.
    return _mm256_fmadd_ps( extra, _mm256_set1_ps( mul ), result );
}
...