My understanding of your code is that you've got `struct Number { int16 a,b,c,d; };` and the value of a number is equal to the floating point number 256.0a + b + c / 256.0 + d / 65536.0. (if the arch had floats, which it doesn't) If this assumption is wrong, my bad.
What I'm proposing is that you do `b = i16.sqrt(256a + b); a = 0;` then do 1-2 iterations of newton's method. For instance, if you have the value 1234.56789 then i16.sqrt(256*a+b) = i16.sqrt(1234) = 35. So your first approximation is x = 35.56789. Do x = (x + 1234.56789/x)/2, and with floats you get x = 35.139035. Do it again, you get x = 35.136418. The true value is also 35.136418, so there ya go.
> I'd be concerned about Newton's method (specifically the squaring of x_n) overflowing the fixed point.
I may have lost the plot. Are you calculating _inverse_ square root? x_n isn't squared for Newton's method of square root, but it is for inverse square root.
(Either way, the square of the square root of a number should be the same as the original number. So if you have overflow, you have other problems.)
Unrelated: What do you need pi and trig functions for? I wrote a ray tracer (actually a path tracer, same diff) about a year ago, and I only used trig/pi to set up the camera for the render. While the render was actually rendering, no trig was happening. It didn't matter for me, so I just used the standard library, but in your case it would probably be fine to use slow brutish methods.