Add pure FixedPoint sqrt implementation and tests

[bytegrrrl]

This implements a pure fixed point square root function per #1543.

Notes:
I've just replaced the existing function with a pure fixed point implementation (can't overload when only return type differs). It should be trivial to make it return an actual FixedPoint instead, but that would requiring updating prior test cases and I wasn't sure if that was okay to do.

There's a small loss in precision depending on the value of fractional_bits and the position of the most significant bit: if the integer portion is very large, we won't have as much (absolute) precision. Ideally you would want the intermediate_type to be twice the size of raw_type to avoid any losses.

[heinezen]

Very nice, thanks a lot :)

I have a few comments about the code, but nothing too serious.

[heinezen]

Same here, I think the std::numeric_limits definitions are more readable. Also, shifting by one is just multiplication by 2 and the compiler can probably optimize that out on its own :D

Suggested change

	for (uint64_t i = 1; i < (1ul << 63); i = (i << 1) ^ i) {
	for (uint64_t i = 1; i < std::numeric_limits<uint64_t>::max(); i = (i * 2) ^ i) {

[bytegrrrl]

Thank you for all the feedback! I think I've incorporated everything. If this all looks good, I'd be happy to work on some of the other math functions too. (Hopefully I'll be able to break my C habits!)

---

I ran into an interesting bug / edge case withTESTEQUALS_FLOAT by the way:

An unsigned fixed point 0.0 and a double 0.0 fail the equality check.
So, part of the macro tests for inequality with lhs < (rhs - epsilon). The result of rhs - epsilon is static_cast back into its raw type for the actual comparison. Now, if rhs - epsilon is negative and we cast into an unsigned type, then this behavior is undefined. On x86, the value underflows so the inequality condition is met. On arm, the value is set to 0. (Which was why I didn't catch the failing test).

Here's an example case that fails on x86:

TESTEQUALS_FLOAT((FixedPoint<uint32_t, 16>::zero()), 0.0, 1e-4);

I think it should be possible to resolve this by checking for lhs + epsilon < rhs instead. But floats are cursed and aren't associative, so this could cause unexpected behavior down the line. I've chosen to just omit the 0 case for now.

[heinezen]

I found a few new issues

[heinezen]

I have a very small suggestion but everything else looks nice.

[heinezen]

Thanks again!

[bytegrrrl]

Thank you for all the guidance! I might take a look at some of the others soon.