О вычислении арифметических корней

С.В.Савич

Посвящаю

Виктору Ивановичу Поленякину

При решении практических расчётных задач часто возникает необходимость вычисления

корней разной степени. Обычно при программировании на ЭВМ для этой цели используются

итерационные

методы

или

зависимости

логарифмических

и

экспоненциальных

функций.

Аналитические

методы

последовательных

приближений,

применяемые

при

вычислении

арифметических

корней,

имеют

универсальный

характер,

однако

обладают

некоторыми

недостатками,

одним

из

которых

является

зависимость

времени

вычисления

от

величины

аргумента

и

от

выбора

первого

приближения.

Значительно

лучшие

характеристики

при

вычислении, например, квадратного корня, показывает метод, основанный на свойстве суммы

членов арифметической прогрессии нечётных чисел, который можно отнести к группе методов

"цифра за цифрой". Особенность этого метода заключается в получении на каждом циклически

повторяющемся шаге одной верной цифры результата.

В ходе анализа данного метода

был сделан вывод о возможности использовать его

концепцию

для

вычисления

корней

n-й

степени,

и

проведено

численное

исследование

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

последовательность нечётных чисел, используемая для вычисления квадратного корня - это не

только арифметическая прогрессия с шагом 2, но, - главное в этой идее, - также ряд первых

конечных разностей (далее - конечные разности) квадратичной функции с единичным шагом

изменения аргумента.

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

Составим таблицу конечных разностей степенной функции y = x

n

, например, с показателем 1; 2;

3; 5; с единичным шагом на интересующем нас интервале изменения аргумента от 0 до 9. Эта

таблица формируется следующим образом. В столбце 3 записываются значения

y=x

1

, а в

столбце 4 - разность между соседними по вертикали значениями столбца 3. По такому же

принципу заполняются пары столбцов 5-6, 7-8, 9-10.

Таблица. Конечные разности функции y=x

n

.

i

x

Показатель степени

n=1

n=2

n=3

n=5

y=x

1

Δ(x

1

)

y=x

2

Δ(x

2

)

y=x

3

Δ(x

3

)

y=x

5

Δ(x

5

)

1

2

3

4

5

6

7

8

9

10

0

0

0

1

0

1

0

1

0

1

1

1

1

1

1

3

1

7

1

31

2

2

2

1

4

5

8

19

32

211

3

3

3

1

9

7

27

37

243

781

4

4

4

1

16

9

64

61

1024

2101

5

5

5

1

25

11

125

91

3125

4651

6

6

6

1

36

13

216

127

7776

9031

7

7

7

1

49

15

343

169

16807

15961

8

8

8

1

64

17

512

217

32768

26281

9

9

9

81

729

59049

Определение: Для функции y = f(x) обозначим Δx - постоянное конечное приращение

(шаг). Тогда Δy = Δf(x) = f(x + Δx) - f(x) называется конечной разностью первого порядка

функции y = f(x).

В соответствии с этим определением выпишем выражения для вычисления конечных

разностей степенной функции y = x

n

с показателем 1; 2; 3; 5; с единичным шагом

Δ(x

1

) = (x + 1)

1

– x

1

= 1

(1)

Δ(x

2

) = (x + 1)

2

– x

2

= 2x + 1

(2)

Δ(x

3

) = (x + 1)

3

– x

3

= 3x

2

+ 3x + 1

(3)

Δ(x

5

) = (x + 1)

5

– x

5

= 5x

4

+ 10x

3

+ 10x

2

+ 5x + 1

(4).

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

разностей функции y = x

n

, начиная с единицы, равна числу слагаемых, возведённому в степень

n;

из

этого

следует,

что,

вычитая

последовательно,

также

начиная

с

единицы,

конечные

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

целую часть корня. При этом вычитания производятся до тех пор, пока разность от вычитания

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

разность,

получаемая

в

результате

последовательного

вычитания

из

подкоренного

числа

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

Для упрощения дальнейших рассуждений примем пока, что подкоренное число - целое.

Попутно отметим, что все цифровые вычислительные машины, независимо от сложности, на

низком уровне всегда обрабатывают только целые числа.

Рассмотрим, из методических соображений, процесс вычисления "вырожденного" корня

1-й степени из, например, 123. Поскольку конечная разность (1) функции y = x

1

не зависит от x

и при единичном шаге всегда равна 1, то вычисление этого корня сводится к многократному

вычитанию единицы из подкоренного числа. В нашем примере можно сделать 123 вычитания.

На самом деле количество вычитаний можно значительно уменьшить, если вычитать конечные

разности поразрядно со сдвигом. Для этого разобьём подкоренное число на группы, по одной

цифре в группе, так как показатель корня равен 1, и будем вычитать из них конечные разности,

начиная с крайней левой группы, до тех пор, пока остаток от вычитания при обработке каждой

группы не станет меньше следующего вычитаемого. После каждого цикла вычитаний число

вычитаний будет очередной цифрой в значении корня. По своей сути этот алгоритм есть не что

иное, как всем известный способ деления "уголком" на единицу.

Вычисление корней более высоких степеней также основывается на алгоритме деления,

но

только

с

тем

отличием,

что

теперь

значение

вычитаемых

конечных

разностей

будет

переменной величиной, и его надо вычислять по формуле, соответствующей степени корня.

Для квадратного корня конечные разности вычисляются с помощью выражения (2), где

из-за дискретности изменения аргумента можно произвести целочисленную замену x = i

Δ(x

2

) = 2i + 1

(5)

Для примера, вычислим с помощью вычитания конечных разностей квадратный корень из

1225. Вычитаем из подкоренного числа значения 2i + 1 при i, изменяющемся от 0 до n с шагом

1,

до

тех

пор,

пока

остаток

не

станет

меньше

следующего

вычитаемого.

Если

вычитать

непосредственно, то придётся выполнить 35 вычитаний. Таким образом, квадратный корень из

1225 равен 35 - если подсчитывать количество вычитаний, или последнему значению i = n,

увеличенному на единицу (34 + 1 = 35). Для уменьшения количества вычитаний и увеличения

быстродействия применим поразрядное вычитание со сдвигом. Для этого разобьём исходное

число влево и вправо (при наличии дробной части) от десятичной запятой на пары - своего рода

квадратичные

разряды,

так

как

показатель

корня

2,

и

будем

вычитать

из

них

конечные

разности, начиная

с крайней

левой

группы, до

тех

пор,

пока

остаток

от

вычитания

при

обработке каждой группы не станет меньше следующего вычитаемого.

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

(состоящего из нескольких цифр) корня значение i тоже будет многозначным. Для того, чтобы

отделить вычисляемую цифру от уже найденных ранее цифр корня, представим i в виде суммы

i = 10A + i

d

,

(6)

где

A

-

целое

число,

составленное

из

уже

вычисленных

цифр

корня,

а

i

d

-

значение

одноразрядного счётчика вычисляемого разряда (индекс d - сокращ. digit).

Подставляя (6) в (5), получаем выражение

Δ(x

2

) = 20A + 2i

d

+ 1

(7)

Изложим подробнее порядок вычисления квадратного корня.

1. Вначале A = 0, так как цифры корня неизвестны. Из первой группы цифр слева

вычитаем

последовательно

конечные

разности,

вычисленные

по

формуле

(7)

при

i

d

,

увеличивающемся от 0 до n с единичным шагом, до тех пор, пока остаток от вычитания не

станет меньше следующего вычитаемого. Полученный остаток обозначим D. Если он равен

нулю, и справа от этой группы все цифры нули, то корень извлечён точно (первая ненулевая

цифра корня равна i

d

+ 1).

В нашем примере можно вычесть конечные разности три раза: 12 - 1 = 11; 11 - 3 = 8;

8 - 5 = 3; 3 < 7 (7 - это следующее вычитаемое), значит, вычитания прекращаем, при этом i

d

= 2;

D = 3. Первая цифра корня будет B = i

d

+ 1 = 3, также она равна количеству вычитаний.

2. Вычисляем выражение 10A + B. Это будет новое значение A. Теперь к остатку D справа

припишем следующую группу C из двух цифр подкоренного числа, получим число 100D + C. В

нашем примере, приписывая следующие две цифры, получим число 325. Из числа 100D + C

снова

начинаем

вычитать

последовательно

конечные

разности,

предварительно

выполнив

выравнивание

порядков.

Эти

конечные

разности

вычисляются

по

формуле

(7)

с

новым

значением A (в нашем примере A = 3), при i

d

, увеличивающемся от 0 до n с единичным шагом.

Вычитания производятся до тех пор, пока остаток от вычитания не станет меньше следующего

вычитаемого. Если этот остаток равен нулю, и оставшиеся цифры в подкоренном числе - нули,

то корень извлечён точно (последняя ненулевая цифра корня равна i

d

+ 1).

В нашем примере производим следующие действия: 325 - 61 = 264; 264 - 63 = 201;

201 - 65 = 136; 136 - 67 = 69; 69 - 69 = 0; остаток от вычитания D = 0, поэтому вычитания

прекращаем, при этом i

d

= 4. Следующая цифра корня будет B = i

d

+ 1 = 5 (или, что то же самое,

5 вычитаний). В общем случае, если остаток D не равен нулю, то повторяем алгоритм с п. 2.

Для

разбираемого

нами

примера

после

выполнения

п.

2

переменные

будут

иметь

следующие значения: D = 0; B = 5. Корень вычислен точно и равен 10A + B = 35.

Примечание: Если в начале вычисления следующей цифры первое вычитаемое (конечная

разность, вычисленная по формуле (7) при i

d

= 0) больше, чем 100D + C, то следующая цифра

корня B = 0, так как ни одного вычитания сделать нельзя. Тогда число 100D + C будет остатком

D; переходим к выполнению п. 2.

Для того, чтобы подчеркнуть сходство с алгоритмом деления "уголком", запишем все

вычисления, как принято при делении, в столбик (пояснения - в фигурных скобках)

12'25,

-1

{(7) при A=0; i

d

=0; символ "-" - операция вычитания}

11

-3

{i

d

=1}

8

-5

{i

d

=2}

325

{3<7(следующ. конечн.разн.КР); 3 вычит.; приписываем 25}

-61

{(7) при A=3; i

d

=0}

264

-63

{i

d

=1}

201

-65

{i

d

=2}

136

-67

{i

d

=3}

69

-69

{i

d

=4}

0

{D=0; стоп; 5 вычитаний}

Ответ: 35

Здесь с целью упрощения не был описан процесс определения положения десятичной

запятой,

т.е.

порядка

результата.

Учитывая,

что

каждый

квадратичный

разряд

аргумента

соответствует одному разряду вычисленного корня, ясно, что решение этого вопроса не должно

иметь принципиальных затруднений.

Анализируя

этот

алгоритм,

можно

заметить,

что

для

вычисления

каждой

конечной

разности приходится всякий раз вычислять выражение (7). Это замедляет и усложняет процесс

вычислений. Упростить вычисления можно, приняв во внимание следующее соображение -

вычислять не саму конечную разность, а приращение к предыдущему её значению. Эта задача

разделяется на две подзадачи:

1. Вычисление приращения в процессе получения цифр одного разряда;

2. Вычисление приращения при переходе к обработке следующего разряда.

Эта проблема была успешно решена следующим образом. Для решения первой подзадачи

запишем правую часть выражения (7), при исходном i

d

и при i

d

, увеличенном на единицу

Δ(x

2

)

n

= 20A + 2i

d

+ 1

(7.0)

Δ(x

2

)

n+1

= 20A + 2(i

d

+ 1) + 1

(7.1)

Вычитая (7.0) из (7.1), получаем выражение для вычисления приращения, которое нужно

прибавить к конечной разности, чтобы получить конечную разность при следующем i

d

ΔΔ(x

2

) = (20A + 2(i

d

+ 1) + 1) - (20A + 2i

d

+ 1) = 2

(8)

Собственно говоря, это приращение есть конечная разность второго порядка функции

y = x

2

с единичным шагом. Как и следовало ожидать, она является константой и равна 2.

Решение второй подзадачи требует некоторых пояснений. В конце цикла вычитаний

конечных

разностей

при

вычислении

очередной

цифры

корня

имеем

значения

A

и

i

d

,

применённые в этом цикле. В начале вычисления следующей цифры будем иметь следующие

параметры - A

нов

= 10A + i

d

+ 1, и i

d нов

= 0. Запишем начальное выражение (7) для "нового"

цикла вычитаний, используя только что определённые значения A

нов

и i

d нов

Δ(x

2

)

нов

= 20(10A + i

d

+ 1) + 1

(9)

Вычитая (7) из (9), получаем выражение для вычисления приращения, которое нужно

прибавить к последней конечной разности только что завершённого цикла вычитаний, чтобы

получить первое вычитаемое следующего цикла (индекс dd - "между разрядами")

Δ(x

2

)

dd

= 180A + 18i

d

+ 20

(10)

Вычисление квадратного корня начинается с вычитания из первой левой группы единицы,

затем,

для

получения

следующих

конечных

разностей

используем

выражения

(8)

и

(10),

учитывая i

d

. Составление подробного алгоритма не составляет труда, с учётом вышеописанной

схемы и вычислительного примера, и, для экономии места, здесь не приводится.

Теперь

рассмотрим

процесс

вычисления

кубического

корня.

Для

кубического

корня

конечные разности вычисляются с помощью выражения (3), используя подстановки x = i и (6)

Δ(x

3

) = 3i

2

+ 3i + 1 = 3(10A + i

d

)

2

+ 3(10A + i

d

) + 1

Затем

производим

алгебраические

преобразования

и

факторизацию,

получаем

в

результате конечное выражение

Δ(x

3

) = 3(10A + i

d

)(10A + i

d

+ 1) + 1

(11)

Порядок вычисления корня 3-й степени, по существу, аналогичен подробно разобранному

алгоритму вычисления квадратного корня. Отличия в описании алгоритма заключаются в

следующем:

1.

Подкоренное

число

разбиваем

влево

и

вправо

от

запятой

на

тройки (кубические

разряды);

2. Конечные разности вычисляем по формуле (11);

3. Приписывая группу цифр C, состоящую из 3-х цифр, получаем число 1000D + C;

4. Примечание изменяется с учётом использования формулы (11) и разбиения на тройки.

Покажем

процесс

вычисления

кубического

корня,

записанный

в

формате

операции

деления. Все расчёты примеров далее выполнены на программируемом калькуляторе HP-41C

фирмы "Хьюлетт-Паккард". Для примера, возьмём число 56789,321

056'789,321

-1

{(11) при A=0; i

d

=0}

55

-7

{i

d

=1}

48

-19

{i

d

=2}

29789

{29<37(следующ. КР); 3 вычит; приписываем 789}

-2791

{(11)при A=3; i

d

=0}

26998

-2977

{i

d

=1}

24021

-3169

{i

d

=2}

20852

-3367

{i

d

=3}

17485

-3571

{i

d

=4}

13914

-3781

{i

d

=5}

10133

-3997

{i

d

=6}

6136

-4219

{i

d

=7}

1917321

{1917<4447(след. КР); 8 вычит.;приписываем 321;зпт}

-434341

{(11)при A=38; i

d

=0}

1482980

-436627

{i

d

=1}

1046353

-438919

{i

d

=2}

607434

-441217

{i

d

=3}

166217000

{166217<443521(след. КР); 4 вычит.; добавл. 000}

-44248321

{(11)при A=384; i

d

=0}

121968679

-44271367

{i

d

=1}

77697312

-44294419

{i

d

=2}

33402893

{33402893<44317477(след. КР); 3 вычит.; и т.д.}

Ответ 38,43...

Из соображений упрощения расчётов запишем выражения для приращений к конечным

разностям. Формулы для вычисления приращений выводятся на тех же принципах, что и для

квадратного корня. Приведём их без выкладок, учитывая ограниченный объём статьи.

1.

Выражение

для

вычисления

приращения,

которое

нужно

прибавить

к

конечной

разности, чтобы получить конечную разность при следующем i

d

ΔΔ(x

3

) = 6(10A + i

d

+1)

(12)

2.

Выражение

для

вычисления

приращения,

которое

нужно

прибавить

к

последней

конечной

разности

только

что

завершённого

цикла

вычитаний,

чтобы

получить

первое

вычитаемое следующего цикла

Δ(x

3

)

dd

= 33(10A + i

d

+ 1)(90A + 9i

d

+ 10)

(13)

Вычисление кубического корня начинается с вычитания из первой левой группы единицы,

затем, для получения следующих конечных разностей используем приращения, вычисленные

по выражениям (12) и (13). Вычитания производятся до тех пор, пока остаток не станет меньше

следующего

вычитаемого.

Подробное

описание

здесь

не

приводится;

оно

может

быть

составлено на основе алгоритма вычисления квадратного корня после разбора приведённого

примера.

Теперь, учитывая всё вышесказанное о квадратных и кубических корнях, выведем по тем

же правилам аналогичные выражения для вычисления корней 5-й степени. Здесь конечные

разности вычисляются с помощью выражения (4), с использованием подстановок x = i и (6)

Δ(x

5

)=5i

4

+10i

3

+10i

2

+5i+1=5(10A+i

d

)

4

+10(10A+i

d

)

3

+10(10A+i

d

)

2

+5(10A+i

d

)+1=

=5(10A+i

d

)((10A+i

d

)((10A+i

d

)(10A+i

d

+2)+2)+1)+1

(14)

Вычисление корня 5-й степени принципиально не отличается от вычисления кубического

корня. Для примера, вычислим корень 5-й степени из 7167031,46875. Сначала разбиваем число

на

группы

по

пять

цифр

(пятистепенные

разряды),

вправо

и

влево

от

запятой,

получим

00071'67031,46875. Далее вычитаем из первой слева группы конечные разности, вычисленные с

использованием

выражения

(14).

Условие

окончания

вычитаний

и

перехода

к

обработке

следующего разряда такое же, как при извлечении квадратного и кубического корней. После

окончания цикла вычитаний к остатку D приписывается следующая группа из пяти цифр,

получаем

100000D

+

C,

присваивается

новое

значение

A

(A

нов

=

10A

+

B),

и

процесс

продолжается до достижения необходимой точности:

00071'67031,46875

-1

{(14) при A=0; i

d

=0}

70

-31

{i

d

=1}

3967031

{39<211(след. КР); 2 вычитания; приписываем 67031}

-884101

{(14) при A=2; i

d

=0}

3082930

-1069531

{i

d

=1}

2013399

-1282711

{i

d

=2}

73068846875

{730688<1526281; 3 вычит. приписываем 46875}

-14114250151

{(14) при A=23; i

d

=0}

58954596724

-14360780281

{i

d

=1}

44593816443

-14610525961

{i

d

=2}

29983290482

-14863515031

{i

d

=3}

15119775451

-15119775451

{i

d

=4}

0

{стоп; 5 вычитаний}

Ответ 23,5

Из тех же соображений, как при рассмотрении кубического корня, приведём без выкладок

формулы приращений для вычисления конечных разностей.

Приращение в процессе вычисления цифр одного разряда

ΔΔ(x

5

) = 10(10A + i

d

+1)((200A + 40(i

d

+ 1))A + 2i

d

(i

d

+ 2) + 3)

(15)

Приращение при переходе к обработке следующего разряда

Δ(x

5

)

dd

=55(10A+i

d

+1)(90A+9i

d

+10)(A(10100A+10(202i

d

+211))+((101i

d

+211)i

d

+111))

(16)

Подведём некоторые итоги. Сразу видно, что с увеличением степени извлекаемого корня

коэффициенты многочленов становятся очень большими, структура выражений усложняется,

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

статье

все

рассуждения

и

вывод

формул

проводились,

опираясь

на

привычную

всем

десятичную

систему

счисления.

Однако,

так

как

машинное

вычисление

основано

на

использовании

двоичной

арифметики,

приведём,

для

сравнения,

перевод

только

что

показанного расчёта корня 5-й степени в двоичный формат.

Основным выражением, определяющим формулу для вычисления конечных разностей с

учётом уже вычисленных цифр корня, является (6). Так как при двоичном счёте i

d

никогда не

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

i = 10

(bin)

A = 2

(dec)

A

(6b)

Подставляя (6b) вместо x в формулу (4) и применяя схему Горнера для многочленов,

получаем выражение для вычисления конечных разностей функции y = x

5

в двоичной системе

счисления с учётом найденных цифр корня (коэффициенты показаны в десятичном виде)

Δ(x

5

) = (((A + 1)2A + 1)4A + 1)10A + 1

(14b)

Для нашего примера вначале переводим подкоренное число 7167031,46875 в двоичную

систему и разбиваем

его влево и вправо от запятой на группы

по пять

двоичных цифр

(пятистепенные разряды). Для наглядности эти группы будут изображены в виде десятичных

эквивалентов в скобках. Вычитаемое вычисляется по (14b). После каждого вычитания при

переходе к обработке следующего разряда к остатку D приписывается следующая группа C из

пяти двоичных цифр, получаем 100000

(bin)

× D + C = 32

(dec)

× D + C, вычисляется новое значение

A (A

нов

= 10

(bin)

× A + B = 2

(dec)

× A + B; где B = 1, если вычитание выполнялось, или B = 0, если

не выполнялось), и расчёт продолжается до достижения требуемой точности.

7167031,46875

(dec)

= 11011010101110000110111,01111

(bin)

00110'11010'10111'00001'10111,01111

(6)'(26)'(23)'(1)'(23),(15)

6

-1

{(14b) при A=0; 1 вычитание}

5

186

{(5×32+26)<211(след.КР (14b) при A=1);0 вычит.}

5975

{186×32+23}

-2101

{(14b) при A=10

(bin)

=2

(dec)

; 1 вычитание}

3874

123969

{3874×32+1}

-61051

{(14b) при A=101

(bin)

=5

(dec)

; 1 вычитание}

62918

2013399

{62918×32+23}

-1282711

{(14b) при A=1011

(bin)

=11

(dec)

; 1 вычит.; зпт}

730688

23382031

{730688×32+15}

-23382031

{(14b) при A=10111

(bin)

=23

(dec)

; 1 вычитание}

0

{стоп}

Ответ: 10111,1

(bin)

= 23,5

(dec)

Применяя этот же подход, можно вывести формулы конечных разностей для вычисления

в двоичной системе кубических и квадратных корней. Вот эти выражения

для кубического корня Δ(x

3

) = (2A + 1)6A + 1

(11b)

для квадратного корня Δ(x

2

) = 4A + 1

(7b)

Таким

образом,

применение

двоичной

арифметики

значительно

упрощает

процесс

вычислений.

Приведённые

примеры

расчётов

показывают,

что,

в

отличие

от

методов

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

-

если

вычисляемый

корень

представляет

собой

конечное

рациональное

число,

он

будет

вычислен точно за конечное число шагов.

Отметим также, что применяемые аналитические и функциональные методы вычисления

корней

основаны

на

использовании

операций

с

плавающей

запятой,

что

не

способствует

повышению

быстродействия.

Разработанные

с

помощью

описанной

в

этой

статье

универсальной

методики

целочисленные

алгоритмы

вычисления

арифметических

корней

разной степени позволяют ускорить выполнение инженерных расчётов.