Быстрое [s]взятие квадратного корня из целого 64-битного числа[/s] нахождение расстояния между двумя целочисленными точками.

Для начала я решил потренироваться на 32-битных числах.

Вариант-1, найден тут:
http://www.cyberforum.ru/assembler/thread214648.html

int isqrt4(unsigned x) {   unsigned m, y, b;
   m = 0x40000000;
   y = 0;
   while(m != 0) {              // Do 16 times.
      b = y | m;
      y = y >> 1;
      if (x >= b) {  x = x - b; y = y | m; }
      m = m >> 2;
   }
   return y;
}

Там же была ассемблерная копия. Я правда, не знаю, что там автор в своём асмокоде намутил, потому что он работает не лучше выхлопа убогого оптимизатора Д7. Я намутил свой асмокод.

function msqrt(x: cardinal): integer;
asm
  push ebx     
  mov ecx, eax
  bsr ecx, ecx
  and ecx, -2
  mov edx, 1
  shl edx, cl
  mov ecx, eax
  xor eax, eax

  @loop:
    lea ebx, [eax+edx]
    shr eax, 1

    cmp ecx, ebx
    jl @
      sub ecx, ebx
      add eax, edx
    @:

    shr edx, 2
  jnz @loop;

  pop ebx
end;

Оказалось раза в полтора быстрее, чем тот вариант.

Потом я решил тупо отправить число сопроцессору. Оказалось ещё в полтора раза быстрее. К тому же этот метод позволяет считать квадратный корень не только из 32-, но и из 64-битных целых.

Ещё вдвое сделал всех вариант с ХММ (тоже там по ссылке). Правда, я не знаю, умеет ли он извлекать корень из 64-битных чисел.

  cvtsi2ss xmm0,eax
  sqrtss   xmm0,xmm0
  cvtss2si eax,xmm0

И, наконец, ещё один мой вариант. Он грубый, но в полтора раза быстрее последнего.
Он основан на том, что sqrt(x) = sqrt(2^(2n)+y) ~ 2^n + y/(2^(n+1)) = x/(2^(n+1)) + 2^(n-1)
это начало ряда тейлора
Эта формула работает только, если номер первого бита чётный. Для нечётного надо иначе.

sqrt(x) = sqrt(2^(2n-1)+y) = sqrt(2^(2n) + 2y) / sqrt(2) ~ (2^n+y/(2^n)) / sqrt(2) = (x/(2^n)+2^(n-1)) * (1/sqrt(2))

асмокод такой:
(на самом деле асм только ради команды bsr)

function badsqrt(x: cardinal): integer;
asm
  bsr edx, eax
  jz @exit

  mov ecx, edx
  test ecx, 1

  jz @even
  @odd:

    shr ecx, 1
    inc ecx
    shr eax, cl
    mov edx, 1
    shl edx, cl
    shr edx, 1
    add eax, edx

    mov edx, $B504F334    //  это 2^32 / sqrt(2)
    mul edx
    mov eax, edx

    ret
  @even:

    shr ecx, 1
    inc ecx
    shr eax, cl
    mov edx, 1
    shl edx, cl
    shr edx, 2
    add eax, edx

  @exit:
end;

Само умножение незначительно влияет на скорость.

Результат такого корня от аргумента, большего 100, отклоняется от истинного значения не более, чем на 10%.

А вот результаты скорости вычисления. Я 16 миллионов раз брал корень из 25 000 000. Проц - атом 1.6 ггц
В правой колонке - время в миллисекундах

pascal: 2016
asm:    1938
my_asm: 1297
fsqrt:  1062
sqrtss: 531
bad:    312
stdrt:  1250

pascal - это алгоритм по ссылке, записанный на Паскале
asm - это асмокод по ссылке
my_asm - это тот же асмокод, но переписанный мной
fsqrt - это через сопроцессор
sqrtss - это через XMM
bad - это грубый sqrt
stdrt - это лобовое решение, trunc(sqrt(x))

Так вот, у меня вопрос - как быстрее всего брать квадратный корень из целого 64-битного? Надо для вычисления расстояний в играх, определения столкновений и всякой такой фигни, точность, короче, не важна.

http://algolist.manual.ru/maths/count_fast/sqrt.php

Был я там. Забраковал сразу из-за деления.

быстрее всего не брать квадратный корень.

x
> быстрее всего не брать квадратный корень.

Я знаю. Для сравнения расстояний, например, я его не беру.
А вот для определения глубины проникновения одного круга в другой уже не обойтись без корня.

а для чего эта глубина?
и почему вычисления именно в целых числах, пусть и 64 битных?

TarasB
> Был я там. Забраковал сразу из-за деления.

Копиманки? А теорию почитать? А мозг включить?

unsigned isqrt(unsigned x)
{
    unsigned x1;
    int s, g0, g1;

    if (x<=1) return x;

    s = 1;

    x1 = x - 1;
    if (x1>65535) {s = s + 8; x1 = x1 >> 16;}
    if (x1>255) {s = s + 4; x1 = x1>>8;}
    if (x1>15) {s = s + 2; x1 = x1>>4;}
    if (x1>3) {s = s + 1;}

    g0 = 1<<s;

   g1 = (g0 + (x>>s)) >> 1;
   while (g1 < g0)
   {
       g0 = g1;
       g1 = (g0 + (x/g0))>>1;
   }

   return g0
}

0iStalker
> (x/g0)
0iStalker
> делений 0

?

TarasB

5 штук в худшем случае, для 32х битного числа,... причем одно из них заменено сдвигом.

0iStalker
Хорошо, сравни по скорости с решением на хмм. У меня нет хорошего С-компилятора, поэтому я не знаю.

TarasB
>sqrtss xmm0,xmm0
>cvtss2si eax,xmm0
>ss
>eax
Очевидно, double'ам здесь нет пути.

я бы наверно сделал так

function Int64Sqrt(const X: int64): cardinal; 
asm
  fild X
  fsqrt
  fistp dword ptr [esp-4]
  mov eax, [esp-4]
end;


// чуть быстрее. без префиксов/постфиксов и X передаётся через регистр, а не стек
function Int64SqrtP(const X: pint64): cardinal; 
asm
  fild qword ptr [eax]
  fsqrt
  fistp dword ptr [esp-4]
  mov eax, [esp-4]
end;

TarasB
> Хорошо, сравни по скорости с решением на хмм. У меня нет хорошего
> С-компилятора, поэтому я не знаю.

Запустил в цикле на 5млн. итераций

isqrt4 - 0.000116 сек на итерацию (на все 562мсек)
isqrt - 0.000072 сек на итерацию (на все 360мсек)
генерация случайного числа по формуле n = rand()*rand() -  0.0000376сек (на все 188мсек)

Запустил в цикле на 16млн. итераций

isqrt4 - 0.000115 сек на итерацию (на все 1844)
isqrt - 0.000071 сек на итерацию (на все 1141мсек)
генерация случайного числа по формуле n = rand()*rand() -  0.0000361сек (на все 578мсек)

DevilDevil
ну это и есть вариант с сопроцессором, я так тоже делал

0iStalker
isqrt4 - это что?
isqrt - это вариант, что ты привёл?

TarasB
> isqrt4 - это что?

то, что бы в #0 сообщении

>isqrt - это вариант, что ты привёл?

Да