exact fixed-point inverse square root
Nils
I'm looking for an algorithm to calculate the inverse square root of a
fixed point integer.
I know about the Newton-Raphson iteration, but this one does not work if
you're looking for an exact result because each step requires a
multiplication which in fixed point math always goes along with a loss
of precision.
Integer square roots can be calculated exact using a two bits per
iteration algorithm.
For my application the Newton-Raphson iteration does not work. Even when
I add more iterations than nessesary I always got biten by the precision
loss. The bitwise calculation of the square-root (adopted to fixed
point) followed by a division works well but is computational expensive.
I have never seen a inverse square root iteration that is not based on
the Newton-Raphson iteration and I wonder why.
Even if it's not practical for real world applications because it's is
to slow it would be nice to know how it works in principle.
I tried to come up with one of my own but failed. Any ideas?
Nils
Terje Mathisen
Isn't this sort of trivial?
A fixed-point implementation is effectively identical to an
integer-only, if you prescale the input, right?
So, by exact do you mean closest possible result, i.e. perfectly
rounded, or do you mean largest possible value such that
res * res * input_value <= 1.0
i.e. truncated rounding?
Calculate an approximate result using table lookup + NR, then square it
and multiply by the original number. The result should be very close to 1.0.
Depending upon the first step you can then adjust the first guess up or
down by one ulp, then check which guess results in the best answer.
You can also do one extra iteration of the NR, using double-wide
precision, then use the final error term to determine the correctly
rounded result.
Terje
--
- <Terje.M...@hda.hydro.com>
"almost all programming can be viewed as an exercise in caching"
glen herrmannsfeldt
> I'm looking for an algorithm to calculate the inverse square root of a
> fixed point integer.
> I know about the Newton-Raphson iteration, but this one does not work if
> you're looking for an exact result because each step requires a
> multiplication which in fixed point math always goes along with a loss
> of precision.
Fixed point multiply should be exact, but the product is twice as
long as the operands. (Well, sum if they are different lengths.)
> For my application the Newton-Raphson iteration does not work. Even when
> I add more iterations than nessesary I always got biten by the precision
> loss. The bitwise calculation of the square-root (adopted to fixed
> point) followed by a division works well but is computational expensive.
It shouldn't be hard to use N-R for fixed point square root if you
get the rounding right. For inverse square root, many values
won't have an terminating binary fraction.
> I have never seen a inverse square root iteration that is not based on
> the Newton-Raphson iteration and I wonder why.
N-R is fast, especially with a good initial guess. The floating point
ones I know of are two iterations for single precision, four for double.
(The first two can be done in single precision if it is faster.)
> Even if it's not practical for real world applications because it's is
> to slow it would be nice to know how it works in principle.
When I was first learning Fortran I wondered why no ISQRT function.
-- glen
Terje Mathisen
> Nils wrote:
>> I have never seen a inverse square root iteration that is not based on
>> the Newton-Raphson iteration and I wonder why.
>
> N-R is fast, especially with a good initial guess. The floating point
> ones I know of are two iterations for single precision, four for double.
> (The first two can be done in single precision if it is faster.)
x86 SSE has gained ~12-bit parallel lookup opcodes which return
approximate values for reciprocals and reciprocal square root.
With a 12-bit starting point you get float after a single iteration and
double with two, and as you noted even for double you can do the first
in float, allowing four (not rounded but pretty good) inverse square
roots in just a few cycles.
>> Even if it's not practical for real world applications because it's
>> is to slow it would be nice to know how it works in principle.
>
> When I was first learning Fortran I wondered why no ISQRT function.
:-)