Квадратный корень.
(c) Alone Coder
На эту тему я потратил лучшую неделю :)
Всё-таки расколоть безбашенный алгоритм
нахождения целой части квадратного корня,
где что-то сдвигается-вычитается-сдвигает-
ся-вычитается... и почему-то получается
правильный результат :)
Надеюсь, все видели компилятор MCODER2
или хотя бы что-то другое, где без коммен-
тариев приведён этот SQR? (А то вдруг я
так радуюсь зря, а читатели смотрят на эти
большие буквы большими глазами, прикиды-
вая, чего это я съел, что такой счастли-
вый? ;)
Вот процедура:
;HL=аргумент N
OR A
LD A,L
LD L,H
LD DE,64
LD H,D
LD B,8
sqr0 SBC HL,DE
JR NC,$+3
ADD HL,DE
CCF
RL D
ADD A,A
ADC HL,HL
ADD A,A
ADC HL,HL
DJNZ sqr0
;D=результат q=[√N].
;конец фильма. спасибо за внимание ;)
Вы поймёте процедуру в этом виде? Пра-
вильно,не поймёте.И я не пойму.И тем более
будете недоверчиво относиться к идее её
использования, особенно не проверив прави-
льность результат на всех числах от 0 до
32767. А увеличивать разрядность результа-
та вы тем более не возьмётесь, не зная,как
программа будет работать после этого.
Первый шаг.
Всем уважающим себя программистам изве-
стно,что квадрат натурального числа q мож-
но представить в виде суммы последователь-
ных нечётных чисел, от единицы до 2q-1.
Только извлекать через это квадратный
корень (методом вычитания возрастающих не-
чётных, пока они не перестанут вычитаться)
слишком, пожалуй, расточительно.
Читаем книжку.
Первый намёк я увидел здесь:
(Ж.Арсак "Программирование игр и головоло-
мок", дословно)
"Можно обобщить предыдущий алгоритм,ис-
пользуя свойства десятичной записи чисел.
Данное число разделяется на куски по две
цифры, начиная справа; затем мы начинаем
вычитать последовательные нечётные числа
из крайнего слева куска:
1909 1809 1509 1009
- 1 - 3 - 5 - 7
---- ---- ---- ----
1809 1509 1009 309
Если это нельзя продолжать дальше, то
последнее вычитаемое число увеличивается
на единицу,сдвигается на один шаг вправо,и
следом за ней приписывается единица. Это -
первое нечётное число, которое следует вы-
читать из предыдущего остатка.
В приведённом выше примере 7+1=8; при-
писывая 1, получаем 81 и продолжаем:
309 228 145
- 81 - 83 - 85
---- ---- ----
228 145 60
Поскольку продолжать дальше нельзя (по-
следнее возможное вычитание из остатка -
это крайнее справа), то последнее из вычи-
таемых чисел нужно увеличить на 1, а затем
разделить на 2, чтобы получить корень. По-
следний остаток и есть остаток квадратного
корня:
85+1=86, 86/2=43,
1909=(43)¤+60.
Этот алгоритм достаточно прост для про-
граммирования при длинных числах,и он даёт
вполне разумное время вычисления.
У вас много возможностей представлять
свои данные.Так как мы оперируем с кусками
из двух цифр, то вы можете задавать свои
данные таблицами целых чисел в интервале
от 0 до 99.
Вы можете представлять свои целые числа
как цепочки символов, где используются
только численные символы (цифры) от 0 до
9. Выбор способа зависит от ваших предпоч-
тений и от возможностей вашей машины опе-
рировать с таблицами и цепочками. Тщатель-
но рассмотрите, какие операции нужно сде-
лать. Вы ничем не ограничены: почему бы
не запрограммировать и сравнить два разных
решения?
Я предложил вам алгоритм без доказате-
льства. Поэтому постарайтесь его прове-
рить...
Я предложил вам алгоритм для десятичной
системы счисления. Можно предложить похо-
жий алгоритм для двоичной системы. Тогда
не возникает цикл вычитаний последователь-
ных нечётных чисел из каждого куска, пос-
кольку в куске только одно нечётное число:
1. Алгоритм упрощается: если можно вычесть
нечётное число - мы его вычитаем, в проти-
вном случае не делаем ничего. Затем сдви-
гаем, добавляем 1 и приписываем 1 в кон-
це... Этот алгоритм намного проще реализо-
вать. Но тогда нужно сначала перейти к ос-
нованию 2, а затем преобразовать двоичный
результат в десятичный. Вам следует посмо-
треть, что более эффективно..."
Несмотря на неудачный перевод и отсутс-
твие объяснений, "почему", то бишь, откуда
следует правильность алгоритма,всё же оче-
видно, что алгоритм - наш. Имеется в виду,
тот же,что реализован выше. То есть,теперь
мы получили его словесное описание в пере-
ложении на систему счисления с произволь-
ным основанием.
Но нам потребуется именно двоичный слу-
чай. (Невыгодность основания системы счис-
ления, большего 2, следует хотя бы из того
факта, что при бесконечном основании зада-
ча извлечения корня сводится опять-таки к
последовательному вычитанию нечётных чи-
сел.)
Объяснение.
В чём же суть метода?
Совершенно понятно,что первую слева ци-
фру результата мы можем получить всегда и
легко, используя всего лишь две левые циф-
ры исходного числа.
Метод состоит в том,что мы получаем ко-
рень "из первых четырёх цифр", используя
уже полученный корень "из первых двух
цифр" и так далее.
Перейдя сразу к двоичному случаю, мы
можем упростить понимание.
Представьте себе числовую ось, на кото-
рой помечены отрезки нечётных длин, увели-
чивающиеся слева направо, из каковых отре-
зков и "состоит" наше исходное число.Пусть
также число N не является точным квадра-
том, тогда оно окажется на оси правее ква-
драта корня q=[√N], отличаясь от него на
некоторый остаток квадратного корня.
Можно также обратить внимание на тот
факт, что последнее прибавленное нечётное
равно 2q-1, а следующее было бы 2q+1, но N
оказалось слишком мало для этого. Поэтому
очевидно, что остаток квадратного корня не
превышает 2q.
Пусть нам стало известно,что предыдущее
приближение корня (то есть целая часть ко-
рня из N/4 ) равно r, а нам нужно получить
ещё один бит результата, то есть выяснить,
чему равно q: 2r или же 2r+1.
Мы можем это узнать всего лишь одним
сравнением. Если N меньше, чем 4r¤+4r+1
(где 4r+1 - это "следующее прибавляемое
нечётное" к квадрату числа 2r, оно же "по-
следнее прибавляемое нечётное" к квадрату
числа 2r+1), то q=2r, иначе q=2r+1. Логич-
но?
Доказательство.
Особо головастых товарищей вполне уст-
роит нижеприведённое доказательство :)
Дано N Є 0
_
Найти q = [√N]
0 ___
┌ / N ┐
Определим q как │√ ─── │
i └ 4^i ┘
(квадратные скобки - взятие целой части,
знак ^ - возведение в степень)
Очевидно, такое q должно обладать свойст-
вами: i
N
(ё) q ¤ є ───
i 4^i
N
(ж) (q +1)¤ > ───
i 4^i
Обозначим через imax число пар бит в двои-
чном представлении N.
imax > log N
4
Опять-таки, очевидно, что
_____
┌ / N ┐
q = │√ ───── │ = 0
imax └ 4^imax┘
Начинаем проводить доказательство индукци-
ей по i. Пусть известно q (q известно)
i imax
Докажем, что q , получаемое по следующей
i-1
формуле, удовлетворяет свойствам (ё) и (ж)
(тогда q и будет решением).
0
i i-1
┌ 2q , если N-4 q ¤< (2(2q )+1)·4
│ i i i
q = ┤
i-1 │ 2q +1 в противном случае (...Є...)
└ i
Давайте запишем эти условия ещё раз в
следующем виде:
верхнее:
i 1 i
(ь) если N-4 q ¤< (q + ─)·4
i i 4
нижнее:
i 1 i
(ъ) если N-4 q ¤ Є (q + ─)·4
i i 4
Докажем верхний случай.
N
1. Докажем (ё) для q : q ¤є ───────
i-1 i-1 4^(i-1)
что, пользуясь q = 2q , можно записать в
виде: i-1 i
N N
4q ¤ є ───────, т.е. q ¤ є ───,
i 4^(i-1) i 4^i
что и известно по свойству (ё) для q .
i
2. Докажем (ж) для q , т.е.
i-1
N
(q +1)¤ > ───────, т.е.
i-1 4^(i-1)
N
q ¤+2q +1 > ───────, т.е. (пользуясь
i-1 i-1 4^(i-1)
выражением q = 2q ):
i-1 i
N
4q ¤+4q +1 > ───────. Умножим на 4^(i-1):
i i 4^(i-1)
i i i-1
4 q ¤+4 q +4 > N, т.е.:
i i
i i i-1
N-4 q ¤ < q ·4 +4
i i
что известно по условию (ь).
Теперь докажем нижний случай.
N
3. Докажем (ё) для q : q ¤є ───────
i-1 i-1 4^(i-1)
что, пользуясь q = 2q +1, можно записать
в виде: i-1 i
N
4q ¤+4q +1 є ───────. Умножим на 4^(i-1):
i i 4^(i-1)
i i i-1
N-4 q ¤ Є q ·4 +4
i i
что известно по условию (ъ).
4. Докажем (ж) для q , т.е.
i-1
N
(q +1)¤ > ───────, т.е. (пользуясь выра-
i-1 4^(i-1)
жением q = 2q +1 ):
i-1 i
N N
(2q +2)¤ > ───────, т.е. (q +1)¤ > ───,
i 4^(i-1) i 4^i
что опять-таки известно по свойству (ж)
для q .
i
На этом доказательство окончено. (Возг-
ласы: "Что-то ты не то доказал, дядя!! Где
НАШ-ТО алгоритм???" ;))
Начинаем упрощать!
i 1 i
┌ 2q , если N-4 q ¤< (q + ─)·4
│ i i i 4
q = ┤
i-1 │ 2q +1 в противном случае (...Є...)
└ i
i
Обозначим M =N-4 q ¤. В частности, M =N.
i i imax
Получаем уже две формулы:
1 i
┌ 2q , если M < (q + ─)·4
│ i i i 4
q = ┤
i-1 │ 2q +1 в противном случае (...Є...)
└ i
1 i
┌ M , если M < (q + ─)·4
│ i i i 4
M = ┤ 1 i
i-1 │ M -(q + ─)·4 в противном случае.
└ i i 4
Формулы для M доказываются просто.
верхняя:
i-1 i-1 i
M = N-4 q ¤ = N-4 (2q )¤ = N-4 q ¤
i-1 i-1 i i
что и является записью для M .
i
нижняя:
i-1 i-1
M = N-4 q ¤ = N-4 (2q +1)¤ =
i-1 i-1 i
i-1 i i i-1
= N-4 (4q ¤+4q +1) = N-4 q ¤-4 q -4 =
i i i i
i i-1 1 i
= M -4 q -4 = M -(q + ─)·4
i i i i 4
что и требовалось доказать.
1
Теперь обозначим R = q + ─. Тогда q =[R ].
i i 4 i i
1
И, в частности, R = ─.
imax 4
Получаем выражения:
1 i
┌ 2[R ]+ ─, если M < R ·4
│ i 4 i i
R = ┤ 5
i-1 │ 2[R ]+ ─, в противном случае.
└ i 4
i
┌ M , если M < R ·4
│ i i i
M = ┤ i
i-1 │ M -R ·4 , в противном случае.
└ i i
Теперь избавимся от операции возведения
i N
в степень. Обозначим NS = M /4 = ─── -q ¤,
i i 4^i i
N
в частности, NS = ──────.
imax 4^imax
NS в общем случае дробное, не более 2·imax
цифр после запятой.
Остаются две формулы,которыми уже можно
пользоваться:
1
┌ 2[R ]+ ─, если NS < R
│ i 4 i i
R = ┤ 5
i-1 │ 2[R ]+ ─, в противном случае.
└ i 4
┌ 4·NS , если NS < R
│ i i i
NS = ┤
i-1 │ 4(NS -R ), в противном случае.
└ i i
Верхнее выражение для NS выведено так:
i-1 i-1 i
NS = M /4 = M /4 = 4·M /4 = 4·NS .
i-1 i-1 i i i
Нижнее - так:
i-1 i i-1
NS = M /4 = (M -4 R )/4 =
i-1 i-1 i i
i-1
= M /4 -4·R = 4·NS -4·R = 4(NS -R ).
i i i i i i
Итак, последовательность вычислений та-
кая: присваиваем NS и R исходные значения
(NS = N и R = 1/4 соответственно);
imax imax
крутим цикл imax раз, каждый раз сравнивая
NS и R. Вычитаем,если вычитается. Сдвигаем
NS на 2 бита,а [R] на 1 бит влево. Дробная
часть R всегда равна 1/4. На каждом шаге
получаем очередной бит результата. Исполь-
зуются только операции сдвига и сложения/
вычитания.
После окончания вычислений результат
лежит в целой части R. В NS лежит остаток
квадратного корня.
Теперь выясним, какая разрядная сетка
нам необходима для вычислений. Выше уже
указывалось, что N имеет 2imax цифр, соот-
ветственно, NS имеет 2imax цифр после за-
пятой. А сколько до запятой? Оказывается,
всего лишь imax цифр. Дело в том, что если
бы целая часть NS (остатка квадратного ко-
рня) оказалась больше, то... В общем,двумя
главами выше вы уже узнали,что остаток ко-
рня не может превышать 2q, где q - сам ко-
рень.Иначе бы извлёкся корень не q, а q+1.
А разрядность корня у нас всего imax бит.
Соответственно, разрядность целой части NS
не превышает разрядности 2q, т.е. imax+1
бит. Причём за переполнением после оконча-
ния последнего прохода цикла можно не сле-
дить,а максимальное значение 2q до после-
i
днего прохода не превышает q, т.е. одним
битом NS можно пожертвовать. В нашем слу-
чае (если N имеет 16 бит), нужно 8 бит NS
в целой части и 16 в дробной. Для R нужно
8 бит в целой части и ещё 2 фиктивных би-
та, чтобы хранить дробную часть 1/4.
Если увеличить разрядную сетку, можно
продолжать вычисления и получать цифры ре-
зультата после запятой!
Надеюсь,я доступно изложил свои рассуж-
дения... Иначе просто не стоило старать-
ся... К слову, этот же алгоритм в десятич-
ном виде использовался на счётах,когда они
ещё были популярны :)
Как видно, сложность этого алгоритма не
выше, чем у операции деления, даже проще,
так как число проходов цикла вдвое меньше.
Существуют и другие методы извлечения
квадратного корня, точные и неточные. Их
описывали в своих статьях Cardinal/BDA,
Baze/3SC, Dark Raven/Magic org. Попробуйте
найти эти статьи. Например, интересно, что
если число хранится в логарифмированном
виде, его всего лишь надо разделить на два
:)
