Do you have any recommendations on single-precision floating point libraries?

taby said:
but neither emulate the single-precision float

What? In what situation would you have to emulate single-precision floats? Are you working on an embedded microcontroller or something?

You are aware that single-precision is 32-bit float, double is 64-bit float, and half is 16-bit float?

taby said:
minimum precision is 16-bits instead of 8

Am I correctly interpreting the tone of this quote in that you are lamenting that you don't get 8-bit floats?

Or are you actually referring to fixedpoint math?

Also, what is the point of doing calculations in double when you then reduce them into floats? You just add lossy conversion operations to it, but you still perform double-sized instructions internally, which is not the same as single-precision float maths. All you're doing is truncating results, with no performance or efficiency gain except for the fact that the results are now half as large in memory.

Also, if you want to be efficient / performant, you should not use if/else if branches to clamp floats, except if you know that those branches are extremely rarely taken. It would be better to use branchless code and also to use SIMD code while you're at it.

What is even the goal here? Performance gains? Memory compression? Adhering to some IEEE specification without hardware variance or something?

And finally: what function do you even want to compute? How could we possibly recommend a floating point math library if you don't tell us what you need? Do you need statistical functions, numerical integrators, trigonometry, quaternions?

Thank you for your input.

Claude came up with a code that works. Truncating to float from double snaps the value, further quantizing it. Now that I have code from Claude, I can experiment to see which exponent and mantissa bit combinations provide the best answer, and ultimately why.\\

Claude's code is:

struct FloatFormat {
	int expBits;   // E: number of exponent bits  (1 .. 16 is sensible)
	int mantBits;  // M: number of stored mantissa bits, 1 + E + M <= 64
};

// ---------------------------------------------------------------------------
// encode: double -> bit pattern of the custom format (right-aligned in the
// returned uint64_t).  Rounding mode: round to nearest, ties to even.
// ---------------------------------------------------------------------------
std::uint64_t encode(double value, FloatFormat fmt)
{
	const int E = fmt.expBits;
	const int M = fmt.mantBits;
	const std::uint64_t mantMask = (std::uint64_t(1) << M) - 1;
	const std::uint64_t maxExpField = (std::uint64_t(1) << E) - 1;
	const std::int64_t  bias = (std::int64_t(1) << (E - 1)) - 1;

	// --- pick the double apart -------------------------------------------
	std::uint64_t in;
	std::memcpy(∈, &value, sizeof in);

	const std::uint64_t sign = in >> 63;
	const std::int64_t  dExp = static_cast<std::int64_t>((in >> 52) & 0x7FF);
	const std::uint64_t dMan = in & ((std::uint64_t(1) << 52) - 1);

	const std::uint64_t signOut = sign << (E + M);

	// --- special values -----------------------------------------------------
	if (dExp == 0x7FF) {
		if (dMan != 0)                                   // NaN -> quiet NaN
			return signOut | (maxExpField << M) | (std::uint64_t(1) << (M - 1));
		return signOut | (maxExpField << M);             // +/- infinity
	}
	if (dExp == 0 && dMan == 0)                          // +/- zero
		return signOut;

	// --- normalise into  value = (-1)^sign * sig * 2^(e - 52) ---------------
	// with sig having its leading 1 in bit 52.
	std::uint64_t sig;
	std::int64_t  e;
	if (dExp == 0) {                       // subnormal double
		sig = dMan;
		e = 1 - 1023;
		while (!(sig & (std::uint64_t(1) << 52))) { sig <<= 1; --e; }
	}
	else {                               // normal double
		sig = dMan | (std::uint64_t(1) << 52);
		e = dExp - 1023;
	}

	// --- place into the target format ---------------------------------------
	std::int64_t expField = e + bias;      // tentative biased exponent
	std::int64_t shift = 52 - M;        // bits to discard from sig

	if (expField <= 0) {                   // target subnormal (or underflow):
		shift += 1 - expField;           // discard additional bits
		expField = 0;
	}

	// --- round to nearest, ties to even --------------------------------------
	std::uint64_t kept, roundBit, sticky;
	if (shift >= 64) {                     // everything rounds away
		kept = 0; roundBit = 0; sticky = (sig != 0);
	}
	else if (shift <= 0) {               // target mantissa wider than source
		kept = sig << (-shift); roundBit = 0; sticky = 0;
	}
	else {
		kept = sig >> shift;
		const std::uint64_t rem = sig & ((std::uint64_t(1) << shift) - 1);
		roundBit = (rem >> (shift - 1)) & 1;
		sticky = (rem & ((std::uint64_t(1) << (shift - 1)) - 1)) != 0;
	}
	if (roundBit && (sticky || (kept & 1)))
		++kept;                            // round up

	// --- propagate a possible carry out of the mantissa ----------------------
	if (expField == 0) {
		// kept may have become exactly 2^M: smallest normal number.
		if (kept >> M) { expField = 1; kept &= mantMask; }
	}
	else {
		// kept holds the implicit bit at position M; rounding may have
		// carried it to position M+1 (e.g. 1.111...1 -> 10.000...0).
		if (kept >> (M + 1)) { ++expField; kept >>= 1; }
		kept &= mantMask;                  // drop the implicit leading 1
	}

	// --- overflow -> infinity -------------------------------------------------
	if (expField >= static_cast<std::int64_t>(maxExpField))
		return signOut | (maxExpField << M);

	return signOut
		| (static_cast<std::uint64_t>(expField) << M)
		| kept;
}

// ---------------------------------------------------------------------------
// decode: bit pattern of the custom format -> double (exact, since any
// format with M <= 52 and E <= 11-ish is a subset of binary64's range;
// for wider formats the conversion itself rounds via ldexp).
// ---------------------------------------------------------------------------
double decode(std::uint64_t bits, FloatFormat fmt)
{
	const int E = fmt.expBits;
	const int M = fmt.mantBits;
	const std::uint64_t maxExpField = (std::uint64_t(1) << E) - 1;
	const std::int64_t  bias = (std::int64_t(1) << (E - 1)) - 1;

	const std::uint64_t sign = (bits >> (E + M)) & 1;
	const std::uint64_t expField = (bits >> M) & maxExpField;
	const std::uint64_t man = bits & ((std::uint64_t(1) << M) - 1);

	const double s = sign ? -1.0 : 1.0;

	if (expField == maxExpField)
		return man ? std::numeric_limits<double>::quiet_NaN()
		: s * std::numeric_limits<double>::infinity();

	if (expField == 0)                       // zero or subnormal
		return s * std::ldexp(static_cast<double>(man),
			static_cast<int>(1 - bias - M));

	return s * std::ldexp(static_cast<double>(man | (std::uint64_t(1) << M)),
		static_cast<int>(expField - bias - M));
}

double truncate_normalized_double(double d)
{
	if (d < 0.0)
	{
		return 0.0;
	}
	else if (d > 1.0)
	{
		return 1.0;
	}

	const FloatFormat precision{ 8, 23 };
	const std::uint64_t b = encode(d, precision);
	return decode(b, precision);
}	

So your goal is to have a more compact representation of floating point numbers? Is that compression just intended for storage / transmission, or is it absolutely required that you can also do math directly on the compressed version? And is the range you want to represent known? What's the precision requirements for different numeric ranges? And do you really care about being able to represent infinities and NaNs?

Your encoding looks super slow with all the if/else and STL calls (which are also notoriously slow).

Wait… what? You make a generalised double-precision floating point decomposer and composer, handling even cases like infinity and NaN, but you only call it on values from 0 to 1?

And why do you reduce doubles in the range from 0 to 1 to 8+23? That's literally the same as just casting the double to a float, just dozens of times slower. Because a float is 1s+8e+23m on all platforms.

And then you return that as a double again anyway. What could possibly be the benefit of all this?

I am further quantizing the value of beta, which produces interesting results. It’s just a coincidence though.

taby wrote:

I am further quantizing the value of beta, which produces interesting results. It’s just a coincidence though.

So, you are still using rounding errors to get closer to the desired results?

... you will never figure out how gravity works.

But maybe you can confirm that simulation theory nonsense is actually true. ;D

Like I said, it’s just a coincidence. :)