<div dir="auto">Oops, that's 8th and 9th</div><br><div class="gmail_quote"><div dir="ltr">On Mon, Jun 11, 2018, 4:38 PM Matt Barber <<a href="mailto:brbrofsvl@gmail.com">brbrofsvl@gmail.com</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="auto">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.</div><br><div class="gmail_quote"><div dir="ltr">On Mon, Jun 11, 2018, 3:27 PM Jonathan Wilkes <<a href="mailto:jon.w.wilkes@gmail.com" target="_blank" rel="noreferrer">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 Mon, Jun 11, 2018 at 10:38 AM, Matt Barber <<a href="mailto:brbrofsvl@gmail.com" rel="noreferrer noreferrer" target="_blank">brbrofsvl@gmail.com</a>> wrote:<br>
><br>
><br>
> On Mon, Jun 11, 2018 at 12:43 PM Jonathan Wilkes <<a href="mailto:jon.w.wilkes@gmail.com" rel="noreferrer noreferrer" target="_blank">jon.w.wilkes@gmail.com</a>><br>
> wrote:<br>
>><br>
>> On Mon, Jun 11, 2018 at 9:00 AM, Matt Barber <<a href="mailto:brbrofsvl@gmail.com" rel="noreferrer noreferrer" target="_blank">brbrofsvl@gmail.com</a>> wrote:<br>
>> ><br>
>> ><br>
>> > On Sun, Jun 10, 2018 at 10:56 PM Jonathan Wilkes<br>
>> > <<a href="mailto:jon.w.wilkes@gmail.com" rel="noreferrer noreferrer" target="_blank">jon.w.wilkes@gmail.com</a>><br>
>> > wrote:<br>
>> >><br>
>> >> On Wed, Jun 6, 2018 at 9:57 PM, Matt Barber <<a href="mailto:brbrofsvl@gmail.com" rel="noreferrer noreferrer" target="_blank">brbrofsvl@gmail.com</a>><br>
>> >> wrote:<br>
>> >> > Even with double precision numbers stored in those tables, you're<br>
>> >> > still<br>
>> >> > restricting your input to the 254 exponent values and 10 bits of<br>
>> >> > mantissa.<br>
>> >> ><br>
>> >> > The rough outline of the algorithm hinges on the idea that sqrt(a*b)<br>
>> >> > =<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).<br>
>> >> > That<br>
>> >> > means if you can pre-compute all possible sqrt(e) except NaNs, infs,<br>
>> >> > and<br>
>> >> > denormals, and enough mantissa values to whatever significant bits<br>
>> >> > you<br>
>> >> > want,<br>
>> >> > you can get really efficient. The reason to precompute the reciprocal<br>
>> >> > first<br>
>> >> > is because you save the divide in the perform function. Instead, you<br>
>> >> > can<br>
>> >> > rely on f*(1/sqrt(f)) = sqrt(f), so it's a simple multiply by the<br>
>> >> > float<br>
>> >> > you<br>
>> >> > already have.<br>
>> >> ><br>
>> >> > The rest of it is just bit paperwork. All the << 23 you see are to<br>
>> >> > get<br>
>> >> > to<br>
>> >> > the exponent bits. The mantissa table holds rsqrt values for 1.0 +<br>
>> >> > 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<br>
>> >> > truncated<br>
>> >> > to 18 bits of precision. One problem is that you'll get very bad<br>
>> >> > values<br>
>> >> > for<br>
>> >> > very high input exponents because you're bashing everything to<br>
>> >> > single-precision exponent range. Might not be the worst thing, but<br>
>> >> > 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<br>
>> >> course.)<br>
>> >><br>
>> ><br>
>> > Sometimes type-punning is easier with uint_32 array[2] than it is with<br>
>> > uint_64, and it may also be more universal – I haven't checked about the<br>
>> > availability of uint_64 everywhere.<br>
>> ><br>
>> > If you wanted the full range of 64-bit exponents you'd make DUMTAB1SIZE<br>
>> > 2048<br>
>> > (2^11).<br>
>> > Not sure what the best thing is for DUMTAB2SIZE – a similar<br>
>> > proportion of bits would be 23, which would be a 64MB table that might<br>
>> > have<br>
>> > some cache-miss lookup overhead and I'm not certain it would actually be<br>
>> > more efficient than just calculating sqrt(double) directly. Would need<br>
>> > some<br>
>> > benchmarking for sure. But if you brought that down to 20, 8MB might be<br>
>> > a<br>
>> > good compromise. I can think about this some more since it involves so<br>
>> > many<br>
>> > magic numbers and it will be annoying to have to recode it.<br>
>><br>
>> I forgot about sqrt inside [expr].<br>
>><br>
>> So what if we keep the current algo with PD_FLOATSIZE == 64? That should<br>
>> be<br>
>> fast for [sqrt~] at both t_float sizes, with an escape hatch of [expr]<br>
>> for higher<br>
>> precision math. (At least it appears [expr] will compute sqrt at<br>
>> double precision<br>
>> by default, and it seems to use t_float properly.)<br>
>><br>
>> Then if a user comes up with a case that requires both higher precision<br>
>> and<br>
>> better performance than [expr~] can provide, we can create a benchmark<br>
>> patch<br>
>> and start trying out various magic numbers for those tables.<br>
>><br>
>> -Jonathan<br>
>><br>
><br>
> If it were up to me I'd try to do the analogous thing for doubles since the<br>
> rest of the math in double precision expands in all the other objects and we<br>
> want to keep it as seamless as possible. Otherwise, if it's really<br>
> documented well that these particular objects have only float precision in<br>
> single- and double-precision PD, it could be not too bad. But if we wanted<br>
> to go for it, bit precision in the exponents is going to be most important.<br>
> 16-20 bits of mantissa precision (0.5-8MB table) would probably work out OK,<br>
> but it would be worth testing.<br>
<br>
Quick sanity check--  from d_math.c in both Purr Data and Pd Vanilla:<br>
<br>
/* sigrsqrt - reciprocal square root good to 8 mantissa bits  */<br>
<br>
That should read "good to 10 mantissa bits," right?<br>
<br>
right-shift by 23 plus bitmask to filter the sign gets us the 8 exponent bits.<br>
<br>
right-shift by 13 plus bitmask to filter out the exponent/sign bits<br>
gets us the 10 most<br>
significant mantissa bits.<br>
<br>
-Jonathan<br>
_______________________________________________<br>
L2Ork-dev mailing list<br>
<a href="mailto:L2Ork-dev@disis.music.vt.edu" rel="noreferrer noreferrer" target="_blank">L2Ork-dev@disis.music.vt.edu</a><br>
<a href="https://disis.music.vt.edu/listinfo/l2ork-dev" rel="noreferrer noreferrer noreferrer" target="_blank">https://disis.music.vt.edu/listinfo/l2ork-dev</a></blockquote></div>
</blockquote></div>