Блог

Мне нужна кроссплатформенная библиотека / алгоритм, который будет преобразовывать между 32-разрядными и 16-разрядными числами с плавающей запятой. Мне не нужно выполнять математические вычисления с 16-разрядными числами; Мне просто нужно уменьшить размер 32-разрядных чисел с плавающей запятой, чтобы их можно было отправлять по сети. Я работаю на C .

Я понимаю, какую точность я бы потерял, но это нормально для моего приложения.

16-разрядный формат IEEE был бы отличным.

Полное преобразование с одинарной точностью до половинной точности. Это прямая копия моей версии SSE, поэтому она не содержит ветвей. Он использует тот факт, что -true == ~0 для предварительной выборки без ветвей (GCC преобразует if операторы в нечестивый беспорядок условных переходов, в то время как Clang просто преобразует их в условные перемещения.)

Обновление (2019-11-04): переработано для поддержки значений одинарной и двойной точности с полностью правильным округлением. Я также помещаю соответствующий if оператор над каждым выбором без ветвей в качестве комментария для ясности. Все входящие NAN преобразуются в базовый тихий NaN для скорости и работоспособности, поскольку нет способа надежно преобразовать встроенное сообщение NaN между форматами.

 #include <cstdint> // uint32_t, uint64_t, etc.
#include <cstring> // memcpy
#include <climits> // CHAR_BIT
#include <limits>  // numeric_limits
#include <utility> // is_integral_v, is_floating_point_v, forward

namespace std
{
  template< typename T , typename U >
  T bit_cast( Uamp;amp; u ) {
    static_assert( sizeof( T ) == sizeof( U ) );
    union { T t; }; // prevent construction
    std::memcpy( amp;t, amp;u, sizeof( t ) );
    return t;
  }
} // namespace std

template< typename T > struct native_float_bits;
template<> struct native_float_bits< float >{ using type = std::uint32_t; };
template<> struct native_float_bits< double >{ using type = std::uint64_t; };
template< typename T > using native_float_bits_t = typename native_float_bits< T >::type;

static_assert( sizeof( float ) == sizeof( native_float_bits_t< float > ) );
static_assert( sizeof( double ) == sizeof( native_float_bits_t< double > ) );

template< typename T, int SIG_BITS, int EXP_BITS >
struct raw_float_type_info {
  using raw_type = T;

  static constexpr int sig_bits = SIG_BITS;
  static constexpr int exp_bits = EXP_BITS;
  static constexpr int bits = sig_bits   exp_bits   1;

  static_assert( std::is_integral_v< raw_type > );
  static_assert( sig_bits >= 0 );
  static_assert( exp_bits >= 0 );
  static_assert( bits <= sizeof( raw_type ) * CHAR_BIT );

  static constexpr int exp_max = ( 1 << exp_bits ) - 1;
  static constexpr int exp_bias = exp_max >> 1;

  static constexpr raw_type sign = raw_type( 1 ) << ( bits - 1 );
  static constexpr raw_type inf = raw_type( exp_max ) << sig_bits;
  static constexpr raw_type qnan = inf | ( inf >> 1 );

  static constexpr auto abs( raw_type v ) { return raw_type( v amp; ( sign - 1 ) ); }
  static constexpr bool is_nan( raw_type v ) { return abs( v ) > inf; }
  static constexpr bool is_inf( raw_type v ) { return abs( v ) == inf; }
  static constexpr bool is_zero( raw_type v ) { return abs( v ) == 0; }
};
using raw_flt16_type_info = raw_float_type_info< std::uint16_t, 10, 5 >;
using raw_flt32_type_info = raw_float_type_info< std::uint32_t, 23, 8 >;
using raw_flt64_type_info = raw_float_type_info< std::uint64_t, 52, 11 >;
//using raw_flt128_type_info = raw_float_type_info< uint128_t, 112, 15 >;

template< typename T, int SIG_BITS = std::numeric_limits< T >::digits - 1,
  int EXP_BITS = sizeof( T ) * CHAR_BIT - SIG_BITS - 1 >
struct float_type_info 
: raw_float_type_info< native_float_bits_t< T >, SIG_BITS, EXP_BITS > {
  using flt_type = T;
  static_assert( std::is_floating_point_v< flt_type > );
};

template< typename E >
struct raw_float_encoder
{
  using enc = E;
  using enc_type = typename enc::raw_type;

  template< bool DO_ROUNDING, typename F >
  static auto encode( F value )
  {
    using flt = float_type_info< F >;
    using raw_type = typename flt::raw_type;
    static constexpr auto sig_diff = flt::sig_bits - enc::sig_bits;
    static constexpr auto bit_diff = flt::bits - enc::bits;
    static constexpr auto do_rounding = DO_ROUNDING amp;amp; sig_diff > 0;
    static constexpr auto bias_mul = raw_type( enc::exp_bias ) << flt::sig_bits;
    if constexpr( !do_rounding ) { // fix exp bias
      // when not rounding, fix exp first to avoid mixing float and binary ops
      value *= std::bit_cast< F >( bias_mul );
    }
    auto bits = std::bit_cast< raw_type >( value );
    auto sign = bits amp; flt::sign; // save sign
    bits ^= sign; // clear sign
    auto is_nan = flt::inf < bits; // compare before rounding!!
    if constexpr( do_rounding ) {
      static constexpr auto min_norm = raw_type( flt::exp_bias - enc::exp_bias   1 ) << flt::sig_bits;
      static constexpr auto sub_rnd = enc::exp_bias < sig_diff
        ? raw_type( 1 ) << ( flt::sig_bits - 1   enc::exp_bias - sig_diff )
        : raw_type( enc::exp_bias - sig_diff ) << flt::sig_bits;
      static constexpr auto sub_mul = raw_type( flt::exp_bias   sig_diff ) << flt::sig_bits;
      bool is_sub = bits < min_norm;
      auto norm = std::bit_cast< F >( bits );
      auto subn = norm;
      subn *= std::bit_cast< F >( sub_rnd ); // round subnormals
      subn *= std::bit_cast< F >( sub_mul ); // correct subnormal exp
      norm *= std::bit_cast< F >( bias_mul ); // fix exp bias
      bits = std::bit_cast< raw_type >( norm );
      bits  = ( bits >> sig_diff ) amp; 1; // add tie breaking bias
      bits  = ( raw_type( 1 ) << ( sig_diff - 1 ) ) - 1; // round up to half
      //if( is_sub ) bits = std::bit_cast< raw_type >( subn );
      bits ^= -is_sub amp; ( std::bit_cast< raw_type >( subn ) ^ bits );
    }
    bits >>= sig_diff; // truncate
    //if( enc::inf < bits ) bits = enc::inf; // fix overflow
    bits ^= -( enc::inf < bits ) amp; ( enc::inf ^ bits );
    //if( is_nan ) bits = enc::qnan;
    bits ^= -is_nan amp; ( enc::qnan ^ bits );
    bits |= sign >> bit_diff; // restore sign
    return enc_type( bits );
  }

  template< typename F >
  static F decode( enc_type value )
  {
    using flt = float_type_info< F >;
    using raw_type = typename flt::raw_type;
    static constexpr auto sig_diff = flt::sig_bits - enc::sig_bits;
    static constexpr auto bit_diff = flt::bits - enc::bits;
    static constexpr auto bias_mul = raw_type( 2 * flt::exp_bias - enc::exp_bias ) << flt::sig_bits;
    raw_type bits = value;
    auto sign = bits amp; enc::sign; // save sign
    bits ^= sign; // clear sign
    auto is_norm = bits < enc::inf;
    bits = ( sign << bit_diff ) | ( bits << sig_diff );
    auto val = std::bit_cast< F >( bits ) * std::bit_cast< F >( bias_mul );
    bits = std::bit_cast< raw_type >( val );
    //if( !is_norm ) bits |= flt::inf;
    bits |= -!is_norm amp; flt::inf;
    return std::bit_cast< F >( bits );
  }
};

using flt16_encoder = raw_float_encoder< raw_flt16_type_info >;

template< typename F >
auto quick_encode_flt16( F amp;amp; value )
{ return flt16_encoder::encode< false >( std::forward< F >( value ) ); }

template< typename F >
auto encode_flt16( F amp;amp; value )
{ return flt16_encoder::encode< true >( std::forward< F >( value ) ); }

template< typename F = float, typename X >
auto decode_flt16( X amp;amp; value )
{ return flt16_encoder::decode< F >( std::forward< X >( value ) ); }
  

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

Половина с плавающей запятой:
float f = ((hamp;0x8000)<<16) | (((hamp;0x7c00) 0x1C000)<<13) | ((hamp;0x03FF)<<13);

С плавающей запятой до половины:
uint32_t x = *((uint32_t*)amp;f);
uint16_t h = ((x>>16)amp;0x8000)|((((xamp;0x7f800000)-0x38000000)>>13)amp;0x7c00)|((x>>13)amp;0x03ff);

std::frexp извлекает значение и показатель степени из обычных чисел с плавающей запятой или удвоений — затем вам нужно решить, что делать с показателями, которые слишком велики, чтобы поместиться в число с плавающей запятой половинной точности (насыщать …?), Соответствующим образом скорректировать и сложить число с половинной точностью. В этой статье приведен исходный код C, чтобы показать вам, как выполнить преобразование.

Учитывая ваши потребности (-1000, 1000), возможно, было бы лучше использовать представление с фиксированной запятой.

 //change to 20000 to SHORT_MAX if you don't mind whole numbers
//being turned into fractional ones
const int compact_range = 20000;

short compactFloat(double input) {
    return round(input * compact_range / 1000);
}
double expandToFloat(short input) {
    return ((double)input) * 1000 / compact_range;
}
  

Это даст вам точность с точностью до 0,05. Если вы измените 20000 на SHORT_MAX, вы получите немного больше точности, но некоторые целые числа на другом конце будут десятичными.

Почему так сложно? Моя реализация не нуждается в какой-либо дополнительной библиотеке, соответствует формату IEEE-754 FP16, управляет как нормализованными, так и денормализованными числами, не имеет ветвей, занимает около 40 тактов для обратного и обратного преобразования и канав NaN или Inf для расширенного диапазона. В этом волшебная сила битовых операций.

 typedef unsigned short ushort;
typedef unsigned int uint;

uint as_uint(const float x) {
    return *(uint*)amp;x;
}
float as_float(const uint x) {
    return *(float*)amp;x;
}

float half_to_float(const ushort x) { // IEEE-754 16-bit floating-point format (without infinity): 1-5-10, exp-15,  -131008.0,  -6.1035156E-5,  -5.9604645E-8, 3.311 digits
    const uint e = (xamp;0x7C00)>>10; // exponent
    const uint m = (xamp;0x03FF)<<13; // mantissa
    const uint v = as_uint((float)m)>>23; // evil log2 bit hack to count leading zeros in denormalized format
    return as_float((xamp;0x8000)<<16 | (e!=0)*((e 112)<<23|m) | ((e==0)amp;(m!=0))*((v-37)<<23|((m<<(150-v))amp;0x007FE000))); // sign : normalized : denormalized
}
ushort float_to_half(const float x) { // IEEE-754 16-bit floating-point format (without infinity): 1-5-10, exp-15,  -131008.0,  -6.1035156E-5,  -5.9604645E-8, 3.311 digits
    const uint b = as_uint(x) 0x00001000; // round-to-nearest-even: add last bit after truncated mantissa
    const uint e = (bamp;0x7F800000)>>23; // exponent
    const uint m = bamp;0x007FFFFF; // mantissa; in line below: 0x007FF000 = 0x00800000-0x00001000 = decimal indicator flag - initial rounding
    return (bamp;0x80000000)>>16 | (e>112)*((((e-112)<<10)amp;0x7C00)|m>>13) | ((e<113)amp;(e>101))*((((0x007FF000 m)>>(125-e)) 1)>>1) | (e>143)*0x7FFF; // sign : normalized : denormalized : saturate
}
  

Пример того, как его использовать и проверить правильность преобразования:

 #include <iostream>

void print_bits(const ushort x) {
    for(int i=15; i>=0; i--) {
        cout << ((x>>i)amp;1);
        if(i==15||i==10) cout << " ";
        if(i==10) cout << "      ";
    }
    cout << endl;
}
void print_bits(const float x) {
    uint b = *(uint*)amp;x;
    for(int i=31; i>=0; i--) {
        cout << ((b>>i)amp;1);
        if(i==31||i==23) cout << " ";
        if(i==23) cout << "   ";
    }
    cout << endl;
}

int main() {
    const float x = 1.0f;
    const ushort x_compressed = float_to_half(x);
    const float x_decompressed = half_to_float(x_compressed);
    print_bits(x);
    print_bits(x_compressed);
    print_bits(x_decompressed);
    return 0;
}
  

Вывод:

 0 01111111    00000000000000000000000
0 01111       0000000000
0 01111111    00000000000000000000000
  

В этой статье я опубликовал адаптированную версию этого алгоритма преобразования FP32<->FP16 с подробным описанием того, как работает магия манипулирования битами. В этой статье я также предоставляю несколько сверхбыстрых алгоритмов преобразования для различных 16-разрядных форматов Posit.

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

Отправьте небольшой заголовок, который состоит только из минимального и максимального значений float32, затем вы можете отправить свою информацию в виде 16-битного значения интерполяции между ними. Поскольку вы также говорите, что точность не является большой проблемой, вы можете даже отправлять 8 бит за раз.

Ваше значение будет примерно таким: во время восстановления:

 float t = _t / numeric_limits<unsigned short>::max();  // With casting, naturally ;)
float val = h.min   t * (h.max - h.min);
  

Надеюсь, это поможет.

-Том

Этот вопрос уже немного устарел, но для полноты картины вы также можете взглянуть на this paper для преобразования с плавающей запятой и с плавающей запятой в половину.

Они используют подход без разветвлений, основанный на таблицах, с относительно небольшими справочными таблицами. Он полностью соответствует стандарту IEEE и даже превосходит IEEE-совместимые процедуры преобразования Phernost без ветвей по производительности (по крайней мере, на моей машине). Но, конечно, его код гораздо лучше подходит для SSE и не так подвержен эффектам задержки памяти.

Это преобразование для 16-в-32-разрядной с плавающей запятой выполняется довольно быстро для случаев, когда вам не нужно учитывать бесконечности или NAN, и может принимать денормальные значения как ноль (DAZ). Т.е. Он подходит для вычислений, зависящих от производительности, но вам следует остерегаться деления наноль, если вы ожидаете встретить денормальные значения.

Обратите внимание, что это наиболее подходит для x86 или других платформ, которые имеют условные перемещения или эквиваленты «set if».

1. Удалите знаковый бит из ввода
2. Выровняйте старший бит мантиссы по 22-му биту
3. Отрегулируйте смещение экспоненты
4. Установите биты равными нулю, если входной показатель равен нулю
5. Повторно вставьте знаковый бит

Для точности от одного до половины применяется обратное, с некоторыми дополнениями.

 void float32(float* __restrict out, const uint16_t in) {
    uint32_t t1;
    uint32_t t2;
    uint32_t t3;

    t1 = in amp; 0x7fff;                       // Non-sign bits
    t2 = in amp; 0x8000;                       // Sign bit
    t3 = in amp; 0x7c00;                       // Exponent

    t1 <<= 13;                              // Align mantissa on MSB
    t2 <<= 16;                              // Shift sign bit into position

    t1  = 0x38000000;                       // Adjust bias

    t1 = (t3 == 0 ? 0 : t1);                // Denormals-as-zero

    t1 |= t2;                               // Re-insert sign bit

    *((uint32_t*)out) = t1;
};

void float16(uint16_t* __restrict out, const float in) {
    uint32_t inu = *((uint32_t*)amp;in);
    uint32_t t1;
    uint32_t t2;
    uint32_t t3;

    t1 = inu amp; 0x7fffffff;                 // Non-sign bits
    t2 = inu amp; 0x80000000;                 // Sign bit
    t3 = inu amp; 0x7f800000;                 // Exponent

    t1 >>= 13;                             // Align mantissa on MSB
    t2 >>= 16;                             // Shift sign bit into position

    t1 -= 0x1c000;                         // Adjust bias

    t1 = (t3 > 0x38800000) ? 0 : t1;       // Flush-to-zero
    t1 = (t3 < 0x8e000000) ? 0x7bff : t1;  // Clamp-to-max
    t1 = (t3 == 0 ? 0 : t1);               // Denormals-as-zero

    t1 |= t2;                              // Re-insert sign bit

    *((uint16_t*)out) = t1;
};
  

Обратите внимание, что вы можете изменить константу 0x7bff , чтобы 0x7c00 она переполнялась до бесконечности.

Исходный код см. на GitHub.

Большинство подходов, описанных в других ответах здесь, либо неправильно округляют при преобразовании из float в half, отбрасывают субнормальные значения, что является проблемой, поскольку 2 **-14 становится вашим наименьшим ненулевым числом, либо делают неудачные вещи с Inf / NaN . Inf также является проблемой, потому что наибольшее конечное число пополам немного меньше 2 ^ 16. OpenEXR был излишне медленным и сложным, последний раз, когда я смотрел на него. Быстрый правильный подход будет использовать FPU для выполнения преобразования, либо в виде прямой инструкции, либо с использованием аппаратного округления FPU, чтобы добиться правильного результата. Любое преобразование с плавающей запятой должно быть не медленнее, чем таблица поиска 2 ^ 16 элементов.

Трудно превзойти следующие:

В OS X / iOS вы можете использовать vImageConvert_PlanarFtoPlanar16F и vImageConvert_Planar16FtoPlanarF. См. Accelerate.framework.

Intel ivybridge добавил для этого инструкции SSE. См. f16cintrin.h. Аналогичные инструкции были добавлены в ARM ISA для Neon. Смотрите vcvt_f32_f16 и vcvt_f16_f32 в arm_neon.h . На iOS вам нужно будет использовать arch arm64 или armv7s, чтобы получить к ним доступ.

Этот код преобразует 32-разрядное число с плавающей запятой в 16-разрядное и обратно.

 #include <x86intrin.h>
#include <iostream>

int main()
{
    float f32;
    unsigned short f16;
    f32 = 3.14159265358979323846;
    f16 = _cvtss_sh(f32, 0);
    std::cout << f32 << std::endl;
    f32 = _cvtsh_ss(f16);
    std::cout << f32 << std::endl;
    return 0;
}
  

Я тестировал с Intel icpc 16.0.2:

 $ icpc a.cpp
  

g 7.3.0:

 $ g   -march=native a.cpp
  

и clang 6.0.0:

 $ clang   -march=native a.cpp
  

Он печатает:

 $ ./a.out
3.14159
3.14062
  

Документация об этих встроенных свойствах доступна по адресу:

https://software.intel.com/en-us/node/524287

https://clang.llvm.org/doxygen/f16cintrin_8h.html

Вопрос старый, и на него уже был дан ответ, но я подумал, что стоит упомянуть библиотеку C с открытым исходным кодом, которая может создавать 16-битные, совместимые с IEEE, с плавающей запятой половинной точности и имеет класс, который действует практически идентично встроенному типу float, но с 16 битами вместо 32. Это «половинный» класс библиотеки OpenEXR. Код находится под разрешительной лицензией в стиле BSD. Я не верю, что у него есть какие-либо зависимости за пределами стандартной библиотеки.

У меня была точно такая же проблема, и я нашел эту ссылку очень полезной. Просто импортируйте файл «ieeehalfprecision.c» в свой проект и используйте его следующим образом :

 float myFloat = 1.24;
uint16_t resultInHalf;
singles2halfp(amp;resultInHalf, amp;myFloat, 1); // it accepts a series of floats, so use 1 to input 1 float

// an example to revert the half float back
float resultInSingle;
halfp2singles(amp;resultInSingle, amp;resultInHalf, 1);
  

Я также меняю некоторый код (см. Комментарий автора (Джеймса Турсы) По ссылке) :

 #define INT16_TYPE int16_t 
#define UINT16_TYPE uint16_t 
#define INT32_TYPE int32_t 
#define UINT32_TYPE uint32_t
  

Я нашел реализацию преобразования из формата с плавающей запятой в формат с плавающей запятой и обратно с использованием AVX2. Существует гораздо более быстрая, чем программная реализация этих алгоритмов. Я надеюсь, что это будет полезно.

преобразование 32-разрядной с плавающей запятой в 16-разрядную с плавающей запятой:

 #include <immintrin.h"

inline void Float32ToFloat16(const float * src, uint16_t * dst)
{
    _mm_storeu_si128((__m128i*)dst, _mm256_cvtps_ph(_mm256_loadu_ps(src), 0));
}

void Float32ToFloat16(const float * src, size_t size, uint16_t * dst)
{
    assert(size >= 8);

    size_t fullAlignedSize = sizeamp;~(32-1);
    size_t partialAlignedSize = sizeamp;~(8-1);

    size_t i = 0;
    for (; i < fullAlignedSize; i  = 32)
    {
        Float32ToFloat16(src   i   0, dst   i   0);
        Float32ToFloat16(src   i   8, dst   i   8);
        Float32ToFloat16(src   i   16, dst   i   16);
        Float32ToFloat16(src   i   24, dst   i   24);
    }
    for (; i < partialAlignedSize; i  = 8)
        Float32ToFloat16(src   i, dst   i);
    if(partialAlignedSize != size)
        Float32ToFloat16(src   size - 8, dst   size - 8);
}
  

преобразование 16-разрядной с плавающей запятой в 32-разрядную с плавающей запятой:

 #include <immintrin.h"

inline void Float16ToFloat32(const uint16_t * src, float * dst)
{
    _mm256_storeu_ps(dst, _mm256_cvtph_ps(_mm_loadu_si128((__m128i*)src)));
}

void Float16ToFloat32(const uint16_t * src, size_t size, float * dst)
{
    assert(size >= 8);

    size_t fullAlignedSize = sizeamp;~(32-1);
    size_t partialAlignedSize = sizeamp;~(8-1);

    size_t i = 0;
    for (; i < fullAlignedSize; i  = 32)
    {
        Float16ToFloat32<align>(src   i   0, dst   i   0);
        Float16ToFloat32<align>(src   i   8, dst   i   8);
        Float16ToFloat32<align>(src   i   16, dst   i   16);
        Float16ToFloat32<align>(src   i   24, dst   i   24);
    }
    for (; i < partialAlignedSize; i  = 8)
        Float16ToFloat32<align>(src   i, dst   i);
    if (partialAlignedSize != size)
        Float16ToFloat32<false>(src   size - 8, dst   size - 8);
}
  

Спасибо коду за десятичную точность до одинарной

На самом деле мы можем попытаться отредактировать тот же код с точностью до половины, однако это невозможно с помощью компилятора gcc C, поэтому выполните следующие действия

 sudo apt install clang
  

Затем попробуйте следующий код

 // A C code to convert Decimal value to IEEE 16-bit floating point Half precision

#include <stdio.h>

void printBinary(int n, int i)
{
 

    int k;
    for (k = i - 1; k >= 0; k--) {
 
        if ((n >> k) amp; 1)
            printf("1");
        else
            printf("0");
    }
}
 
typedef union {
    
    __fp16 f;
    struct
    {
        unsigned int mantissa : 10;
        unsigned int exponent : 5;
        unsigned int sign : 1;
 
    } raw;
} myfloat;
 

// Driver Code
int main()
{
    myfloat var;
    var.f = 11;
    printf("%d | ", var.raw.sign);
    printBinary(var.raw.exponent, 5);
    printf(" | ");
    printBinary(var.raw.mantissa, 10);
    printf("n");
    return 0;
}
  

Скомпилируйте код в вашем терминале

 clang code_name.c -o code_name
./code_name
  

Здесь

__fp16

поддерживается ли в компиляторе clang C 2-байтовый тип данных с плавающей запятой