[L2Ork-dev] Precision of [sqrt~] for double-precision

[Jonathan Wilkes]

On Wed, Jun 6, 2018 at 9:57 PM, Matt Barber <brbrofsvl at gmail.com> wrote:
> Even with double precision numbers stored in those tables, you're still
> restricting your input to the 254 exponent values and 10 bits of mantissa.
>
> The rough outline of the algorithm hinges on the idea that sqrt(a*b) =
> sqrt(a) * sqrt(b). So, since floating point is represented
> exponent*mantissa, you can get sqrt(e*m) by doing sqrt(e)*sqrt(m). That
> means if you can pre-compute all possible sqrt(e) except NaNs, infs, and
> denormals, and enough mantissa values to whatever significant bits you want,
> you can get really efficient. The reason to precompute the reciprocal first
> is because you save the divide in the perform function. Instead, you can
> rely on f*(1/sqrt(f)) = sqrt(f), so it's a simple multiply by the float you
> already have.
>
> The rest of it is just bit paperwork. All the << 23 you see are to get to
> the exponent bits. The mantissa table holds rsqrt values for 1.0 + intervals
> of 1./1024. (that initial 1.0 is implied in normalized float
> representation).
>
> So basically, you're accepting double-precision values for inputs truncated
> to 18 bits of precision. One problem is that you'll get very bad values for
> very high input exponents because you're bashing everything to
> single-precision exponent range. Might not be the worst thing, but something
> to consider.

When PD_FLOATSIZE == 64, could we type-pun between double and uint_64
and make DUMTAB1SIZE 512? (Getting the relevant bitmath right, of course.)

-Jonathan