java - Java中s15.16定点数的平方根

Question

我想编写一个函数来计算 s15.16 定点数的平方根。我知道它是一个有 15 位整数和 16 位小数的有符号数。无论如何，没有任何图书馆就可以做到吗？任何其他语言也可以。

Answer 1

我假设您问这个问题是因为您所在的平台不提供浮点，否则您可以通过浮点平方根实现 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 相乘得到平方根，即 sqrt(x) = rsqrt(x) * x。以下代码显示了如何以这种方式计算正确舍入的 16.16 定点平方根（因为平方根的参数必须是正数，这对于 s15.16 定点同样适用）。通过最小化残差 x - sqrt(x)*sqrt(x) 来执行舍入。

很抱歉，平方根函数本身是 32 位 x86 内联汇编代码，但我上次需要它大约是在 10 年前，这就是我所拥有的。希望大家能从相当广泛的评论中提炼出相关的操作。我包括了用于起始近似值的表格的生成以及详尽地测试函数的测试框架。

#include <stdlib.h>
#include <math.h>

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);
}

Answer 2

使用你最喜欢的整数平方根算法，简单的观察是 √(2 ^-16 a) = 2 ^-8 √a。