Java中s15.16定點數的平方根

Question

我想編寫一個函數來計算s15.16定點數的平方根。我知道它是一個有15位數字和16位數字的有符號數字。無論如何沒有任何圖書館就可以做到這一點？其他任何語言都可以。

Answer 1

使用your favourite integer square root algorithm，用簡單的觀察，即√（2 ^-16一個）= 2 ^-8√A。

Answer 2

我假設你問這個問題是因爲你所在的平臺沒有提供浮點數，否則你可以通過浮點平方根實現15.16個定點平方根，如下所示（這是C代碼，I假設Java代碼看起來非常相似）：

int x, r; 
r = (int)(sqrt (x/65536.0) * 65536.0 + 0.5); 

如果你的目標平臺提供了一種快速整數乘法（尤其是，無論是雙寬度結果或乘法高位指令乘），你可以騰出一些一個小表的內存，使用牛頓 - 拉夫遜迭代加上一個基於表格的起始逼近通常是要走的路。通常，近似於倒數平方根，因爲它具有更方便的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); 
}