ГЛАВА и

НЕКОТОРЫЕ ЭЛЕМЕНТАРНЫЕ ФУНКЦИИ

11.1. Целочисленный квадратный корень

Под "целочисленным квадратным корнем" подразумевается функция 4х . Чтобы расширить область применения и чтобы избежать решения вопроса о том, как поступать с отрицательными аргументами функции, будем считать jc беззнаковой величиной: О < х < 2" -1.

Метод Ньютона

Для чисел с плавающей точкой, по сути, универсальным методом вычисления квадратного корня является метод Ньютона. Этот метод начинается с некоего значения о, которое является начальной оценкой >/а . Затем выполняется серия уточнений значения квадратного корня по формуле

8п + -

Приведенный метод имеет квадратичную сходимость, т.е. если в некоторый момент g„ имеет точность п битов, то g„+i имеет точность 2п битов. Рабочая программа должна иметь средства для определения достижения необходимой точности и прекращения вычислений.

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

Теорема. Пусть „i = +[fl/„J)/2 , где g„ - целое число, большее 0. Тогда

а) если g„> 4а ,то 4а < g,, < g„;

б) если g„= 4а ,то 4а й„i < 4а +\.

Иными словами, если у нас есть завышенная целочисленная оценка g„ величины 4а , то следующее приближение будет строго меньше предыдущего, но не меньше

4а . Таким образом, если начать со слишком большого приближения, то будет получена монотонно убывающая последовательность приближений. Если же начальное приближение g„= 4а , то следующее приближение оказывается либо равным предыдущему, либо на 1 большим. Таким образом, мы получаем способ определения того, что последовательность приближений сошлась: если начать с > >/а , то точное решение g„

будет достигнуто в тот момент, когда „,., > g„.

Случай а-О следует рассматривать отдельно в связи с тем, что при этом может возникнуть ситуация деления О на 0.

Доказательство, а) Поскольку g„ - целое число,

Sn+l -

g.+ -

g« + -

Так как g„ > >/a и g„ - целое число, g„ > >/a . Определим e следующим образом: g„ =(1 + Б)>/а . Тогда e>0 и

2g„ (l + e)fl + fl 2(l + e)>/ 2 + 2e + 8- f-

2(l+e) 2+28

2(l+e)

= g„4i = g«.i <

g+fl 2g„

2g„

= g„.i<g„. g„+i<g»-

6) Поскольку g„= ,>/a-l<g„<>/a, так что g, < fl< (g„ +1) Следовательно,

2g»

>/Jg„.,:sLgJ+i=W+i-

Сложной частью при использовании метода Ньютона для вычисления -Jx

является

получение первого приближения. Процедура в листинге 11.1 устанавливает первое приближение go равным наименьщей степени 2, которая больше или равна л[х. Например, gQ=2 для jc=4, а для х=5 в качестве первого приближения принимается go=4.

Листинг 11.1. Метод Ньютона поиска целочисленного квадратного корня

int isqrt(unsigned х) {

unsigned xl;

int s, gO, gl;

if (x <= 1) return x;

S = 1;

xl = X - 1;

if (xl > 65535)

if (xl > 255)

if (xl > 15)

if (xl > 3)

gO = 1 << s;

S = s + 8; XI = xl >> 16;} S = S + 4; xl = Xl >> 8;" s = s + 2; xl = xl >> 4; s = s + 1; } go = 2**s

gl = (gO + (X >> s)) >> 1; gi = (gO + x/g0)/2 while (gl < gO) { Повторяем, пока приближения

gO = gl; строго уменьшаются

gl = (go + (x/gO)) » 1;

return gO;

Поскольку первое приближение представляет собой степень 2, для получения нет необходимости в реальном делении - можно обойтись сдвигом вправо.

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

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

if (х <= 1) return х;

S = 16 - nlz(x - 1)/2;

Еще один вариант в случае отсутствия команды вычисления количества ведущих нулевых битов - вычисление s посредством бинарного дерева поиска. Этот метод позволяет получить несколько лучшее приближение go: наименьшую степень 2, которая больше или равна yfx . Для некоторых значений jc это дает меньшее значение go, но достаточно большое, чтобы выполнялся критерий сходимости из рассмотренной ранее теоремы. Отличие этих схем показано ниже.

Соответствующая процедура показана в листинге 11.2. Она особенно удобна при работе с малыми значениями jc (О < х < 24), так как при этом не требуется выполнение деления.

Листинг 11.2. Целочисленный квадратный корень с вычислением первого приближения путем бинарного поиска

int isqrt(unsigned х) {

int s, gO, gl; if (x <= 4224)