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

[Matt Barber]

Oops, that's 8th and 9th

On Mon, Jun 11, 2018, 4:38 PM Matt Barber <brbrofsvl at gmail.com> wrote:

> Good to 8 because of the sqrt operation. If you compare results of sqrt
> with 10 mantissa bits on and then with 23 mantissa bits on, they will
> differ at the 7th bit, or 8th when you count the implied bit.
>
> On Mon, Jun 11, 2018, 3:27 PM Jonathan Wilkes <jon.w.wilkes at gmail.com>
> wrote:
>
>> 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
>> _______________________________________________
>> L2Ork-dev mailing list
>> L2Ork-dev at disis.music.vt.edu
>> https://disis.music.vt.edu/listinfo/l2ork-dev
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://disis.music.vt.edu/pipermail/l2ork-dev/attachments/20180611/f36cb6bd/attachment-0001.html>