Работа этого алгоритма похожа на школьный метод поиска квадратного корня. Ниже показан пример поиска 4\19 в соответствии с этим алгоритмом на 8-битовом компьютере.

1011 ООН хо Изначально, х = 179 (ОхВЗ) - 1 Ы

Ь2 0010 0000 у2

Х2 ООН 0000 у2 ЬЗ 0001 1000 уз

0000 1010 Х4 0000 1101 у4

Результат вычислений равен 13; в регистре х остается остаток 10.

При использовании обычного трюка со знаковым сдвигом вправо на 31 бит можно избежать выполнения проверки if x>=b. Можно доказать, что старший бит Ь всегда нулевой (на самом деле * < 5 • 2), что упрощает предикат х>=ь (см. раздел 2.11). В результате группу операторов if можно заменить следующими:

t = (int)(х I -(х - b)) » 31; -1, если х >= Ь, иначе О X = X - (Ь Ь t) ; у = у I (т ь t) ;

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

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

а кроме того, как обычно, возникает проблема поиска хорошего начального приближения ль-

Однако существует аппаратный алгоритм, приведенный в листинге 11.5, который также легко реализуется программно.

Листинг 11.5. Аппаратный алгоритм поиска целочисленного кубического корня

int icbrt(unsigned х)

int S;

unsigned у, b;

S = 30; у = 0;

while(s >= 0) { 11 итераций У = 2*y;

b = (3*y*(y +1) +1) « s; s = s - 3; if (X >= b) { X = X - b;

У = У + 1;

return у;

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

Этого умножения легко избежать путем определенной оптимизации. Введем еще одну беззнаковую переменную у2, которая будет содержать значение и обновляться при каждом изменении у. Непосредственно перед присвоением у = о вставим присвоение у2 = о, а перед у = 2*у- присвоение у2 = 4*у2. Присвоение значения переменной b заменим выражением Ь= (з* (у2+у)+1) «s, а непосредственно перед у = у+1 вставим у2 = у2+2*у+1. Полученная в результате программа не содержит умножений, за исключением умножения на малую константу, которое можно заменить сдвигом и сложением. Эта программа содержит три сложения с 1, которые можно заменить командой или. Такая программа будет работать быстрее, если только умножение в вашем компьютере не выполняется за два такта (или быстрее).

Будьте осторожны: в [19] указывается, что этот код неработоспособен при непосредственной адаптации его для 64-битовых машин. Присвоение b в этом случае может вызвать переполнение. Этой проблемы можно избежать, убрав сдвиг на s позиций в присвоении Ь, вставив после присвоения b выражение bs = b«s и заменив две строки if (x>=b) {x=x-b... строками if (x>=bsbbb== (bs>>s) ) {x=x-bs....

в случае вычисления кубического корня метод Ньютона работает существенно хуже в силу более сложной итеративной формулы

Количество умножений, которые надо выйолнить при использовании данного метода при показателе степени п > 1, равно

[log,nj + nbits(n)-l.

Это не всегда минимально необходимое число умножений. Например, при и=27 метод бинарного разложения приводит к вычислению

что требует выполнения семи умножений. Однако если воспользоваться схемой

Вычисление л:" бинарным разложением п

Хорошо известен метод вычисления х", где и - неотрицательное целое число, использующий двоичное представление числа и. Этот метод применим к вычислению выражений вида хХХ...х, где • представляет собой любой ассоциативный оператор, такой, как сложение, умножение (включая умножение матриц) или конкатенация строк (с использованием записи (аЬ)=ababab). В качестве примера,рассмотрим вычисление y=x. Поскольку 13 в двоичном представлении имеет вид 1101 (т.е. 13 равно 8+4+1), то

x"=X***=xxX.

Таким образом, д: можно вычислить, как показано ниже.

(, <-x-

y<:-t,-t-X

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

Если показатель степени - неотрицательная целая переменная, данный метод можно реализовать в виде функции, показанной в листинге 11.6.

Листинг 11.6. Вычисление дг" бинарным разложением л

int iexp(int х, unsigned n) {

int p. У;

у = 1; Инициализация результата р = X; и переменной р while(1) {

if (n Ь 1)

У = Р*У; Если п нечетно, многким на р п = п >> 1; Позиция очередного бита п

if (п == 0)

return у; Если больше нет битов п р = р*р; Степень для очередного бита п