Казалось бы, что может быть проще fsqrt, но мы пойдем другим путем 1) Самый известный целочисленный алгоритм для вычисления квадратного корня из числа поражает своей простотой unsigned sqrt_cpu_int(usigned X) { unsigned div = 1, result = 0; while (X > 0) { X -= div; div += 2; result += X < 0 ? 0 : 1; } return result; } Недостаток данного алгоритма - количество итераций будет увеличиваться с ростом Х 2) вычисление квадратного корня по разложению в ряд Тейлора. Пусть Х - любое число; f(X) - некоторая функция, зависящая от X; a - известное число, близкое к Х; f(a) - известное значение функции. Разложим f(X) в ряд Тейлора: f(X)=f(a)+(X-a)f '(a)+((X-a)² *f" (a))/2! + ... + ((X-a)n *f n (a))/n! Пусть X - число, из которого нужно извлечь корень. Тогда f(X)=√[o]X[/o]; a - известное число близкое к X; f(a)=√[o]a[/o] - известное число близкое к √[o]X[/o], и f(X)=√[o]a[/o] +(X-a)/(2√[o]a[/o])+...=(2a+X-a)/(2√[o]a[/o])=(a+X)/(2√[o]a[/o]) Величина √[o]X[/o] может быть найдена, если задаться величиной √[o]a[/o] и затем вычислить f(X). f(X)² можно сравнить с исходным числом Х. Если точность окажется недостаточной, тогда число а заменяется на f(X)², √[o]a[/o] на f(X) и вычисление повторяется 3) поиск целочисленного квадратного корня методом Ньютона начинается с некоторого значения g0, которое является начальной оценкой √[o]X[/o]. Затем выполняется серия уточнений значения квадратного корня по формуле gn+1=(gn+X/gn)/2. Для уменьшения количество итераций можно на первом этапе более точно подобрать начальное значения для переменной g0 { usigned x1; int a, g0, g1; if ( x <= 1 ) return x; a = 1; x1 = x - 1; if (x1 > 65535) {a = a + 8; x1 = x1 >> 16;} if (x1 > 255) {a = a + 4; x1 = x1 >> 8;} if (x1 > 15) {a = a + 2; x1 = x1 >> 4;} if (x1 > 3) {a = a + 1;} g0 = 1 << a; // g0 = 2ª g1 = (g0 + (x >> a)) >> 1; while (g1 < g0) // повторяем, пока приближения строго уменьшаются { g0 = g1; g1 = (g0 + (x/g0)) >> 1} } 4) Так как деление достаточно медленная операция, можно от нее отказаться. Для вычисления квадратного корня используем умножение. 32-разрядное число Х будет иметь 16-разрядный квадратный корень. На каждом шаге происходит уточнение 1 бита значения корня. Для ускорения сделано "безбранчевое" вычисление прибавить или отнять значение в соответствующем разряде Code (Text): mov ebx,4000h mov ecx,8000h mov eax,40000000h @@: cmp eax,X je @f; преждевременный выход из цикла sbb edx,edx; если CF=1 тогда edx=0FFFFFFFFh иначе edx=0 sub ecx,ebx and edx,ebx; если CF=1 тогда edx=ebx иначе edx=0 lea ecx,[ecx+edx*2]; если CF=1 тогда ecx=ecx+ebx иначе ecx=ecx-ebx shr ebx,1; переходим к следующему разряду mov eax,ecx mul eax; получаем "eax в квадрате" test ebx,ebx jnz @b @@: mov eax,ecx; eax:=sqrt(X)
Mikl___ Меня тоже дико восхитил способ 1 ) Думаю есть способы его улучшения. По своей сути этот метод - способ пробежаться по ряду квадратов целых чисел 1, 4, 9, 16, 25, 36... и имхо делать эту пробежку можно более эффективно чем вычитанием.
Y_Mur Это сумма арифметической прогрессии из нечетных чисел, а как вам 4-ый способ? [offtop]рисунок сгенерировался на сайте http://dxdy.ru/math/ по строкам class="latex" alt="$\sum\limits_{n=1}^{k}(2n-1)=k^2$" title="$\sum\limits_{n=1}^{k}(2n-1)=k^2$" /> Вот бы и на WASM.RU/FORUM такое же![/offtop]
Mikl___ Метод 4 хороший... суть - частный метод решения уравнения y=f(x) на отрезке [x2,x1] методом деления отрезка пополам. Пример кода (рекурсивное решение): Code (Text): //x1 > x2 long sqrt_long(long x_2, long x1, long x2) { if(x1 == x2 + 1) return x1; // long x = (x1+x2)>>1;//(x1+x2)/2 // if(x*x >= x_2) return sqrt_long(x_2, x, x2); // return sqrt_long(x_2, x1, x); } ... long x_2 = ... long x = sqrt_long(x_2,65536,0); printf("sqrt(%d)=%d",x_2,x);
Я так понимаю, 3) и 4) (c) Уоррен? Там и ещё один метод есть, без умножений/делений, как и 1), но логарифмической сложности - функции isqrt4/isqrt5 в http://www.hackersdelight.org/HDcode/isqrt.c.txt .
gorodon много итераций. можно попробовать ограничить диапазон учитывая, что старший бит X в X^2 = Y, где значащие биты Y в диапазоне [N, 0], будет на N/2 позиции. те подбирать надо биты только на участке [N/2 - 1, 0] АДД а лучше даже не половинным делением, а тупым побитным перебором с проверкой двигаясь от N/2 - 1 к 0 на 64 бита [63, 0] это будет максимум 30 итераций (а может и нет. не проверял) ЗЫ извинямс, если ктото уже АДД2: Code (Text): unsigned sqrt(unsigned Y){ unsigned N = get_last_bit_pos(Y); unsigned X = 1 << (N >> 1); unsigned bit; if(N == 0) return Y; for(bit = X >> 1; bit > 0; bit >>= 1){ X |= bit; if(X * X > Y) X &= ~bit; // тут можно еще на равенство проверять // if(X * X == Y) return X; } return X; }
3) Да, Уоррен ("Алгоритмические трюки для программистов" -- М. :Издательский дом "Вильямс" 2003 стр 198 Листинг 11.1 Метод Ньютона поиска целочисленного квадратного корня) 4) мой, но не скрою, написал после прочтения Уоррена "Аппаратная реализация поиска целочисленного квадратного корня" и "Бинарный поиск целочисленного квадратного корня"
Если уж речь зашла об Уоррене, то вот метод который меня "убил" то же на ассемблере Code (Text): mov ebp,X mov ecx,40000000h xor eax,eax @@: mov ebx,eax or ebx,ecx shr eax,1 cmp ebp,ebx sbb edx,edx not edx mov edi,edx and edx,ebx and edi,ecx sub ebp,edx or eax,edi shr ecx,2 test ecx,ecx jnz @b ;eax=sqrt(X) Не понятно как, ни деления, ни умножения, но работает! Причем алгоритм был описан аж в 1945 в Von Neumann J. "First Draft of a Reaport on the EDVAC"
Дык это все тот же алгоритм 4 - побитовой обработки... Квадрат от каждого бита записан в m (10^2=100 100^2=10000... (bin)) - потом нужные биты "съезжают" на нужное число разрядов вниз, уже находясь в y. Это как дифференцирование реализовывать с помощью интеграторов... - из той же серии оптимизации...
Чуть исправил код #9 Code (Text): mov ebp,X inc ebp mov ecx,40000000h; здесь для уменьшения количества итераций ;можно подбирать значение, как в методе Ньютона xor eax,eax @@: lea ebx,[eax+ecx] shr eax,1 cmp ebx,ebp sbb edx,edx mov edi,edx and edx,ebx sub ebp,edx and edi,ecx or eax,edi shr ecx,2 jnz @b
Code (Text): mov ebp,X inc ebp bsr ecx,ebp and ecx,-2 mov ebx,1 shl ebx,cl;для уменьшения количества итераций xor eax,eax @@: lea ecx,[eax+ebx] shr eax,1 cmp ecx,ebp sbb edx,edx mov edi,edx and edx,ecx sub ebp,edx and edi,ebx or eax,edi shr ebx,2 jnz @b
