s15.16 固定小数点数の平方根を計算する関数を書きたいと思います。15 桁の整数と 16 桁の小数を持つ符号付きの数値であることはわかっています。ライブラリなしでそれを行う方法はありますか?他の言語でも構いません。
2 に答える
あなたが使用しているプラットフォームが浮動小数点を提供していないため、この質問をしていると思います。それ以外の場合は、次のように浮動小数点平方根を介して 15.16 固定小数点平方根を実装できます (これは C コードです。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);
}
お気に入りの整数平方根アルゴリズムを使用して、√(2 -16 a) = 2 -8 √aという単純な観察を行います。