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

инж. С.В.Савич

Посвящаю
Виктору Ивановичу Поленякину,
Георгию Иосифовичу Кирьянову

В процессе решения практических расчётных задач довольно часто возникает необходимость вычисления корней разной степени. Обычно при программировании на ЭВМ для этой цели используются стандартные библиотечные функции вычисления логарифма и экспоненты или итерационные методы. Аналитические методы последовательных приближений, часто применяемые при вычислении арифметических корней, имеют универсальный характер, однако обладают некоторыми недостатками, одним из которых является зависимость времени вычисления от величины аргумента и от выбора первого приближения. Значительно лучшие характеристики при вычислении, например, квадратного корня, показывает метод, описанный в статье “Оригинальный метод извлечения квадратного корня” (www.hijos.ru/2012/04/25/), который можно отнести к группе методов “цифра за цифрой”. Особенность этого метода, основанного на свойстве суммы членов арифметической прогрессии нечётных чисел, заключается в получении на каждом циклически повторяющемся шаге одной верной цифры результата.

В ходе анализа данного метода возникла идея распространить его концепцию на процесс вычисления корней n-й степени, а также провести численное исследование получаемых алгоритмов. Основанием для такого подхода является то обстоятельство, что последовательность нечётных чисел, используемая для вычисления квадратного корня — это не только арифметическая прогрессия с шагом 2, но, — главное в этой идее, — также ряд первых конечных разностей (далее — конечные разности) для квадратичной функции с единичным шагом изменения аргумента.

Упомянем коротко метод “цифра за цифрой”, применяемый для “ручного” вычисления квадратного корня, заключающийся в том, что из подкоренного числа, разбитого на пары, в определённом порядке последовательно вычитается выражение (20A+B)\times B (здесь A — целое число, составленное из уже найденных цифр корня, B — искомая следующая цифра корня, определяемая подбором), при условии, чтобы значение этого выражения не превышало текущее уменьшаемое. Детально этот способ описан в учебной и справочной литературе, сейчас он имеет лишь историческое значение.

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

Определение: Для функции y = f(x) обозначим \Delta x — постоянное конечное приращение (шаг). Тогда \Delta y = \Delta f(x) = f(x + \Delta x) - f(x) называется конечной разностью первого порядка функции y=f(x).

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

Выпишем выражения для конечных разностей степенной функции y = x^n, например, с показателем 1; 2; 3; 5; с единичным шагом

(1)   \begin{equation*} \Delta (x^1) = (x + 1)^1 - x^1 = 1 \end{equation*}

(2)   \begin{equation*} \Delta (x^2) = (x + 1)^2 - x^2 = 2x + 1 \end{equation*}

(3)   \begin{equation*} \Delta (x^3) = (x + 1)^3 - x^3 = 3x^2 + 3x + 1 \end{equation*}

(4)   \begin{equation*} \Delta (x^5) = (x + 1)^5 - x^5 = 5x^4 + 10x^3 + 10x^2 + 5x + 1 . \end{equation*}

Теперь составим таблицу конечных разностей этих функций на интересующем нас интервале изменения аргумента от 0 до 9:

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

Анализируя эту таблицу, можно увидеть, что сумма последовательных конечных разностей функции y=x^n, начиная с единицы, равна числу слагаемых, возведённому в степень n; из этого следует, что, вычитая последовательно, также начиная с единицы, конечные разности из подкоренного числа, можно, путём подсчёта количества вычитаний, определить целую часть корня. При этом вычитания производятся до тех пор, пока остаток от вычитания не станет меньше следующего вычитаемого.

Для упрощения дальнейших рассуждений примем пока, что подкоренное число — целое. Попутно отметим, что все компьютеры, независимо от сложности, на низком уровне всегда обрабатывают только целые числа. Рассмотрим, из методических соображений, процесс вычисления “вырожденного” корня 1-й степени из, например, 123. Поскольку конечная разность (1) для функции y = x^1 при любом x равна 1, то вычисление этого корня сводится к многократному вычитанию единицы из подкоренного числа. В нашем примере можно сделать 123 вычитания. На самом деле количество вычитаний можно значительно уменьшить, если вычитать конечные разности поразрядно со сдвигом. Для этого разобьём подкоренное число на группы, по одной цифре в группе, так как показатель корня равен 1, и будем вычитать из них конечные разности, начиная с крайней левой группы, до тех пор, пока остаток от вычитания при обработке каждой группы не станет меньше следующего вычитаемого. После каждого цикла вычитаний число вычитаний будет очередной цифрой в значении корня. По своей сути этот алгоритм есть не что иное, как всем известный способ деления “уголком” на единицу.

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

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

(5)   \begin{equation*} \Delta (x^2) = 2i + 1 \end{equation*}

Для примера, вычислим с помощью вычитания конечных разностей квадратный корень из 1225. Вычитаем из подкоренного числа значения 2i + 1 при i, изменяющемся от 0 до n с шагом 1, до тех пор, пока остаток не станет меньше следующего вычитаемого. Если вычитать непосредственно, то придётся выполнить 35 вычитаний. Таким образом, квадратный корень из 1225 равен 35 — если подсчитывать количество вычитаний, или последнему значению i = n, увеличенному на единицу (34 + 1 = 35). Для уменьшения количества вычитаний и увеличения быстродействия применим поразрядное вычитание со сдвигом. Для этого разобьём исходное число влево и вправо (при наличии дробной части) от десятичной запятой на пары — своего рода квадратичные разряды, так как показатель корня 2, и будем вычитать из них конечные разности, начиная с крайней левой группы, до тех пор, пока остаток от вычитания при обработке каждой группы не станет меньше следующего вычитаемого.

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

(6)   \begin{equation*} i = 10A + i_d, \end{equation*}

где A — число, составленное из уже вычисленных цифр корня, а i_d — значение одноразрядного счётчика вычисляемого разряда (индекс d — сокращ. digit).

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

(7)   \begin{equation*} \Delta (x^2) = 20A + 2i_d + 1. \end{equation*}

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

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), при_id, увеличивающемся от 0 до n с единичным шагом. Вычитания производятся до тех пор, пока остаток от вычитания не станет меньше следующего вычитаемого. Если этот остаток равен нулю, и оставшиеся цифры в подкоренном числе — нули, то корень извлечён точно (последняя ненулевая цифра корня равна i_d + 1).

В нашем примере производим следующие действия:

    \[325-61=264; 264-63=201; 201-65=136; 136-67=69; 69-69=0;\]

0<71 (следующее вычитаемое), поэтому вычитания прекращаем, при этом 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.

Для того, чтобы подчеркнуть сходство с алгоритмом деления “уголком”, запишем все вычисления в столбик (пояснения — в фигурных скобках)

Ответ: 35.

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

Анализируя этот алгоритм, можно заметить, что для вычисления каждой конечной разности приходится всякий раз вычислять выражение (7). Это замедляет и усложняет процесс вычислений. Упростить вычисления можно, приняв во внимание следующее соображение — вычислять не саму конечную разность, а приращение к предыдущему её значению. Эта задача разделяется на две подзадачи:

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

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

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

    \[\Delta(x^2)_n = 20A + 2i_d + 1 \eqno(7.0)\]

    \[\Delta(x^2)_{n+1} = 20A + 2(i_d + 1) + 1\eqno        (7.1)\]

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

(8)   \begin{equation*} \Delta\Delta(x^2) = (20A + 2(i_d + 1) + 1) - (20A + 2i_d + 1) = 2. \end{equation*}

Собственно говоря, выражение (8) — это вычисление конечной разности второго порядка для квадратичной функции. Как и следовало ожидать, она является константой и равна 2.

Решение второй подзадачи требует некоторых пояснений. В конце цикла вычитаний конечных разностей при вычислении очередной цифры корня имеем значения A и i_d, применённые в этом цикле. В начале вычисления следующей цифры будем иметь следующие параметры — A_{new} = 10A + i_d + 1, и i_{d\ new} = 0. Запишем начальное выражение (7) для “нового” цикла вычитаний, используя только что определённые значения A_{new} и i_{d\ new}

(9)   \begin{equation*} \Delta(x^2)_{new} = 20(10A + i_d + 1) + 1. \end{equation*}

Вычитая (7) из (9), получаем выражение для вычисления приращения, которое нужно прибавить к последней конечной разности только что завершённого цикла вычитаний, чтобы получить первое вычитаемое следующего цикла (индекс dd — “между разрядами”)

(10)   \begin{equation*} \Delta(x^2)_{dd} = 180A + 18i_d + 20. \end{equation*}

Вычисление квадратного корня начинается с вычитания из первой левой группы единицы, затем, для получения следующих конечных разностей используем выражения (8) и (10), учитывая i_d. Составление подробного алгоритма не составляет труда, с учётом вышеописанной схемы и вычислительного примера, и, для экономии места, здесь не приводится.

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

    \[\Delta(x^3) = 3i^2 + 3i + 1 = 3(10A + i_d)^2 + 3(10A + i_d) + 1.\]

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

(11)   \begin{equation*} \Delta(x^3) = 3(10A + i_d)(10A + i_d + 1) + 1. \end{equation*}

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

1. Подкоренное число разбиваем влево и вправо от запятой на тройки (кубические разряды);

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

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

4. Текст примечания изменяется с учётом разбиения на тройки.

Покажем процесс вычисления кубического корня, записанный в столбик. Все расчёты примеров далее выполнены на программируемом калькуляторе HP-41C фирмы “Хьюлетт-Паккард”. Для примера, возьмём число 56789,321 (пояснения — в фигурных скобках)

Ответ 38.43\ldots

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

1. Выражение для вычисления приращения, которое нужно прибавить к конечной разности, чтобы получить конечную разность при следующем i_d

(12)   \begin{equation*} \Delta\Delta(x^3) = 6(10A + i_d +1). \end{equation*}

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

(13)   \begin{equation*} \Delta(x^3)_{dd} = 33(10A + i_d + 1)(90A + 9i_d + 10). \end{equation*}

Вычисление кубического корня начинается с вычитания из первой левой группы единицы, затем, для получения следующих конечных разностей используем приращения, вычисленные по выражениям (12) и (13). Вычитания производятся с учётом условия “до тех пор, пока остаток не станет меньше следующего вычитаемого”.

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

Теперь, учитывая всё вышесказанное о квадратных и кубических корнях, выведем по тем же правилам аналогичные выражения для вычисления корней 5-й степени. Здесь конечные разности вычисляются с помощью выражения (4), с использованием подстановок x = i и (6)

    \[\Delta(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=\]

(14)   \begin{equation*} =5(10A+i_d)((10A+i_d)((10A+i_d)(10A+i_d+2)+2)+1)+1. \end{equation*}

Вычисление корня 5-й степени принципиально не отличается от вычисления кубического корня. Для примера, вычислим корень 5-й степени из 7167031,46875. Сначала разбиваем число на группы по пять цифр (пятистепенные разряды), вправо и влево от запятой, получим 00071'67031,46875. Далее вычитаем из первой слева группы конечные разности, вычисленные с использованием выражения (14). Условие окончания вычитаний и перехода к обработке следующего разряда такое же, как при извлечении квадратного и кубического корней. После окончания цикла вычитаний к остатку D приписывается следующая группа из пяти цифр, получаем 100000D + C, присваивается новое значение A (A_{new} = 10A + B), и процесс продолжается до достижения необходимой точности. Далее приведём все вычисления, записанные в столбик (пояснения — в фигурных скобках)

Ответ 23,5.

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

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

(15)   \begin{equation*} \Delta\Delta(x^5) = 10(10A + i_d +1)((200A + 40(i_d + 1))A + 2i_d(i_d + 2) + 3). \end{equation*}

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

(16)   \begin{equation*} \Delta(x^5)_{dd}=55(10A+i_d+1)(90A+9i_d+10)(A(10100A+10(202i_d+211))+((101i_d+211)i_d+111)). \end{equation*}

Подведём некоторые итоги. Сразу видно, что с увеличением степени извлекаемого корня коэффициенты многочленов становятся очень большими, структура выражений усложняется, затрудняя выполнение вычислений. Отчасти это связано с тем, что, в целях наглядности, в этой статье все рассуждения и вывод формул проводились, опираясь на привычную всем десятичную систему счисления. Однако, так как машинное вычисление базируется на двоичной системе, приведём, для сравнения, “перевод” только что показанного расчёта корня 5-й степени в двоичный формат.

Основным выражением, определяющим формулу для вычисления конечных разностей с учётом уже вычисленных цифр корня, является (6). Так как при двоичном счёте i_d никогда не будет больше нуля, то для двоичной системы это выражение приобретает вид

    \[i = 10_{(bin)}A = 2_{(dec)}A.\eqno(6b)\]

Подставляя (6b) вместо x в формулу (4) и выполняя алгебраические преобразования с применением схемы Горнера для многочленов, получаем выражение для вычисления конечных разностей функции x^5 в двоичной системе счисления с учётом найденных цифр корня

    \[\Delta(x^5) = (((A + 1)2A + 1)4A + 1)10A + 1.\eqno(14b)\]

Для нашего примера вначале переводим подкоренное число 7167031,46875 в двоичную систему и разбиваем его влево и вправо от запятой на группы по пять цифр (пятистепенные разряды). Для наглядности эти группы будут изображены в виде десятичных эквивалентов в скобках. Вычитаемое вычисляется по (14b). После каждого вычитания при переходе к обработке следующего разряда к остатку D приписывается следующая группа C из пяти двоичных цифр, получаем 100000_{(bin)} \times D + C = 32_{(dec)} \times D + C, вычисляется новое значение A (A_{new} = 10_{(bin)} \times A + B = 2_{(dec)} \times A + B; где B=1, если вычитание выполнялось, или B=0, если не выполнялось), и процесс продолжается до достижения заданной точности.

    \[7167031,46875_{(dec)} = 11011010101110000110111,01111_{(bin)}\]

    \[00110'11010'10111'00001'10111,01111\]

    \[(6)'(26)'(23)'(1)'(23),(15)\]

Ответ: 10111,1_{(bin)}=23,5_{(dec)}.

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

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

    \[\Delta(x^3) = (2A+1)6A+1,\eqno(11b)\]

для квадратного корня

    \[\Delta(x^2) = 4A+1.\eqno(7b)\]

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

Оставьте свой отзыв

Добавить изображение