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. Попробуйте 
найти эти статьи. Например, интересно, что 
если  число  хранится  в логарифмированном 
виде, его всего лишь надо разделить на два 
:)