ProjectPhysX
commented
2 years ago Hi all,
I'm researching Posits as a potential storage format for the DDFs in the lattice Boltzmann method. Therefore I needed much faster conversion algorithms than previously available, especially to be suitable for GPUs. I wrote these conversion algorithms using just basic bit manipulation and to be branchless. They are much shorter than the existing SoftPosit conversion, and also much faster. Tests with them integrated into the LBM showed 22x performance gain compared to the existing SoftPosit implementation.
Below is the validation plot, showing digits = -log10(fabs(log10(x_converted_back_and_forth/x))), proving the conversion is bitwise identical and handles rounding/overflow correctly. Note that I have ditched the infinity and NaR definitions to spare some operations.
The purpose of these algorithms really is to be as fast as possible for applications where lots of conversions are needed, for example compressing numbers for storage in memory.
Regards, Moritz Lehmann
typedef unsigned int uint;
typedef unsigned short ushort;
float as_float(const uint x) {
return *(float*)&x;
}
uint as_uint(const float x) {
return *(uint*)&x;
}
int as_int(const float x) {
return *(int*)&x;
}
ushort float_to_p160(const float x) {
const uint b = as_uint(x);
const int e = ((b&0x7F800000)>>23)-127; // exponent-bias
int m = (b&0x007FFFFF)>>9; // mantissa
const int v = abs(e); // shift
const int r = (e<0 ? 0x0002 : 0xFFFE)<<(13-v); // generate regime bits
m = ((m>>(v-(e<0)))+1+(e<-13)*0x2)>>1; // rounding: add 1 after truncated position; in case of lowest numbers, saturate
return (b&0x80000000)>>16 | (e>-16)*((r+m)&0x7FFF) | (e>13)*0x7FFF; // sign | regime+mantissa ("+" handles rounding overflow) | saturate
}
float p160_to_float(const ushort x) {
const uint sr = (x>>14)&1; // sign of regime
ushort t = x<<2; // remove sign and first regime bit
t = sr ? ~t : t; // positive regime r>=0 : negative regime r<0
const int r = 142-(as_int((float)t)>>23); // evil log2 bit hack to count leading zeros for regime
const uint m = (x<<(r+10))&0x007FFFFF; // extract mantissa and bit-shift it in place
const int rs = sr ? r : -r-1; // negative regime r<0 : positive regime r>=0
return as_float((x&0x8000)<<16 | (r!=158)*((rs+127)<<23 | m)); // sign | regime | mantissa
}
ushort float_to_p161(const float x) {
const uint b = as_uint(x);
const int e = ((b&0x7F800000)>>23)-127; // exponent-bias
int m = (b&0x007FFFFF)>>10; // mantissa
const int ae = abs(e);
const int v = ae>>1; // shift, ">>1" is the same as "/2"
const int e2 = ae&1; // "&1" is the same as "%2"
const int r = ((e<0 ? 0x0002 : 0xFFFE<<e2)+e2)<<(13-v-e2); // generate regime bits, merge regime+exponent and shift in place
m = ((m>>(v-(e<0)*(1-e2)))+(e>-28)+(e<-26)*0x3)>>1; // rounding: add 1 after truncated position; in case of lowest numbers, saturate
return (b&0x80000000)>>16 | (e>-31)*((r+m)&0x7FFF) | (e>26)*0x7FFF; // sign | regime+exponent+mantissa ("+" handles rounding overflow) | saturate
}
float p161_to_float(const ushort x) {
const uint sr = (x>>14)&1; // sign of regime
ushort t = x<<2; // remove sign and first regime bit
t = sr ? ~t : t; // positive regime r>=0 : negative regime r<0
const int r = 142-(as_int((float)t)>>23); // evil log2 bit hack to count leading zeros for regime
const uint e = (x>>(12-r))&1; // extract mantissa and bit-shift it in place
const uint m = (x<<(r+11))&0x007FFFFF; // extract mantissa and bit-shift it in place
const int rs = (sr ? r : -r-1)<<1; // negative regime r<0 : positive regime r>=0, "<<1" is the same as "*2"
return as_float((x&0x8000)<<16 | (r!=158)*((rs+e+127)<<23 | m)); // sign | regime+exponent | mantissa
}EDIT: Flush to zero now works correctly. EDIT 2: simplified evil log2 bit hack part
milankl
https://twitter.com/ProjectPhysX/status/1444007905058070531?s=20 suggests a new conversion between Float32 and Posit16