<div dir="ltr"><div class="gmail_default" style="font-family:verdana,sans-serif"><br></div><br><div class="gmail_quote"><div dir="ltr">On Sun, Jun 10, 2018 at 10:56 PM Jonathan Wilkes <<a href="mailto:jon.w.wilkes@gmail.com">jon.w.wilkes@gmail.com</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">On Wed, Jun 6, 2018 at 9:57 PM, Matt Barber <<a href="mailto:brbrofsvl@gmail.com" target="_blank">brbrofsvl@gmail.com</a>> wrote:<br>
> Even with double precision numbers stored in those tables, you're still<br>
> restricting your input to the 254 exponent values and 10 bits of mantissa.<br>
><br>
> The rough outline of the algorithm hinges on the idea that sqrt(a*b) =<br>
> sqrt(a) * sqrt(b). So, since floating point is represented<br>
> exponent*mantissa, you can get sqrt(e*m) by doing sqrt(e)*sqrt(m). That<br>
> means if you can pre-compute all possible sqrt(e) except NaNs, infs, and<br>
> denormals, and enough mantissa values to whatever significant bits you want,<br>
> you can get really efficient. The reason to precompute the reciprocal first<br>
> is because you save the divide in the perform function. Instead, you can<br>
> rely on f*(1/sqrt(f)) = sqrt(f), so it's a simple multiply by the float you<br>
> already have.<br>
><br>
> The rest of it is just bit paperwork. All the << 23 you see are to get to<br>
> the exponent bits. The mantissa table holds rsqrt values for 1.0 + intervals<br>
> of 1./1024. (that initial 1.0 is implied in normalized float<br>
> representation).<br>
><br>
> So basically, you're accepting double-precision values for inputs truncated<br>
> to 18 bits of precision. One problem is that you'll get very bad values for<br>
> very high input exponents because you're bashing everything to<br>
> single-precision exponent range. Might not be the worst thing, but something<br>
> to consider.<br>
<br>
When PD_FLOATSIZE == 64, could we type-pun between double and uint_64<br>
and make DUMTAB1SIZE 512? (Getting the relevant bitmath right, of course.)<br><br></blockquote><div><br></div><div class="gmail_default" style="font-family:verdana,sans-serif">​Sometimes type-punning is easier with uint_32 array[2]​ than it is with uint_64, and it may also be more universal – I haven't checked about the availability of uint_64 everywhere.</div><div class="gmail_default" style="font-family:verdana,sans-serif"><br></div><div class="gmail_default" style="font-family:verdana,sans-serif">If you wanted the full range of 64-bit exponents you'd make DUMTAB1SIZE 2048 (2^11). Not sure what the best thing is for DUMTAB2SIZE – a similar proportion of bits would be 23, which would be a 64MB table that might have some cache-miss lookup overhead and I'm not certain it would actually be more efficient than just calculating sqrt(double) directly. Would need some benchmarking for sure. But if you brought that down to 20, 8MB might be a good compromise. I can think about this some more since it involves so many magic numbers and it will be annoying to have to recode it.</div><div class="gmail_default" style="font-family:verdana,sans-serif"><br></div></div></div>