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

SQRTPD вам не поможет?
А если проц новый то вообще вычисли расстояние между жвумя точками дело плёвое
SUBPD -> DPPD -> SQRTSD

Сразу код на асме можно?

Код чего? Как найти расстояние меджу двумя точками в 2-д пространстве если их координаты заданны int64 (а может сразу Float64 хранить?)?

Bishop
> Код чего? Как найти расстояние меджу двумя точками в 2-д пространстве если их
> координаты заданны int64

Нет, если координаты заданы int32 (но корень там придётся брать уже из int64, так как сумма квадратов).
Да, мне нужен асмокод этого с использованием спец.команд. Интересно сравнить как бы.

Ты мне скажи в чём координаты хранишь изначально.

Сори, понял, щас напишу

код написан для FPC из Trunk позднее 4го апреля (тогда применили мой патч для поддержки инструкций SSE3+)

//Предпологаю что A и B не выравняны никак и являються 8-ми байтными областями памяти
//(с возможностью ПРОЧИТАТЬ их как 16-ти байтные области не получив ошибку сегментации)
//Т.к. непомню calling conversion x86 пишу под х64:
//[RCX] - адрес А; [RDX] - адрес B; [RAX] - Result (нижние 32 бита)
Function Dist(Const A, B: Int32_Vec2): Int32; Assembler; NoStackFrame;
Asm
 MOVUPS XMM0, [RCX]
 MOVUPS XMM1, [RDX]
 SUBPD XMM0, XMM1
 DPPD XMM0, XMM0, 000110001b
 SQRTSD XMM0, XMM0
 CVTSD2SI EAX, XMM0
End;   

P.S. Тьфу, совсем плохо уже соображаю, вот верный код

//Предпологаю что A и B не выравняны никак и являються 8-ми байтными областями памяти
//Т.к. непомню calling conversion x86 пишу под х64:
//[RCX] - адрес А; [RDX] - адрес B; [RAX] - Result (нижние 32 бита)
Function Dist(Const A, B: Int32_Vec2): Int32; Assembler; NoStackFrame;
Asm
 CVTDQ2PD XMM0, [RCX]
 CVTDQ2PD XMM1, [RDX]
 SUBPD XMM0, XMM1
 DPPD XMM0, XMM0, 000110001b
 SQRTSD XMM0, XMM0
 CVTSD2SI EAX, XMM0
End;   

Bishop
> //[RCX] - адрес А; [RDX] - адрес B; [RAX] - Result (нижние 32 бита)

Для х86 она такая же, только вместо RAX пишется EAX
А что вообще делает DPPD и что за число 000110001? (равное 49).

Тебе поможет документация Intel по asm командам.
Это инструкция добавлена в SSE4.1 и осуществляет скалярное произведение. imm8 значение 00110001 говорит что нужно взять скал\рное произедение 2-х компонентного вектора (теоретически можно и однокомпонентного) и записать результат в младшие 64-бита SSE регистра (можно в старшие или и туда и туда)

Еще есть интрукция DPPS она может до 4-х компонентных векторов, т.к. работает с Float32, тогда как DPPD с Float64

Увы, у меня эта команда не пошла...

Тогда я решил написать такую функцию:

function Dist(dx, dy: integer): integer;  
var
  a, b: integer;
begin
  a := abs(dx);
  b := abs(dy);
  if a>b then begin
    b := b shr 1 - a*(1-sqrt(3)/2); //???
    if b >= 0 then result := a+b else result := a;
  end else begin
    a := a shr 1 - b*(1-sqrt(3)/2); //???
    if a >= 0 then result := a+b else result := b;
  end;
end;         

В метрике, определяемой этой функцией, круги выглядят, как правильные 12-угольники, функция может быть на 3.5% больше истинного результата, и не может быть меньше. Приближение хорошее.
Но есть проблема.
Вот это вот a*(1-sqrt(3)/2) мне не нравится.
Оптимальный вариант - это умножение на 575416509 (это (1-sqrt(3)/2)*2^32), и взятие старших 32 битов длинного результата (сделать imul и использовать edx из пары edx:eax, в которую пишется результат умножения).
К сожалению, в Д7 нету никакого оператора для этого простого действия.
Пришлось делать асмотрах.
Как без еретичного асма сделать это? Вариант (a*8780) shr 16 годится только для не очень больших чисел, а у меня большие бывают.

TarasB
если устраивает такая точность можно сделать таблицей из 64 или 256 (2^(2n)) значений

Aslan
> если устраивает такая точность можно сделать таблицей из 64 или 256 (2^(2n))
> значений

И как я к таблице буду обращаться? Таблица-то нужна с шагом, пропорциональным значению.
То есть не 2,4,6,8,10,12,14, ..., а 1,2,5,10,20,50,100...

Кому интересно, вот рабочий асмотрах-вариант:

function Dist(dx, dy: integer): integer;   
asm
  push ebx
  mov ecx, edx    // ecx = dy
  cdq
  xor eax, edx  
  mov ebx, eax
  sub ebx, edx    // ebx = abs(dx)
  mov eax, ecx    // ebx = abs(dx), eax = dy
  cdq
  xor eax, edx
  sub eax, edx    // ebx = abs(dx), eax = abs(dy)
  cmp eax, ebx
  jl @1
  // eax>ebx
  mov ecx, eax    // ecx = abs(dy)
  mov edx, 575416509
  imul edx         // edx = abs(dy) * (1-sqrt(3)/2)
  shr ebx, 1
  mov eax, ecx
  sub ebx, edx
  jl @2
  add eax, ebx
  @2:
  pop ebx
  ret
  @1:
  // ebx>eax
  mov ecx, eax
  mov eax, ebx    // eax,ebx = abs(dx), ecx = abs(dy)
  mov edx, 575416509
  imul edx         // edx = abs(dx) * (1-sqrt(3)/2)
  shr ecx, 1
  mov eax, ebx
  sub ecx, edx    // ecx = abs(dy)/2 - abs(dx)*(1-sqrt(3)/2)
  jl @3
  add eax, ecx
  @3:
  pop ebx
end;

Bishop
> Как найти расстояние меджу двумя точками в 2-д пространстве
Если тебя устраивает 3.5% погрешность, то устроят и Кармаковские 6% (1.5*100/sqrt(2)-100):

int P_AproxDistance(int dx, int dy)
{
  dx = abs(dx);
  dy = abs(dy);
  if (dx < dy)
    return dx+dy-(dx>>1);
  return dx+dy-(dy>>1);
}

entryway
Этот вариант был в 27 посте (я его со времён 27 поста как раз до такого вида и дооптимизировал)...
Только нафига dx+dy-(dy>>1) вместо dx+dy>>1?
Выдаёт 8-угольник, для него расхождение не 6, а 12 процентов. Нет, он мне не годится.
А 12-угольник на 25% всего медленнее.

TarasB
int64 разве нет в твоем делфи чтобы сделать (a*8780) shr 16? Ну и разрядов еще парочку сократить можно.

Такое еще загуглил первой ссылкой:
http://www.flipcode.com/archives/Fast_Approximate_Distance_Functions.shtml