Square root of s15.16 fixed point number in Java - Stack Overflow
to customize your list.
Join the Stack Overflow Community
Stack Overflow is a community of 6.7 million programmers, just like you, helping each other.
J it only takes a minute:
I want to write a function to calculate the square root of a s15.16 fixed point number. I know its a signed number with 15 digit int and 16 digit fraction. Is there anyway to do it without any libraries? Any other languages is fine too.
9,68032856
1,78832948
I assume you are asking this question because the platform you are on does not provide floating-point, otherwise you can implement 15.16 fixed-point square root via the floating-point square root as follows (this is C code, I assume Java code will look very similar):
r = (int)(sqrt (x / 65536.0) * 65536.0 + 0.5);
If your target platform provides a fast integer multiply (in particular, either a multiply with double-width result or a multiply-high instruction), and you can spare some memory for a small table, use of Newton-Raphson iterations plus a table-based starting approximation is typically the way to go. Typically, one approximates the reciprocal square root because it has a more convenient NR iteration. This gives rsqrt(x) = 1 / sqrt(x). By multiplying it with x one then gets the square root, i.e. sqrt(x) = rsqrt(x) * x. The following code shows how to compute a correctly rounded 16.16 fixed-point square root in this fashion (since the argument to the square root must be positive, this works just as well for s15.16 fixed-point). Rounding is performed by minimizing the residual x - sqrt(x)*sqrt(x).
I apologize that the square root function itself is 32-bit x86 inline assembly code but I last needed this about 10 years ago and this is all I have. I hope you can extract relevant operations from the fairly extensive comments. I included the generation of the table for the starting approximation as well as a test framework that tests the function exhaustively.
#include &stdlib.h&
#include &math.h&
unsigned int tab[96];
__declspec(naked) unsigned int __stdcall fxsqrt (unsigned int x)
edx, [esp + 4]// x
ecx, 31// 31
eax,// bsr(x)
$// if (!x) return x, avoid out-of-bounds access
// save per calling convention
// save per calling convention
ecx,// leading zeros = lz = 31 - bsr(x)
// compute table index
ecx, 0// lz & 0xfffffffe
edx,// z = x && (lz & 0xfffffffe)
edx, 25// z && 25
// retrieve initial approximation from table
edx, [tab+4*edx-128];// r = tab[(z && 25) - 32]
// first Newton-Raphson iteration
ebx, [edx*2+edx]// 3 * r
// f = (((unsigned __int64)z) * r) && 32
ebx, 22// r = (3 * r) && 22
ebx,// r = r - f
// second Newton-Raphson iteration
// prod = ((unsigned __int64)r) * z
eax,// s = prod && 32
// prod = ((unsigned __int64)r) * s
eax, 0x// 0x
eax,// s = 0x - (prod && 32)
// prod = ((unsigned __int64)r) * s
eax,// r = prod && 32
// prod = ((unsigned __int64)r) *
// restore per calling convention
// restore per calling convention
eax, [esp + 4]// x
eax, 17// x && 17
// denormalize
ecx, 1// lz && 1
edx, 3// r = (unsigned)(prod && 32); r && 3
edx,// r = (r && (lz && 1)) && 3
// remainder can be negative
ecx, [edx+edx]// 2*r
ecx,// 2*r*r
eax,// rem = (x && 17) - (2*r*r))
ecx, [edx+edx+1]// 2*r+1
ecx,// ((int)(2*r+1)) & rem))
ecx, [edx+1]// r++
edx,// if (((int)(2*r+1)) & rem) r++
eax,// result in EAX per calling convention
4// pop function argument and return
int main (void)
unsigned int i,
// build table of reciprocal square roots and their (rounded) cubes
for (i = 0; i & 96; i++) {
r = (unsigned int)(sqrt (1.0 / (1.0 + (i + 0.5) / 32.0)) * 256.0 + 0.5);
tab[i] = ((r * r * r + 4) & 0x00ffffff8) * 256 +
// exhaustive test of 16.16 fixed-point square root
r = (unsigned int)(sqrt (i / 65536.0) * 65536.0 + 0.5);
if (r != fxsqrt (i)) {
printf ("error @ %08x: ref = %08x
res=%08x\n", i, r, fxsqrt (i));
} while (i);
9,68032856
Use , with the simple observation that √(2-16a) = 2-8√a.
103k19130234
Your Answer
Sign up or
Sign up using Google
Sign up using Facebook
Sign up using Email and Password
Post as a guest
Post as a guest
By posting your answer, you agree to the
Not the answer you're looking for?
Browse other questions tagged
rev .24947
Stack Overflow works best with JavaScript enabled