Точное преобразование 32-разрядного целого числа без знака в число с плавающей точкой в ​​диапазоне (-1; 1)

Question

Согласно подобным статьям , половина чисел с плавающей точкой находится в интервале [-1,1].Не могли бы вы предложить, как использовать этот факт, чтобы заменить наивное преобразование 32-разрядного целого числа без знака в число с плавающей запятой (при сохранении равномерного распределения)?

Наивный код:

uint32_t i = /* randomly generated */;
float f = (float)i / (1ui32<<31) - 1.0f;

Проблема здесь в том, что сначала число i преобразуется в float, теряя до 8 младших битов точности.Только тогда число масштабируется до [0; 2) интервала, а затем до [-1; 1) интервала.

Пожалуйста, предложите решение на C или C ++ для CPU x86_64 или CUDA, если вы его знаете.

Обновление: решение с double хорошо для x86_64, но слишком медленно в CUDA.Извините, я не ожидал такого ответа.Любые идеи, как этого добиться без использования чисел с плавающей запятой двойной точности?

Answer 1

Вы можете выполнить вычисление, используя вместо этого double, чтобы не потерять точность значения uint32_t, а затем присвоить результат float.

.

float f = (double)i / (1ui32<<31) - 1.0;

Answer 2

Необходимо понять, что в то время как (float)i теряет 8-битную точность (так что она имеет 24-битную точность), результат также имеет только 24-битную точность. Таким образом, эта потеря точности не обязательно является плохой вещью (это на самом деле более сложно, потому что если i меньше, он потеряет менее 8 бит. Но все будет хорошо).

Так что нам просто нужно исправить диапазон, чтобы изначально неотрицательное значение было сопоставлено с INT_MIN..INT_MAX.

Это выражение работает: (float)(int)(value^0x80000000)/0x80000000.

Вот как это работает:

1. Часть (int)(value^0x80000000) переворачивает бит знака, поэтому 0x0 отображается на INT_MIN, а 0xffffffff отображается на INT_MAX.
2. Затем происходит преобразование в float. Здесь происходит некоторое округление, и мы теряем точность (но это не проблема).
3. Затем просто разделите на 0x80000000, чтобы попасть в диапазон [-1..1]. Поскольку это деление просто корректирует экспонентную часть, это деление не теряет никакой точности.

Итак, есть только одно округление, остальные операции не теряют точности. Эти цепочки операций должны иметь тот же эффект, что и вычисление результата с бесконечной точностью, затем округление до float (это теоретическое округление имеет тот же эффект, что и округление на шаге 2.)

Но, чтобы быть абсолютно уверенным, я с помощью грубой силы проверил все 32-битные значения, что это выражение приводит к тому же значению, что и (float)((double)value/0x80000000-1.0).

Answer 3

В случае, если вы отбросите ограничение равномерного распределения, оно выполнимо только для 32-битной целочисленной арифметики:

//---------------------------------------------------------------------------
float i32_to_f32(int   x)
    {
    int exp;
    union _f32          // semi result
        {
        float f;        // 32bit floating point
        DWORD u;        // 32 bit uint
        } y;
    // edge cases
    if (x== 0x00000000) return  0.0f;
    if (x< -0x1FFFFFFF) return -1.0f;
    if (x> +0x1FFFFFFF) return +1.0f;
    // conversion
    y.u=0;                              // reset bits
    if (x<0){ y.u|=0x80000000; x=-x; }  // sign (31 bits left)
    exp=((x>>23)&63)-64;                // upper 6 bits -> exponent -1,...,-64 (not 7bits to avoid denormalized numbers)
    y.u|=(exp+127)<<23;                 // exponent bias and bit position
    y.u|=x&0x007FFFFF;                  // mantissa
    return y.f;
    }
//---------------------------------------------------------------------------
int f32_to_i32(float x)
    {
    int exp,man,i;
    union _f32          // semi result
        {
        float f;        // 32bit floating point
        DWORD u;        // 32 bit uint
        } y;
    // edge cases
    if (x== 0.0f) return  0x00000000;
    if (x<=-1.0f) return -0x1FFFFFFF;
    if (x>=+1.0f) return +0x1FFFFFFF;
    // conversion
    y.f=x;
    exp=(y.u>>23)&255; exp-=127;        // exponent bias and bit position
    if (exp<-64) return 0.0f;
    man=y.u&0x007FFFFF;                 // mantissa
    i =(exp<<23)&0x1F800000;
    i|= man;
    if (y.u>=0x80000000) i=-i;          // sign
    return i;
    }
//---------------------------------------------------------------------------

Я решил использовать только 29 бит + знак = ~ 30 бит целого числа, чтобы избежать хаоса денормализованных чисел, который мне лень кодировать (это даст вам 30 или даже 31 бит, но гораздо медленнее и сложнее).

Но распределение не является ни линейным, ни равномерным:

в красном - это float в диапазоне <-1,+1>, а в синем * integer в диапазоне <-1FFFFFFF,+1FFFFFFF>.

С другой стороны, в обоих преобразованиях вообще нет округления ...

PS. Я думаю, что может быть способ несколько линеаризовать результат с использованием предварительно вычисленного LUT для 6-битного показателя (64 значения).

Answer 4

Есть какие-нибудь идеи, как этого добиться без использования чисел с плавающей запятой двойной точности?

Не вдаваясь в подробности о float:

Shift u до тех пор, пока не будет установлен старший значащий бит, значение преобразования float уменьшается вдвое.

«Сохраняя равномерное распределение»

50% значений uint32_t будут находиться в [0,5... 1,0)
25% от значений uint32_t будет в [0,25 ... 0,5)
12,5% от значений uint32_t будет в [0,125 ... 0,25)
6,25% от значений uint32_t будут находиться в [0,0625 ... 0,125)
...

#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>

float ui32to0to1(uint32_t u) {
  if (u) {
    float band = 1.0f/(1llu<<32);
    while ((u & 0x80000000) == 0) {
      u <<= 1;
      band *= 0.5f;
    }
    return (float)u * band;
  }
  return 0.0f;
}

Некотором тестовом коде, чтобы показать функциональную эквивалентность double.

int test(uint32_t u) {
  volatile float f0 = (float) ((double)u / (1llu<<32));
  volatile float f1 = ui32to0to1(u);
  if (f0 != f1) {
    printf("%8lX %.7e %.7e\n", (unsigned long) u, f0, f1);
    return 1;
  }
  return 0;
}

int main(void) {
  for (int i=0; i<100000000; i++) {
    test(rand()*65535u ^ rand());
  }
  return 0;
}

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

Для повышения эффективности цикл должен только повторяться с 32 до FLT_MANT_DIG, что обычно составляет 24 *. 1033 *

float ui32to0to1(uint32_t u) {
  float band = 1.0f/(1llu<<32);
  for (int i = 32; (i>FLT_MANT_DIG && ((u & 0x80000000) == 0)); i--) {
    u <<= 1;
    band *= 0.5f;
  }
  return (float)u * band;
}

---

Этот ответ сопоставляет [от 0 до 2 ³² -1] с [0,0 до 1,0)

Для сопоставления с [от 0 до 2 ³² -1] до (-1,0-1,0).Может составлять -0,0.

if (u >= 0x80000000) {
  return ui32to0to1((u - 0x80000000)*2);
} else
  return -ui32to0to1((0x7FFFFFFF - u)*2);
}

Answer 5

Я предлагаю (если вы хотите избежать деления и использовать точное начальное значение с плавающей запятой, равное 1,0 * 2 ^ -32):

float e = i * ldexp(1.0,-32) - 1.0;