Фантазии о Вселенной и мой личный сайт
Страница запроса

Вычисление квадратного корня

model tiny
codeseg

    mov     bp, 47714       ; задаем входные данные
    xor     bx, bx          ; нижний предел
    mov     cx, 0xFFFF      ; верхний предел
rt: push    bx
    inc     bx
    cmp     bx, cx
    pop     bx
    je      ex              ; если BX = CX - 1, то выход
    mov     ax, bx
    add     ax, cx          ; AX = (BX + CX)
    pushf                   ; сохраняется CF
    shr     ax, 1           ; AX = AX div 2
    popf
    jnc     sk              ; если CF = 0, ничего не менять
    or      ah, 0x80        ; иначе записать в старший бит 1
sk: push    ax
    xor     dx, dx
    mul     ax              ; DX:AX = AX*AX
    cmp     dx, bp          ; DX < BP?
    pop     ax
    jnc     bg              ; если нет, то искомое число должно быть меньше
    mov     bx, ax          ; иначе поднимаем нижний предел
    jmp     rt
bg: mov     cx, ax          ; если искомое число меньше
    jmp     rt              ; то опускаем верхний предел
ex: ; AH = целое, AL = дробное
    ; CALL ... печать числа ...
    int     0x20            ; В десятичной системе дробь = AL / 256
Программа выглядит страшно и ужасно, хотя, для кого как... В действительности, в ее основе лежит очень легкий алгоритм двоичного поиска. Скорость выполнения этой программы довольно низкая, учитывая то, что здесь используется log2(65536) = 16 проходов для получения окончательного результата. Сама же программа занимает 48 байт. Честно сказать, все это дело можно было решить тремя инструкциями FPU (Float-Point Unit):
FLD [Число]
FSQRT
FSTP [Его_Корень]
с гораздо большей производительностью, но... алгоритм бинарного поиска можно использовать не только в таких ситуациях, как вычислить квадратный корень или, например, найти результат деления двух ГИГАНТСКИХ чисел, по 7000 знаков каждый, но еще и при организации поиска по упорядоченному списку, например. Так что некоторая польза от этого алгоритма точно будет.

Первые 3 строчки - инициализация, где BP - число, из которого нужно вычислить корень, а BX - нижний, а СX - верхний предел. Поиск ведется половинчатыми делениями между пределами. То есть, находится среднее арифметическое: AX = (BX + CX) div 2
0          31767       65535
BX [ ------ AX ------ ] CX
Умножая AX*AX / 65536, получаем некое число... Если оно меньше, чем BP, значит, AX нужно сместить в сторону CX, чтобы увеличить получаемое число. Если же оно больше, чем BP, то AX смещаем в сторону BX, то есть в сторону уменьшения AX. Для этого нужно установить новые пределы.
  • Если верхний предел CX будет равен AX, то следующее значение AX будет лежать в пределах от BX до текущего AX.
  • Так и наоборот, если BX = AX, то следущий AX будет лежать между AX текущим и верхним пределом.
Постепенно BX и CX будут сближаться друг с другом, пока BX не будет на единицу меньше, чем СX. Здесь же BX < AX < CX, а, по сути, почти равно по теореме о двух сходящихся функциях.

С технической стороны объяснения требуют два момента - при нахождении среднего арифметического и при умножении числа.

Найти среднее арифметическое BX и CX просто:
MOV AX, BX		; AX = BX
ADD AX, CX		; AX = AX + CX = BX + CX
SHR AX, 1		; сдвиг AX на единицу вправо, т.е. фактически,
			; деление надвое
Сложность заключается в том, что при сложении 2 чисел, сумма которых превышает 65535, получится число, меньшее 65536. Старший бит уйдет в так называемых "флаг переноса", Carry Flag (CF) и будет установлен в 1.

Например, складывая 0x895A + 0xABCD = 0x3527 а старший бит CF = 1. Итак, запомнив состояние флагов с помощью инструкции PUSHF, делим полученный результат пополам. Далее, восстанавливаем ранее сохраненные флаги инструкцией POPF, потому что те могут измениться в процессе деления, смотрим, установлен ли CF в единицу. Если да, то это говорит о том, что он должен оказаться в самом старшем бите полученного результата от деления... Почему? Перейдем в двоичную систему счисления.
	   10001001 01011010	  0x895A
	+  10101011 11001101	+ 0xABCD
	--------------------	  ------
	1  00110101 00100111    0x1 3527
	CF    AH       AL

	1  x0011010 10010011	  0x1A93
	\--^
SHR AX, 1 сместил ровно на один бит вправо весь регистр AX, таким образом, поделив его надвое, без учета того бита, что был в CF. Так что если CF = 1, старший бит, помеченный как "x" будет установлен тоже в 1. Если CF = 0, то "x" = 0 тоже.

С этим все понятно... Надеюсь. Теперь с числами. В программе имеет место быть вот что: AH - это целая часть получаемого числа, AL - дробная. Если перевести в нормальный язык математики, то a = AX / 256. Чтобы представить, допустим, a = 4,564 в регистре AX, нужно a*256 = AX, что совершенно очевидно. Здесь AX будет 4,564 * 256 = 1168. При возведении в степень AX*AX = (a*256)*(a*256) = a*a*65536. Целая часть квадрата [a] окажется за пределами регистра AX, то есть в регистре DX. Есть функция MUL AX, AX, она возводит в квадрат число AX и помещает в регистр DX верхние 16 бит, которые получились при умножении, а в AX - нижние 16, которые, в данном случае, будут дробью.

Далее число DX сравнивается с BP (в котором записано целочисленное значение числа, из которого нужно вычесть квадратный корень). Если DX меньше, то AX смещается в сторону верхнего предела, чтобы DX стало больше, если DX больше - то AX уходит в сторону, соответственно, нижнего предела, ближе к BX.

Таким образом, итерируясь, AX приходит к нужному значению. После чего вызывается подрограмма печати этого числа... Таков алгоритм выполнения бинарного поиска в этой программе.

***
Существует и другой подход к вычислению квадратного корня. Его можно использовать и на практике, он тоже вылеплен примерно "из того же теста". Это итерационный математический алгоритм:

a = 1/2 * (a + b/a)

где [b] = число, из которого необходимо вычислить квадратный корень, и [a] = результат. При запуске алгоритма желательно, чтобы [a] = [b]. Докажу правомерность этого выражения:
a = 1/2 * (a + b/a)
2*a = a + b/a;
2*a - a = b/a;
a = b/a;
a*a = b;
a = b1/2 = sqrt(b);
Суть в том, чтобы a и b/a стали оба равны [a], а так как они оба будут равны [a], то их сумма получится равной 2*a, потому эта сумма делится надвое, получая [a]. Мудреж... Все по полочкам. Допустим, увеличивается число A, то, в таком случае, уменьшается значение B/A т.к. они обратно пропорциональны. Между двумя этими значениями A и B/A ищется среднее арифметическое по формуле (A + B/A) * 0.5, получая, в конечном итоге, новую A, лежающую между двумя своими предыдущими значениями. На втором шаге снова повторяется тоже самое, но только в этом случае расстояние между A и B/A сокращается... на следующем шаге еще больше сокращается. Значит, если lim |A - B/A| = 0, что говорит о том, что B/A = A. И отсюда B = A*A. Выводы можно сделать очевидные...

Этот чудесный алгоритм (который выдумал Ньютон, а я перевыдумал), можно расширить. Например, если записать A = 1/2*(A + B/[A*A]), то при нужных преобразованиях получаем A = B / [A*A] = A*A*A = B, т.е. кубический корень. Обобщим: A = 1/2*(A + B/An-1) теперь можно получать корень n-й степени. Чудесно... Хотя калькулятор взять намного проще и нажать пару кнопок :-)