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

[Jonathan Wilkes]

On Mon, Jun 11, 2018 at 9:00 AM, Matt Barber <brbrofsvl at gmail.com> wrote:
>
>
> On Sun, Jun 10, 2018 at 10:56 PM Jonathan Wilkes <jon.w.wilkes at gmail.com>
> wrote:
>>
>> 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.)
>>
>
> 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.
>
> 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.

I forgot about sqrt inside [expr].

So what if we keep the current algo with PD_FLOATSIZE == 64? That should be
fast for [sqrt~] at both t_float sizes, with an escape hatch of [expr]
for higher
precision math. (At least it appears [expr] will compute sqrt at
double precision
by default, and it seems to use t_float properly.)

Then if a user comes up with a case that requires both higher precision and
better performance than [expr~] can provide, we can create a benchmark patch
and start trying out various magic numbers for those tables.

-Jonathan

>
>
> _______________________________________________
> L2Ork-dev mailing list
> L2Ork-dev at disis.music.vt.edu
> https://disis.music.vt.edu/listinfo/l2ork-dev