Inferno
#04
22 июня 2003 |
|
For Coderz - Алгоритм нахождения целой части квадратного корня.
Квадратный корень. (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. Попробуйте найти эти статьи. Например, интересно, что если число хранится в логарифмированном виде, его всего лишь надо разделить на два :)
Другие статьи номера:
Похожие статьи:
В этот день... 12 декабря