- 0iStalker
- Модератор
Вариант для 64бит, -
unsigned __int64 isqrt64(unsigned __int64 x) { unsigned __int64 x1; __int64 s, g0, g1; if ( x<=1) return x; s = 1; x1 = x - 1; if ( x1>4294967295) {s = s + 16; x1 = x1 >> 32;} 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
так ты с хмм-вариантом сравни
искрт4 тормозной же, как выяснилось
И этот хмм умеет работать с 64 битами (типа брать число из регистровой пары еах:едх)?
TarasB
>ну это и есть вариант с сопроцессором, я так тоже делал
да я понял
просто так удобно и быстро
но если копать в сторону bad_sqrt, то наверное тоже самое можно сделать и для int64
кстати потесть скорость (по сравнению с предыдущим вариантом)
function bad_fast_sqrt(x: cardinal): cardinal; asm bsr ecx, eax jnz @1 ret @1: test ecx, 1 jz @even @odd: shr ecx, 1 mov edx, 1 shr eax, cl shl edx, cl shr eax, 1 mov ecx, $B504F334 add eax, edx mul ecx mov eax, edx ret @even: shr ecx, 1 mov edx, 1 shr eax, cl shl edx, cl shr eax, 1 shr edx, 1 add eax, edx end;
Ты переставил инструкции, чтобы лучше параллелилось? Ну, ускорение на 4 процента (время на 32 млн вычислений уменьшилось с 625 до 609).
0iStalker
твой вариант медленно сходится, надо перед методом Ньютона получить первые 2-3 бита, что делается подстановкой 10(2), 11(2), 100(2), 101(2), 110(2), 111(2)
автор хочет вообще без делений
TarasB
твой первый алгоритм не понял как работает, можешь объяснить принцип?
TarasB
> так ты с хмм-вариантом сравни
Без учёта времени на получение случайного 32х битного числа - xmm быстрее в 4.3 (124 vs 541мс) раза... Но это наверняка еще от процессора зависит.
Aslan
> твой первый алгоритм не понял как работает, можешь объяснить принцип?
Нибубу, алгоритм не мой, это алгоритм какого-то чувака времён 45 года.
Я только после долгой отладки и подстановки разных чисел начал так смутно понимать, что он делает, но объяснить это смутное понимание не могу.
0iStalker
> xmm быстрее
Ну вот.
Осталось понять, как ему скормить регистровую пару.
TarasB
странно
должно было процентов на 30
а так?
какая разница с stdrt:
function IntSqrt(const X: cardinal): cardinal; asm push eax fild dword ptr [esp] fsqrt fistp dword ptr [esp] pop eax end;
DevilDevil
> а так?
> какая разница с stdrt:
дык это и есть мой fsqrt.
Только вместо push eax я писал mov [esp-4], eax
Или так хуже?
На всякий случай замечу, что соотношение производительности CPU и FPU у Atom cущественно отличается от такового для "обычных" процессоров Intel и AMD, причем у последних это отношение примерно одинаково.
Так что если речь идет об оптимизации производительности, то результат обязательно нужно проверять на "десктопных" процессорах.
andriano
> На всякий случай замечу, что соотношение производительности CPU и FPU у Atom
> cущественно отличается от такового для "обычных" процессоров Intel и AMD
В какую сторону?
andriano
> Так что если речь идет об оптимизации производительности, то результат
> обязательно нужно проверять на "десктопных" процессорах.
Селерон 600 МГц покатит? На нём аналогичный результат, только цифры побольше.
TarasB
> Селерон 600 МГц
Да вы, батенька, счастливый обладатель раритета ! :) Софтовый рендер пишете?
Насколько я понял, в том алгоритме идет последовательный подбор битов, начиная со старшего, т.е. если число 32-битное, то начинаем с 2^15, берем квадрат, если больше аргумента, то записываем этот 15й бит 1,
и так далее, только взятик квадрата оптимизированно, так что без умножений, см. мой вариант
unsigned int isqrt(unsigned int a) { // 0<=x<2^32 => 0<=y=sqrt(x)<2^16 unsigned int m=( 1<<30); // 2^(2k) unsigned int n=( 1<<15); // 2^k unsigned int x=0; // результат unsigned int y=0; // x^2 if( a>=m) { x=n; y=m; } for( int k=14;k>=0;--k) { m>>=2; // если (x+n)^2<=a, x=x+n unsigned int y2=y+( x<<k)+m; if( y2<=a) { x+=n; y=y2; } n>>=1; } return x; }
Вообще, если не ошибаюсь, на современных CPU умножение,деление одинаковы по скорости со сложением, даже floating point
Я тут немного подумал, и решил, что мне нужен не сам корень, а расстояние между двумя точками. То есть функция sqrt(sqr(a)+sqr(b))
Я для неё взял такой аналог:
(abs(a-b)+(abs(a)+abs(b))*3) div 4
это полусумма ||a||1 и ||a||inf, то есть (max(a,b) + abs(a)+abs(b)) div 2, то есть ((a+b+abs(a-b)) div 2 + abs(a)+abs(b)) div 2 то есть то, что я сказал.
Короче, в 10 раз быстрее, чем по формуле Пифагора. Отклонение от истины 0.89 (2/sqrt(5)).
Что-то мне не понятен ваш энтузиазм. Придумали SSE для вычисления квадратного корня. Чем SSE не угодило?
SQRTPS xmm, xmm/m команда вычисляет квадратный корень для каждого из четырех чисел во входном операнде и записывает результаты в выходной операнд.
SQRTSS xmm, xmm/m команда вычисляет квадратный корень из младшего элемента входного операнда и записывает результат в младший элемент в выходной операнд. Остальные элементы выходного операнда не меняются.
RCPPS xmm, xmm/m команда определяет приближенное обратное значение для каждого из четырех чисел входного операнда и записывает результаты в XMM-регистр.
RSQRTPS xmm, xmm/m команда вычисляет приближенное обратное значение для квадратного корня из каждого из четырех чисел входного операнда и записывает результаты в XMM-регистр.
RCPSS xmm, xmm/m команда определяет приближенное обратное значение для числа, находящегося в младшем элементе входного операнда и записывает результат в младший элемент выходного операнда. Остальные элементы выходного операнда не меняются.
RSQRTSS xmm, xmm/m команда вычисляет приближенное обратное значение для квадратного корня из числа , находящегося в младшем элементе входного операнда и записывает результат в младший элемент выходного операнда. Остальные элементы выходного операнда не меняются.
Доступны начиная с Pentium 3.
KKH
> Что-то мне не понятен ваш энтузиазм. Придумали SSE для вычисления квадратного
> корня. Чем SSE не угодило?
Как подружить с 64 битами?
Вообще нету ли готовой команды
x -> y -> sqrt(sqr(x)+sqr(y))?
Тема в архиве.