Round double to float hardware-matching precision for default SSE rounding mode

double-round-to-float

constexpr uint64_t mask(int i)
{
	return (1ull << i) - 1;
}

float Round1(double tmp)
{
	uint64_t i;
	memcpy(&i, &tmp, sizeof i);

	uint32_t sgn = (i >> 32) & 0x8000'0000u;
	uint64_t exp = (i >> 52) & mask(11);
	uint64_t man = i & mask(52);

	// Handle NaN
	if (exp == 2047 && man != 0)
	{
		uint32_t fp_man = uint32_t(man >> (52 - 23));
		uint32_t fp = sgn | 0x7f80'0000 | fp_man;
		return AsFloat(fp);
	}

	// unbias exponent
	int exp_u = exp ? (int)(int64_t((exp - 1023ul) << 52) >> 52) : 0;

	uint32_t fp_exp = exp ? (unsigned char)(exp_u + 127) : 0;

	// Handle subnormal results
	if (exp_u < -126)
	{
		int denorm_shift_count = std::min(-(exp_u + 126), 25);
		man |= 1ull << 52;
		man >>= denorm_shift_count;
		// Subnormal values use exp == 0, mantissa != 0
		// If we shifted all the way to leave just zeroes in the mantissa, we'll return zero or negative zero.
		fp_exp = 0;
	}
	// Handle overflow by returning infinity
	else if (exp_u > 127)
	{
		fp_exp = 255;
		man = 0;
	}

	// Reconstitute float var
	uint32_t fp_man = uint32_t(man >> (52 - 23));

	// Round mantissa

	// Handle the case that it's a tie-breaking round
	if ((man & mask(52 - 23)) == (1 << 28))
	{
		fp_man += fp_man & 1;
	}
	// Handle the case that we need to round up
	else if (man & (1 << 28))
	{
		fp_man++;
	}

	fp_man &= (uint32_t)mask(23);

	return AsFloat(sgn | (fp_exp << 23) | fp_man);
}