Java中s15.16定点数的平方根
我想编写一个函数来计算s15.16定点数的平方根。 我知道它有一个15位数字和16位数分数的带符号数字。 无论如何没有任何图书馆吗？ 任何其他语言也没关系。
我假设你问这个问题，因为你所在的平台没有提供浮点，否则你可以通过浮点平方根实现15.16定点平方根，如下所示（这是C代码，我假设Java代码将看起来非常相似）：
int x, r; r = (int)(sqrt (x / 65536.0) * 65536.0 + 0.5); 如果您的目标平台提供快速整数乘法（特别是乘以双倍宽度结果或乘高指令），并且您可以为小型表节省一些内存，则使用Newton-Raphson迭代加上表格 – 基于起始近似通常是要走的路。 通常，一个近似倒数平方根，因为它具有更方便的NR迭代。 这给出了rsqrt（x）= 1 / sqrt（x）。 通过将其乘以x 1，然后得到平方根，即sqrt（x）= rsqrt（x）* x。 下面的代码显示了如何以这种方式计算正确舍入的16.16定点平方根（因为平方根的参数必须为正，这对于s15.16定点也同样如此）。 通过最小化残差x-sqrt（x）* sqrt（x）来执行舍入。
我很抱歉，平方根函数本身就是32位x86内联汇编代码，但我最后需要这个大约10年前，这就是我所拥有的。 我希望你能从相当广泛的评论中提取相关的操作。 我包括用于起始近似的表的生成以及用于详尽测试函数的测试框架。
#include #include unsigned int tab[96]; __declspec(naked) unsigned int __stdcall fxsqrt (unsigned int x) { __asm { mov edx, [esp + 4] ;// x mov ecx, 31 ;// 31 bsr eax, edx ;// bsr(x) jz $done ;// if (!x) return x, avoid out-of-bounds access push ebx ;// save per calling convention push esi ;// save per calling convention sub ecx, eax ;// leading zeros = lz = 31 - bsr(x) // compute table index and ecx, 0xfffffffe ;// lz & 0xfffffffe shl edx, cl ;// z = x << (lz & 0xfffffffe) mov esi, edx ;// z mov eax, edx ;// z shr edx, 25 ;// z >> 25 // retrieve initial approximation from table mov edx, [tab+4*edx-128];// r = tab[(z >> 25) - 32] // first Newton-Raphson iteration lea ebx, [edx*2+edx] ;// 3 * r mul edx ;// f = (((unsigned __int64)z) * r) >> 32 mov eax, esi ;// z shl ebx, 22 ;// r = (3 * r) << 22 sub ebx, edx ;// r = r - f // second Newton-Raphson iteration mul ebx ;// prod = ((unsigned __int64)r) * z mov eax, edx ;// s = prod >> 32 mul ebx ;// prod = ((unsigned __int64)r) * s mov eax, 0x30000000 ;// 0x30000000 sub eax, edx ;// s = 0x30000000 - (prod >> 32) mul ebx ;// prod = ((unsigned __int64)r) * s mov eax, edx ;// r = prod >> 32 mul esi ;// prod = ((unsigned __int64)r) * z; pop esi ;// restore per calling convention pop ebx ;// restore per calling convention mov eax, [esp + 4] ;// x shl eax, 17 ;// x << 17 // denormalize shr ecx, 1 ;// lz >> 1 shr edx, 3 ;// r = (unsigned)(prod >> 32); r >> 3 shr edx, cl ;// r = (r >> (lz >> 1)) >> 3 // round to nearest; remainder can be negative lea ecx, [edx+edx] ;// 2*r imul ecx, edx ;// 2*r*r sub eax, ecx ;// rem = (x << 17) - (2*r*r)) lea ecx, [edx+edx+1] ;// 2*r+1 cmp ecx, eax ;// ((int)(2*r+1)) < rem)) lea ecx, [edx+1] ;// r++ cmovl edx, ecx ;// if (((int)(2*r+1)) < rem) r++ $done: mov eax, edx ;// result in EAX per calling convention ret 4 ;// pop function argument and return } } int main (void) { unsigned int i, r; // 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 + r; } // exhaustive test of 16.16 fixed-point square root i = 0; do { 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)); break; } i++; } while (i); }
使用您最喜欢的整数平方根算法 ，简单观察√（2 -16 a）= 2-8√a。