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

[Jonathan Wilkes]

On Mon, Jun 11, 2018 at 10:38 AM, Matt Barber <brbrofsvl at gmail.com> wrote:
>
>
> On Mon, Jun 11, 2018 at 12:43 PM Jonathan Wilkes <jon.w.wilkes at gmail.com>
> wrote:
>>
>> 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
>>
>
> If it were up to me I'd try to do the analogous thing for doubles since the
> rest of the math in double precision expands in all the other objects and we
> want to keep it as seamless as possible. Otherwise, if it's really
> documented well that these particular objects have only float precision in
> single- and double-precision PD, it could be not too bad. But if we wanted
> to go for it, bit precision in the exponents is going to be most important.
> 16-20 bits of mantissa precision (0.5-8MB table) would probably work out OK,
> but it would be worth testing.

Quick sanity check--  from d_math.c in both Purr Data and Pd Vanilla:

/* sigrsqrt - reciprocal square root good to 8 mantissa bits  */

That should read "good to 10 mantissa bits," right?

right-shift by 23 plus bitmask to filter the sign gets us the 8 exponent bits.

right-shift by 13 plus bitmask to filter out the exponent/sign bits
gets us the 10 most
significant mantissa bits.

-Jonathan