= 3; = 5; 1)

else s else s

= 4; = 6;

if (x <= 24)

if (x <= 3) return (x + 3) else if (x <= 8) return 2; else return (x >> 4) + 3; else if (x <= 288)

if (x <= 80) s else if (x <= 1088) x else if (x <= 1025*1025 -if (x <= 257*257 - 1)

if (x <= 129*129 - 1) else if (x <= 513*513 - 1) else if (x <= 4097*4097 - 1)

if (x <= 2049*2049 - 1) s else if (x <= 16385*16385 -if (x <= 8193*8193 - 1) else if (x <= 32769*32769 -gO = 1 « s; go = 2**s

Далее все так же, как и в листинге 11.1

S = 7; else S = 9; else

= 11; else s

= 8; = 10;

= 12;

s = 13; else s = 14; 1) s = 15; else s = 16;

Время выполнения алгоритма из листинга 11.1 на компьютере с базовым набором RISC-команд в наихудшем случае составляет около 26 + (D + 6)n тактов, где D - время

выполнения команды деления в циклах, an - количество итераций цикла while. Время выполнения алгоритма из листинга 11.2 составляет в наихудшем случае 27+(£) + 6)п

тактов в предположении (в обоих случаях), что команда ветвленш выполняется за один такт. В приведенной ниже таблице показаны среднее количество итераций цикла в обоих алгоритмах для равномерно распределенного в указанных диапазонах значения х.

Если считать, что время деления равно 20 тактам, а х равномерно распределено от О до 9999, то оба алгоритма требуют около 81 такта процессорного времени.

Бинарный поиск

Поскольку алгоритм, основанный на методе Ньютона, начинается с бинарного поиска начального приближения, почему бы не воспользоваться бинарным поиском самого квадратного корня? Этот способ может начать работу с двух границ: по всей видимости, О и 2*. Затем в качестве приближения квадратного корня берется середина между границами. Если квадрат этого значения больше аргумента х, то средняя точка становится новой верхней границей; в противном случае средняя точка становится новой нижней границей. Итерации продолжаются до тех пор, пока нижняя и верхняя границы не будут отличаться на 1. При этом нижняя граница будет представлять собой искомый результат.

Таким образом можно избежать деления, но потребуется значительное количество умножений- 16, если в качестве начальных границ воспользоваться значениями О и 2*. (Этот метод на каждой итерации вычисляет по одному биту точного значения.) В листин-

ге 11.3 показана одна из вариаций данного алгоритма, в которой используются несколько улучшенные по сравнению со значениями О и 2" начальные границы. Кроме того, процедура в листинге 11.3 для большинства RISC-компьютеров экономит один такт на цикл, изменяя аяЬтак, что в программе используется сравнение Ь>а вместо Ь-аИ.

Листинг 11.3. Бинарный поиск целочисленного квадратного корня

int isqrt(unsigned х)

unsigned а, b, m; Границы и средняя точка а = 1;

b = (х >> 5) +8; См. пояснения в тексте if (b > 65535) b = 65535; do {

m = (a+b) >> 1; if (m*m > x) b = m - 1; else a = m + 1; } while (b >= a); return a-1;

В начале каждой итерации должны выполняться предикаты а< -Jx +1 к Ь> Начальное значение b должно быть легко вычислимым и близким к -Jx

. Вот некоторые из возможных начальных значений: дг, хА + \, х + 8 + 2, х + 16+4, х-!-32 + 8, x-i-64+ 16 и т.д. Выражения в начале этого списка лучше подходят для небольших значений X, выражения ближе к концу списка- для х побольше. (Значение x-f-2 + l также применимо, но не имеет смысла, так как выражение х-4 + 1 везде дает лучшую или такую же границу.)

Семь вариаций приведенного в листинге 11.3 алгоритма можно получить, заменяя а на а + \,Ьш f»-l, используя вместо m = (fl + f») + 2 выражение т = {а + Ь + \)2 или комбинируя перечисленные подстановки.

Время выполнения процедуры из листинга 11.3 составляет примерно 6+(М+7.5)п

тактов, где М - время вьшолнения умножения в тактах, an - количество выполняемых итераций цикла. Приведенная ниже таблица содержит средние количества выполнения циклов для равномерно распределенных в указанных интервалах значений дг.

Если считать, что время вьшолнения умножения - 5 тактов, а х равномерно распределено от О до 9999, то время вьшолнения алгоритма составит около 94 тактов. Максимальное время выполнения (л =16) равно 206 тактам.

Если компьютер оснащен командой подсчета ведущш нулевых битов, то начальные границы можно определить следующим образом:

b = (1 « (33 - nlz(x))/2) - 1; а = (b + 3)/2;

ИЛИ i> = 2*""°"*-1. Это очень хорошие границы для небольших д: (так, при 0<д:<15 требуется только одна итерация), но для больших значений х дает только небольшое улучшение по сравнению с границами из листинга 11.3. При х из диапазона от О до 9999 среднее количество итераций примерно равно 5.45, что дает время выполнения около 74 тактов (при использовании тех же предположений, что и раньше).

Аппаратный алгоритм

Существует алгоритм, использующий при вычислении квадратного корня сдвиги и вычитания, который очень похож на алгоритм аппаратного деления из листинга 9.2. Будучи аппаратно встроен в 32-битовый компьютер, этот алгоритм использует 64-битовый регистр, инициализируемый аргументом х, за которым следует 32 нулевых бита. При каждой итерации 64-битовый регистр смещается на два бита влево, а текущий результат у (изначально равный 0) - на один бит. Затем из левой половины 64-битового регистра вычитается значение 2у +1. Если результат вычитания неотрицателен, он заменяет левую половину 64-битового регистра, а к значению у прибавляется 1 (для этого не требуется сумматор, так как у в этот момент заканчивается нулевым битом). Если результат вычитания отрицателен, то и 64-битовый регистр, и у остаются неизменными. Итерации выполняются 16 раз.

Этот алгоритм был описан в 1945 году [36].

Пожалуй, неожиданным результатом оказывается то, что время работы этого процесса составляет примерно половину времени, необходимого для выполнения деления 64-5-32=>32 при помощи упомянутого аппатного алгоритма, так как в этом случае вьшолняется в два раза меньше итераций примерно той же сложности, что и в алгоритме деления.

При программном использовании данного алгоритма, вероятно, лучше всего избежать использования сдвига в двойном слове, которое требует около четырех команд сдвига. Алгоритм из листинга 11.4 [19] осуществляет это путем сдвига у и маскирования бита т справа. В среднем ему требуется выполнение 149 базовых RISC-команд. Два выражения у I m могут быть заменены сложением у+т.

Листинг 11.4. Аппаратный алгоритм поиска целочисленного квадратного корня

int isqrt(unsigned х)

unsigned m, у, b; m = 0X40000000; у = 0;

while(m != 0) { Выполняется 16 раз b = у I m; У = У » 1; if (X >= b) {

X = X - b; у = у I m;

m = m >> 2;

return У;