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

[Matt Barber]

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.
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://disis.music.vt.edu/pipermail/l2ork-dev/attachments/20180611/0048745f/attachment.html>