Квадратный корень s15.16 номер фиксированной точки в Java

Question

Я хочу написать функцию для вычисления квадратного корня из числа неподвижных точек s15.16. Я знаю его подписанное число с 15-значным int и 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); 

Если ваша целевая платформа обеспечивает быстрое умножение целого чисел (в частности, либо многократно с двойной шириной результата или многосвязным высокой инструкцией), и вы можете сэкономить некоторые память для небольшой таблицы, использование итераций 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); 
}