Методы вычисления квадратных корней

Методы вычисления квадратных корней — это вычислительные алгоритмы для вычисления приближённых значений главных (или неотрицательных) квадратных корней (обычно обозначаемых как ${\displaystyle {\sqrt {S}}}$, ${\displaystyle {\sqrt[{2}]{S}}}$ или ${\displaystyle S^{1/2}}$) вещественного числа. Арифметически это означает, что если дано число ${\displaystyle S}$, процедура находит число, которое при умножении на себя даёт ${\displaystyle S}$. Алгебраически это означает процедуру нахождения неотрицательного корня уравнения ${\displaystyle x^{2}-S=0}$. Геометрически это означает построение стороны квадрата с заданной площадью.

Любое положительное вещественное число имеет два корня^[1]. Главное значение квадратного корня большинства положительных вещественных чисел является иррациональным числом с бесконечной последовательностью десятичных цифр. Как результат, десятичное представление любого такого квадратного корня может быть вычислено только приближённо с конечной точностью (знаков после запятой). Однако, даже если мы берём корень от полного квадрата целого числа, так что результат имеет конечное представление, некоторые процедуры, используемые для вычисления корня, могут вернуть лишь ряд приближений с возрастающей точностью.

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

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

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

Процедуры поиска квадратных корней (в частности, корня из 2) известны по меньшей мере со времён древнего Вавилона (XVII век до нашей эры). Метод Герона из Египта первого века был первым проверяемым алгоритмом для вычисления квадратного корня. Современные аналитические методы начались разрабатываться после принятия арабских цифр в Западной Европе в Раннем Ренессансе. В настоящие дни почти все вычислительные устройства имеют функцию быстрого и точного вычисления квадратного корня в виде встроенной конструкции языка программирования, библиотечной функции или аппаратного оператора, которые основываются на описанных ниже процедурах.

Начальная оценка

Многие итеративные алгоритмы вычисления квадратного корня требуют задания начального случайного значения. Это значение должно быть ненулевым положительным числом. Оно должно быть между 1 и ${\displaystyle S}$, числом, квадратным корень которого мы ищем, поскольку квадратный корень должен быть в этих пределах. Если начальное значение очень далеко от корня, алгоритму потребуется больше итераций. Если начать с ${\displaystyle x_{0}=1}$ (или с ${\displaystyle S}$), будет отработано лишних примерно ${\displaystyle {\tfrac {1}{2}}\vert \log _{2}S\vert }$ итераций просто для получения порядка корня. Поэтому полезно иметь грубую оценку корня, которая может иметь слабую точность, но зато легко вычисляется. В общем случае чем точнее оценка, тем быстрее сходимость. Для метода Ньютона (называемого также вавилонским или методом Герона), начальное значение несколько большее корня даёт более быструю сходимость, по сравнению с начальным значением, меньшим корня.

Вообще говоря, оценка рассматривается на произвольном интервале, в котором известно, что в нём содержится корень (таком как ${\displaystyle [x_{0},S/x_{0}]}$). Получение лучшей оценки вовлекает либо получение более узких границ интервала, либо лучшего функционального приближения к ${\displaystyle f(x).}$ Последнее обычно означает использование для аппроксимации многочленов более высокого порядка, хотя не все аппроксимации используют многочлены. Общие методы оценки бывают скалярные, линейные, гиперболические и логарифмические. Десятичная система счисления обычно используется для оценки в уме или на бумаге. Двоичная система счисления более пригодна для компьютерных оценок. При оценке экспонента и мантисса обычно обрабатываются отдельно.

Десятичная оценка

Обычно число ${\displaystyle S}$ выражается в экспоненциальном виде как ${\displaystyle a\times 10^{2n}}$, где ${\displaystyle 1\leqslant a<100}$, а n — целое число, тогда оценкой возможного квадратного корня может быть ${\displaystyle {\sqrt {a}}\times 10^{n}}$, где ${\displaystyle 1\leqslant {\sqrt {a}}<10}$.

Скалярные оценки

Скалярные методы делят весь диапазон на интервалы и оценка в каждом интервале представлена одним числом. Если диапазон рассматривается как один интервал, то арифметическое среднее (5,5) или геометрическое среднее (${\displaystyle {\sqrt {10}}\approx 3{,}16)\times 10^{n}}$ являются приемлемыми оценками. Абсолютная и относительная оценка для этих оценок будет отличаться. В общем случае отдельное число будет очень неточно. Более точные оценки разбивают диапазон на два и более интервалов, но скалярная оценка продолжает оставаться очень грубой.

Для двух интервалов, разбитых геометрически, квадратный корень ${\displaystyle {\sqrt {S}}={\sqrt {a}}\times 10^{n}}$ можно оценить как^[2].

${\displaystyle {\sqrt {S}}\approx {\begin{cases}2\cdot 10^{n}&{\text{если }}a<10,\\6\cdot 10^{n}&{\text{если }}a\geqslant 10.\end{cases}}}$

Эта оценка имеет максимальную абсолютную погрешность ${\displaystyle 4\cdot 10^{n}}$ в точке = 100 и максимальную относительную ошибку в 100 % в точке = 1.

Например, для ${\displaystyle S=125348}$ с разложением ${\displaystyle 12{,}5348\times 10^{4}}$, оценка будет ${\displaystyle {\sqrt {S}}\approx 6\cdot 10^{2}=600}$. ${\displaystyle {\sqrt {125348}}=354{,}0}$, с абсолютной ошибкой 246 и относительной ошибкой почти 70 %.

Линейная оценка

Лучшей оценкой и стандартным методом является линейное приближение функции ${\displaystyle y=x^{2}}$ на малой дуге. Если, как и выше, степень выделена из числа ${\displaystyle S}$, а интервал сокращён до ${\displaystyle [1,100]}$, можно использовать секущую или касательную где-то вдоль дуги для аппроксимации, но прямая регрессии метода наименьших квадратов будет более точной.

Прямая, получающаяся методом наименьших квадратов, минимизирует среднее расстояние между оценкой и значением функции. Её уравнение — ${\displaystyle y=8{,}7x-10}$. После преобразования ${\displaystyle x=0{,}115y+1{,}15}$ и округления коэффициентов для упрощения вычислений получим

${\displaystyle {\sqrt {S}}\approx (a/10+1{,}2)\cdot 10^{n}}$

Это лучшая оценка в среднем, которую можно получить одной попыткой линейной аппроксимации функции ${\displaystyle y=x^{2}}$ в интервале ${\displaystyle [1,100]}$. Оценка имеет максимальную абсолютную ошибку 1,2 в точке a=100 и максимальную относительную ошибку в 30 % в точках S=1 и 10^[3].

Чтобы разделить на 10, вычитаем единицу из показателя степени ${\displaystyle a}$ или, образно говоря, передвигаем десятичную запятую на одну позицию влево. Для этой формулы любая добавленная константа, равная 1 плюс маленькое приращение, даёт удовлетворительную оценку, так что запоминать точное число нет необходимости. Аппроксимация (округлённая или не округлённая) с помощью одной прямой, стягивающей область ${\displaystyle [1,100]}$ по точности даёт не более одного верного знака. Относительная ошибка более чем 1/2², так что даёт менее 2 битов информации. Точность сильно ограничена, поскольку область охватывает два порядка, что достаточно большая величина для такого рода оценок.

Существенно лучшую оценку можно получить при помощи кусочно-линейной аппроксимации, то есть с помощью нескольких отрезков, которые приближают поддугу исходной дуги. Чем больше отрезков используется, тем лучше приближение. Наиболее употребительно применение касательных. Критичным моментом является как делить дугу и где располагать точки касания. Действенным методом деления дуги от y=1 до y=100 является геометрический — для двух интервалов границей интервалов является квадратный корень исходного интервала, 1*100, то есть ${\displaystyle [1,{\sqrt[{2}]{100}}]}$ и ${\displaystyle [{\sqrt[{2}]{100}},100]}$. Для трёх интервалов будут кубические корни — ${\displaystyle [1,{\sqrt[{3}]{100}}],[{\sqrt[{3}]{100}},({\sqrt[{3}]{100}})^{2}]}$, и ${\displaystyle [({\sqrt[{3}]{100}})^{2},100]}$, и так далее. Для двух интервалов ${\displaystyle {\sqrt[{2}]{100}}=10}$ является очень удобным числом. Легко получить касательные прямые в точках касания ${\displaystyle x={\sqrt {1*{\sqrt {10}}}}}$ и ${\displaystyle x={\sqrt {10*{\sqrt {10}}}}}$. Их уравнения: ${\displaystyle y=3{,}56x-3{,}16}$ и ${\displaystyle y=11{,}2x-31{,}6}$. Обращая уравнения, получим, что квадратные корни равны ${\displaystyle x=0{,}28y+0{,}89}$ и ${\displaystyle x=0{,}089y+2{,}8}$. Тогда для ${\displaystyle S=a\cdot 10^{2n}}$:

${\displaystyle {\sqrt {S}}\approx {\begin{cases}(0{,}28a+0{,}89)\cdot 10^{n}&{\text{если }}a<10,\\(0{,}089a+2{,}8)\cdot 10^{n}&{\text{если }}a\geqslant 10.\end{cases}}}$

Максимальные абсолютные значения оказываются в правых границах интервалов, в точках a=10 и 100, и равны 0,54 и 1,7 соответственно. Максимальные относительные ошибки появляются на концах интервалов, в точках a=1, 10 и 100, и равны 17 %. 17 % или 0,17. Они больше, чем 1/10, так что метод даёт точность менее одной значащей цифры.

Гиперболическая оценка

В некоторых случаях может оказаться действенной гиперболическая оценка, поскольку гипербола также является выпуклой кривой и может лежать вдоль дуги Y = x² лучше, чем прямая. Гиперболическая оценка вычислительно более сложная, поскольку для неё нужно деление на число с плавающей запятой. Почти оптимальной гиперболической аппроксимацией к x² на интервале ${\displaystyle [1,100]}$ является ${\displaystyle y=190/(10-x)-20}$. После преобразования получим ${\displaystyle x=-190/(y+20)+10}$. Тогда для ${\displaystyle S=a\cdot 10^{2n}}$:

${\displaystyle {\sqrt {S}}\approx \left({\frac {-190}{a+20}}+10\right)\cdot 10^{n}}$

Деление с плавающей запятой должно быть с точностью до одного десятичного знака, поскольку вся оценка даёт такую точность, и такое деление можно выполнить в уме. Гиперболическая оценка в среднем лучше, чем скалярная или линейная оценка. Её максимальная абсолютная ошибка составляет 1,58 в точке 100, а максимальная относительная ошибка составляет 16,0 % в точке 10. Для худшего случая a=10 оценка равна 3,67. Если начать с 10 и применять итерации Нютона-Рафсона напрямую, требуется две итерации, которые дают 3,66, прежде чем достичь точности гиперболической оценки. Для более типичного случая наподобие 75 гиперболическая оценка даёт 8,00 и требуется 5 итераций Ньютона-Рафсона с начальным значением 75, чтобы получить более точный результат.

Арифметическая оценка

Метод, аналогичный кусочно-линейной аппроксимации, но использующий лишь арифметические операции вместо алгебраических уравнений, использует таблицу умножения в обратную сторону — квадратный корень чисел между 1 и 100 где-то между 1 и 10, так что, поскольку мы знаем, что 25 является точным квадратом (5 × 5) и 36 является точным квадратом (6 × 6), то квадратный корень из числа, которое больше 25, но меньше 36, начинается с цифры 5. Аналогично для чисел между другими квадратами. Этот метод даёт правильный первый знак, но точность его всего одна цифра — первая цифра квадратного корня из 35, например, равна 5, но сам корень из 35 почти равен 6.

Лучше делить интервал между двумя квадратами пополам. Так что корень любого числа между 25 и половины пути до 36 (что есть 30,5) оценивается как 5, остальные числа, большие 30,5 вплоть до 36 оцениваются как 6^[4]. Процедура требует очень мало арифметики для нахождения середины двух произведений из таблицы. Вот таблица таких чисел:

a	nearest square	${\displaystyle k={\sqrt {a}}}$ est.
1 to 2,5	1 (= 12)	1
2,5 to 6,5	4 (= 22)	2
6,5 to 12,5	9 (= 32)	3
12,5 to 20,5	16 (= 42)	4
20,5 to 30,5	25 (= 52)	5
30,5 to 42,5	36 (= 62)	6
42,5 to 56,5	49 (= 72)	7
56,5 to 72,5	64 (= 82)	8
72,5 to 90,5	81 (= 92)	9
90,5 to 100	100 (= 102)	10

Конечной операцией будет умножение оценки k на степень десятки, делённой пополам, так что для ${\displaystyle S=a\cdot 10^{2n}}$,

${\displaystyle {\sqrt {S}}\approx k\cdot 10^{n}}$

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

Метод можно распространить до 3 значащих цифр в большинстве случаев, интерполируя между ближайшими квадратами. Если ${\displaystyle k^{2}\leqslant a<(k+1)^{2}}$, то ${\displaystyle {\sqrt {a}}}$ примерно равен k плюс дробь, равная разности a и ${\displaystyle k^{2}}$, делённой на разность между двумя квадратами:

${\displaystyle {\sqrt {a}}\approx k+R}$ где ${\displaystyle R={\frac {(a-k^{2})}{(k+1)^{2}-k^{2}}}}$

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

${\displaystyle {\sqrt {S}}={\sqrt {a}}\cdot 10^{n}\approx (k+R)\cdot 10^{n}}$

Число k есть десятичная цифра, а R есть дробь, которую следует превратить в десятичную. Дробь имеет обычно одну цифру в числителе и одну или две цифры в знаменателе, так что преобразование в десятичную дробь можно провести в уме.

Пример: найти квадратный корень из 75. ${\displaystyle 75=75\times 10^{2\cdot 0}}$, так что a равно 75, а n равно 0. Исходя из таблицы умножения квадратный корень мантиссы должен быть 8 с дробью, поскольку ${\displaystyle 8\times 8=64}$, а ${\displaystyle 9\times 9=81}$, слишком велико. Так что k равно 8 с дробью является десятичным представлением R. Дробь R имеет ${\displaystyle 75-k^{2}=11}$ в числителе и ${\displaystyle 81-k^{2}=17}$ в знаменателе. 11/17 чуть меньше, чем 12/18, что равно 2/3 или 0,67, так что 0,66 является хорошим предположением (здесь можно ограничиться и предположением поскольку ошибка мала). Так что оценка корня равна ${\displaystyle 8+0{,}66=8{,}66}$. √75 до трёх значащих цифр будет 8,66, так что оценка до трёх значащих цифр хорошая. Не все оценки с помощью такого метода столь точны, но они довольно близки.

Двоичная оценка

Когда работа ведётся в двоичной системе счисления (скажем, в процессоре компьютера), ${\displaystyle S}$ выражается как ${\displaystyle a\times 2^{2n}}$, где ${\displaystyle 0{,}1_{2}\leqslant a<10_{2}}$, квадратный корень ${\displaystyle {\sqrt {S}}={\sqrt {a}}\times 2^{n}}$ можно оценить величиной

${\displaystyle {\sqrt {S}}\approx (0{,}485+0{,}485\cdot a)\cdot 2^{n}}$

что является регрессией методом наименьших квадратов по 3 старшим битам. ${\displaystyle {\sqrt {a}}}$ имеет максимальную абсолютную ошибку 0,0408 в точке ${\displaystyle a}$=2 и максимальную относительную ошибку в 3,0 % в точке ${\displaystyle a}$=1. Для вычислений удобна округлённая оценка (поскольку коэффициенты являются степенями 2)

${\displaystyle {\sqrt {S}}\approx (0{,}5+0{,}5\cdot a)\cdot 2^{n}}$^[5]

которая имеет максимальную абсолютную ошибку 0,086 в точке 2 и максимальную относительную ошибку в 6,1 % в точках ${\displaystyle a=0{,}5}$ и ${\displaystyle a=2{,}0}$.

Для ${\displaystyle S=125348=1\;1110\;1001\;1010\;0100_{2}=1{,}1110\;1001\;1010\;0100_{2}\times 2^{16}\,}$ двоичное приближение даёт ${\displaystyle {\sqrt {S}}\approx (0{,}5+0{,}5\cdot a)\cdot 2^{8}=1{,}0111\;0100\;1101\;0010_{2}\cdot 1\;0000\;0000_{2}=1{,}456\cdot 256=372{,}8.}$ Поскольку ${\displaystyle {\sqrt {125348}}=354{,}0}$, оценка даёт абсолютную ошибку в 19 и относительную ошибку 5,3 %. Относительная ошибка чуть меньше 1/2⁴, так что приближение даёт точность до 4+ бит.

Оценку для ${\displaystyle a}$ с точностью до 8 бит можно получить путём просмотра таблицы по старшим 8 битам ${\displaystyle a}$, учитывая, что старший бит задаётся неявно в большинстве представлений чисел с плавающей запятой, а младшие биты после 8 бит должны быть округлены. Таблица содержит 256 байт заранее вычисленных 8-битных квадратных корней. Например, для индекса 11101101₂, что в десятичной системе равно 1,8515625₁₀, значение в таблице равно 10101110₂, что в десятичной системе равно 1,359375₁₀, квадратному корню числа 1,8515625₁₀ с точностью до 8 бит (2+ десятичных знака).

Вавилонский метод

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

Возможно первым алгоритмом, используемым для аппроксимации ${\displaystyle {\sqrt {S}}}$, является метод, известный как вавилонский метод, несмотря на то, что нет никаких прямых свидетельств, за исключением гипотетических умозаключений, что вавилонские математики использовали этот метод^[6]. Метод известен также как метод Герона, по имени греческого математика первого столетия Герона, который дал первое явное описание метода в своей работе 60 года Метрика^[7].Основная методика заключается в том, что если x больше квадратного корня неотрицательного вещественного числа S то ${\displaystyle {\tfrac {S}{x}}}$ будет меньше корня и наоборот. Так что среднее этих двух чисел резонно ожидать более близким к корню (формальное доказательство этого факта основывается на неравенстве о среднем арифметическом, геометрическом и гармоническом, которое показывает, что это среднее всегда больше квадратного корня, что обеспечивает сходимость). Метод эквивалентен использованию метода Ньютона для решения уравнения ${\displaystyle x^{2}-S=0}$.

Точнее, если x является начальным приближением для ${\displaystyle {\sqrt {S}}}$, а ${\displaystyle \varepsilon }$ ошибка в нашей оценке, такая что ${\displaystyle S=(x+\varepsilon )^{2}}$, мы можем раскрыть скобки и получим

${\displaystyle \varepsilon ={\frac {S-x^{2}}{2x+\varepsilon }}\approx {\frac {S-x^{2}}{2x}},}$ поскольку ${\displaystyle \varepsilon \ll x}$.

Следовательно, мы можем компенсировать ошибку и обновить нашу старую оценку

${\displaystyle x+\varepsilon \approx x+{\frac {S-x^{2}}{2x}}={\frac {S+x^{2}}{2x}}={\frac {{\frac {S}{x}}+x}{2}}\equiv x_{\text{обновлённый}}}$

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

1. Начинаем с любого положительного начального значения ${\displaystyle x_{0}}$ (чем ближе к истинному квадратному корню числа S, тем лучше).
2. Положим ${\displaystyle x_{n+1}}$ равным среднему между ${\displaystyle x_{n}}$ и ${\displaystyle {\tfrac {S}{x_{n}}}}$ (используем среднее арифметическое для аппроксимации среднего геометрического).
3. Повторяем шаг 2 пока не достигнем нужной точности.

Алгоритм можно представить следующим образом:

${\displaystyle x_{0}\approx {\sqrt {S}},}$

${\displaystyle x_{n+1}={\frac {1}{2}}\left(x_{n}+{\frac {S}{x_{n}}}\right),}$

${\displaystyle {\sqrt {S}}=\lim _{n\to \infty }x_{n}.}$

Алгоритм работает также хорошо и для p-адических чисел, но не может быть использован для отождествления вещественных квадратных корней с p-адичными квадратными корнями. Можно, например, построить последовательность рациональных чисел, полученных этим методом, которая сходится к +3 в случае вещественных чисел, но к -3 в 2-адичных числах.

Пример

Для вычисления √S, где S = 125348, с точностью до шести значащих цифр используем метод грубой оценки, описанный выше

${\displaystyle {\begin{aligned}{\begin{array}{rlll}x_{0}&=6\cdot 10^{2}&&=600{,}000\\[0,3em]x_{1}&={\frac {1}{2}}\left(x_{0}+{\frac {S}{x_{0}}}\right)&={\frac {1}{2}}\left(600{,}000+{\frac {125348}{600{,}000}}\right)&=404{,}457\\[0.3em]x_{2}&={\frac {1}{2}}\left(x_{1}+{\frac {S}{x_{1}}}\right)&={\frac {1}{2}}\left(404{,}457+{\frac {125348}{404{,}457}}\right)&=357{,}187\\[0.3em]x_{3}&={\frac {1}{2}}\left(x_{2}+{\frac {S}{x_{2}}}\right)&={\frac {1}{2}}\left(357{,}187+{\frac {125348}{357{,}187}}\right)&=354{,}059\\[0.3em]x_{4}&={\frac {1}{2}}\left(x_{3}+{\frac {S}{x_{3}}}\right)&={\frac {1}{2}}\left(354{,}059+{\frac {125348}{354{,}059}}\right)&=354{,}045\\[0.3em]x_{5}&={\frac {1}{2}}\left(x_{4}+{\frac {S}{x_{4}}}\right)&={\frac {1}{2}}\left(354{,}045+{\frac {125348}{354{,}045}}\right)&=354{,}045\end{array}}\end{aligned}}}$

Поэтому ${\displaystyle {\sqrt {125348}}\approx 354{,}045}$.

Сходимость

Предположим, что x₀ > 0 и S > 0. Тогда для любого n x_n > 0. Относительная ошибка^[англ.] x_n определена как

${\displaystyle \varepsilon _{n}={\frac {x_{n}}{\sqrt {S}}}-1>-1}$

а тогда

${\displaystyle x_{n}={\sqrt {S}}\cdot (1+\varepsilon _{n}).}$

Теперь можно показать, что

${\displaystyle \varepsilon _{n+1}={\frac {\varepsilon _{n}^{2}}{2(1+\varepsilon _{n})}}\geqslant 0.}$

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

${\displaystyle \varepsilon _{n+2}\leqslant \min \left\{{\frac {\varepsilon _{n+1}^{2}}{2}},{\frac {\varepsilon _{n+1}}{2}}\right\}}$

а отсюда следует гарантированная сходимость и эта сходимость квадратичная.

Сходимость в худшем случае

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

${\displaystyle {\begin{aligned}S&=1;&x_{0}&=2;&x_{1}&=1{,}250;&\varepsilon _{1}&=0{,}250.\\S&=10;&x_{0}&=2;&x_{1}&=3{,}500;&\varepsilon _{1}&<0{,}107.\\S&=10;&x_{0}&=6;&x_{1}&=3{,}833;&\varepsilon _{1}&<0{,}213.\\S&=100;&x_{0}&=6;&x_{1}&=11{,}333;&\varepsilon _{1}&<0{,}134.\end{aligned}}}$

А тогда в любом случае

${\displaystyle \varepsilon _{1}\leqslant 2^{-2}.\,}$

${\displaystyle \varepsilon _{2}<2^{-5}<10^{-1}.\,}$

${\displaystyle \varepsilon _{3}<2^{-11}<10^{-3}.\,}$

${\displaystyle \varepsilon _{4}<2^{-23}<10^{-6}.\,}$

${\displaystyle \varepsilon _{5}<2^{-47}<10^{-14}.\,}$

${\displaystyle \varepsilon _{6}<2^{-95}<10^{-28}.\,}$

${\displaystyle \varepsilon _{7}<2^{-191}<10^{-57}.\,}$

${\displaystyle \varepsilon _{8}<2^{-383}<10^{-115}.\,}$

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

Метод Бакхшали

Этот метод для поиска приближения квадратного корня был написан в древнеиндийской рукописи, называемой манускриптом Бакхшали. Метод эквивалентен двум итерациям вавилонского метода с начальным значением x₀. Таким образом, алгоритм является квадратично сходящимся, что означает, что число верных знаков приближения увеличивается примерно в четыре раза с каждой итерацией^[8]. Представление алгоритма в современной нотации следующее: Следует вычислить ${\displaystyle {\sqrt {S}}}$, пусть x₀² будет начальным приближением к корню S. Последовательно выполняются итерации

${\displaystyle {\begin{aligned}a_{n}&={\frac {S-x_{n}^{2}}{2x_{n}}},\\b_{n}&=x_{n}+a_{n},\\x_{n+1}&=b_{n}-{\frac {a_{n}^{2}}{2b_{n}}}=(x_{n}+a_{n})-{\frac {a_{n}^{2}}{2(x_{n}+a_{n})}}.\end{aligned}}}$

Это можно использовать для построения рационального приближения к квадратному корню, начав с целого числа. Если ${\displaystyle x_{0}=N}$ — это целое число, выбранное так, что ${\displaystyle N^{2}}$ близко к S, и ${\displaystyle d=S-N^{2}}$ — это разность, абсолютная величина которой минимизируется, то первую итерацию можно записать следующим образом:

${\displaystyle {\sqrt {S}}\approx N+{\frac {d}{2N}}-{\frac {d^{2}}{8N^{3}+4Nd}}={\frac {8N^{4}+8N^{2}d+d^{2}}{8N^{3}+4Nd}}={\frac {N^{4}+6N^{2}S+S^{2}}{4N^{3}+4NS}}={\frac {N^{2}(N^{2}+6S)+S^{2}}{4N(N^{2}+S)}}.}$

Метод Бакхшали может быть обобщён для вычисления произвольного корня, включая дробные корни^[9].

Пример

Используем тот же пример, что был приведён для вавилонского метода. Пусть ${\displaystyle S=125348.}$ Тогда первая итерация даёт

${\displaystyle {\begin{aligned}x_{0}&=600\\a_{1}&={\frac {125348-600^{2}}{2\times 600}}&&=&-195{,}543\\b_{1}&=600+(-195{,}543)&&=&404{,}456\\x_{1}&=404{,}456-{\frac {(-195{,}543)^{2}}{2\times 404{,}456}}&&=&357{,}186\end{aligned}}}$

Аналогично вторая итерация даёт

${\displaystyle {\begin{aligned}a_{2}&={\frac {125348-357{,}186^{2}}{2\times 357{,}186}}&&=&-3{,}126\\b_{2}&=357{,}186+(-3{,}126)&&=&354{,}060\\x_{2}&=354{,}06-{\frac {(-3{,}1269)^{2}}{2\times 354{,}06}}&&=&354{,}046\end{aligned}}}$

Цифра за цифрой

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

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

Палочки Непера включают дополнительные средства для выполнения этого алгоритма. Алгоритм вычисления n-го корня цифра за цифрой^[англ.] является обобщением этого метода.

Основной принцип

Рассмотрим сначала случай нахождения квадратного корня из числа Z, являющегося квадратом двузначного числа XY, где X — это цифра десятков, а Y — цифра единиц. Имеем:

${\displaystyle Z=(10X+Y)^{2}=100X^{2}+20XY+Y^{2}}$

Сначала определим значение X. X — это наибольшая цифра, такая что X² не превосходит Z, от которого отброшены две последние цифры.

На следующей итерации соединяем пару цифр, умножая X на 2 и помещая результат в позицию десятков, а затем пытаемся найти, чему же равно Y.

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

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

${\displaystyle N=(a_{1}+a_{2}+a_{3}+\dotsb +a_{n})^{2}.}$

Путём многократного использования тождества

${\displaystyle (x+y)^{2}=x^{2}+2xy+y^{2},}$

правую часть можно представить в виде

${\displaystyle {\begin{aligned}&(a_{1}+a_{2}+a_{3}+\dotsb +a_{n})^{2}\\=&\,a_{1}^{2}+2a_{1}a_{2}+a_{2}^{2}+2(a_{1}+a_{2})a_{3}+a_{3}^{2}+\dotsb +a_{n-1}^{2}+2\left(\sum _{i=1}^{n-1}a_{i}\right)a_{n}+a_{n}^{2}\\=&\,a_{1}^{2}+[2a_{1}+a_{2}]a_{2}+[2(a_{1}+a_{2})+a_{3}]a_{3}+\dotsb +\left[2\left(\sum _{i=1}^{n-1}a_{i}\right)+a_{n}\right]a_{n}.\end{aligned}}}$

Это выражение позволяет нам найти квадратный корень последовательным подбором значений ${\displaystyle a_{i}}$. Предположим, что числа ${\displaystyle a_{1},\ldots ,a_{m-1}}$ уже подобраны, тогда m-й член задаётся выражением ${\displaystyle Y_{m}=[2P_{m-1}+a_{m}]a_{m},}$, где ${\displaystyle P_{m-1}=\sum _{i=1}^{m-1}a_{i}}$ является найденным приближением к квадратному корню. Теперь каждый новый подбор ${\displaystyle a_{m}}$ должен удовлетворять рекурсии

${\displaystyle X_{m}=X_{m-1}-Y_{m},}$

так что ${\displaystyle X_{m}\geqslant 0}$ для всех ${\displaystyle 1\leqslant m\leqslant n,}$ при начальном значении ${\displaystyle X_{0}=N.}$ Если ${\displaystyle X_{n}=0,}$ найден точный квадратный корень. Если нет, то сумма ${\displaystyle a_{i}}$ даёт подходящую аппроксимацию к квадратному корню и ${\displaystyle X_{n}}$ будет ошибкой аппроксимации.

Например, в десятичной системе мы имеем

${\displaystyle N=(a_{1}\cdot 10^{n-1}+a_{2}\cdot 10^{n-2}+\cdots +a_{n-1}\cdot 10+a_{n})^{2},}$

где ${\displaystyle 10^{n-i}}$ являются указателями положения цифр, а коэффициенты ${\displaystyle a_{i}\in \{0,1,2,\ldots ,9\}}$. На каждом m-м шаге вычисления квадратного корня находится приближённый квадратный корень. Величина ${\displaystyle P_{m-1}}$ и суммируемые члены ${\displaystyle Y_{m}}$ задаются формулами

${\displaystyle P_{m-1}=\sum _{i=1}^{m-1}a_{i}\cdot 10^{n-i}=10^{n-m+1}\sum _{i=1}^{m-1}a_{i}\cdot 10^{m-i-1},}$

${\displaystyle Y_{m}=[2P_{m-1}+a_{m}\cdot 10^{n-m}]a_{m}\cdot 10^{n-m}=\left[20\sum _{i=1}^{m-1}a_{i}\cdot 10^{m-i-1}+a_{m}\right]a_{m}\cdot 10^{2(n-m)}.}$

Поскольку указатели положения ${\displaystyle Y_{m}}$ имеют чётную степень 10, нам нужно работать только с парой старших цифр в оставшемся члене ${\displaystyle X_{m-1}}$ на любом m-м шаге. Раздел ниже систематизирует эту процедуру.

Очевидно, что подобный метод может быть использован для вычисления квадратного корня в любой системе счисления, не обязательно в десятичной. Например, нахождение цифра за цифрой квадратного корня в двоичной системе довольно эффективно, поскольку значение ${\displaystyle a_{i}}$ ищется в малом наборе цифр {0,1}. Это делает вычисление более быстрым, поскольку на каждом шаге значение ${\displaystyle Y_{m}}$ либо равно ${\displaystyle Y_{m}=0}$ для ${\displaystyle a_{m}=0}$, либо ${\displaystyle Y_{m}=2P_{m-1}+1}$ для ${\displaystyle a_{m}=1}$. Факт, что имеется всего две возможности для ${\displaystyle a_{m}}$ также делает проще процесс выбора значения ${\displaystyle a_{m}}$ на m-м шаге вычислений. Это потому, что нам нужно лишь проверить, что ${\displaystyle Y_{m}\leqslant X_{m-1}}$ для ${\displaystyle a_{m}=1.}$ Если это условие выполняется, мы берём ${\displaystyle a_{m}=1}$; а если не выполняется, то берём ${\displaystyle a_{m}=0.}$ Также факт, что умножение на 2 осуществляется сдвигом влево, помогает при вычислениях.

Десятичная система счисления

Запишем исходное число в десятичном виде. Числа, записываются по аналогии алгоритму деления столбиком, и, как и в длинном делении, квадратный корень будет писаться в верхней строке. Теперь разобьём цифры на пары, начиная с запятой, в обе стороны от неё. Десятичная запятая квадратного корня будет на десятичной запятой квадрата. Одна цифра квадратного корня записывается над парой цифр квадрата.

Начиная с крайне левой позиции выполняем следующую процедуру для каждой пары цифр:

1. Сносим вниз старшую пару ещё неиспользованных цифр (если все цифры использованы, пишем "00") и записываем их справа от остатка предыдущего шага (на первом шаге остатка нет). Другими словами, умножаем остаток на 100 и добавляем две цифры. Это будет текущим значением c.
2. Находим p, y и x следующим образом:
   • Пусть ${\displaystyle {\boldsymbol {p}}}$ будет частью уже найденного корня при игнорировании десятичной запятой. (Для первого шага, p = 0.)
   • Определяем наибольшую цифру ${\displaystyle {\boldsymbol {x}}}$, такую что ${\displaystyle x(20p+x)\leqslant c}$. Мы будем использовать новую переменную ${\displaystyle {\boldsymbol {y}}=x(20p+x)}$.
     • Заметим, что ${\displaystyle 20p+x}$ — это просто удвоенное p с дописанной к нему справа цифрой x.
     • Заметим, что x можно найти случайным угадыванием, чему должно равняться c/(20•p), с пробным вычислением y, затем x берётся выше или ниже в зависимости от результата.
   • Помещаем цифру ${\displaystyle x}$ в качестве следующей цифры корня, то есть над парой цифр квадрата, которые только что спустили вниз. Теперь следующее значение p получается из старого значения по формуле ${\displaystyle p\times 10+x}$.
3. Вычитаем y из c для образования нового остатка.
4. Если остаток равен нулю и нет больше цифр, которые можно спустить вниз, алгоритм останавливается. В противном случае возвращаемся на шаг 1 и выполняем следующую итерацию.

Примеры

Находим квадратный корень из 152,2756.

          1  2. 3  4 
       /
     \/  01 52,27 56

         01                   1*1 <= 1 < 2*2                 x = 1
         01                     y = x*x = 1*1 = 1
         00 52                22*2 <= 52 < 23*3              x = 2
         00 44                  y = (20+x)*x = 22*2 = 44
            08 27             243*3 <= 827 < 244*4           x = 3
            07 29               y = (240+x)*x = 243*3 = 729
               98 56          2464*4 <= 9856 < 2465*5        x = 4
               98 56            y = (2460+x)*x = 2464*4 = 9856
               00 00          Алгоритм останавливается: Ответ 12,34

Двоичная система счисления

Этот раздел использует формализм раздела «Вычисление цифра за цифрой» с небольшими изменениями, что ${\displaystyle N^{2}=(a_{n}+\dotsb +a_{0})^{2}}$, а каждое ${\displaystyle a_{m}}$ равно ${\displaystyle 2^{m}}$ или ${\displaystyle 0}$.
Теперь мы пробегаем по всем ${\displaystyle 2^{m}}$ от ${\displaystyle 2^{n}}$ вниз до ${\displaystyle 2^{0}}$ и строим приближённое решение ${\displaystyle P_{m}=a_{n}+a_{n-1}+\ldots +a_{m}}$ в виде суммы всех ${\displaystyle a_{i}}$, для которых мы найдём значение.
Чтобы определить, равно ли ${\displaystyle a_{m}}$ значению ${\displaystyle 2^{m}}$ или ${\displaystyle 0}$, мы берём ${\displaystyle P_{m}=P_{m+1}+2^{m}}$. Если ${\displaystyle P_{m}^{2}\leqslant N^{2}}$ (то есть квадрат нашего приближения включая ${\displaystyle 2^{m}}$ не превосходит исходного квадрата), то полагаем ${\displaystyle a_{m}=2^{m}}$, в противном случае полагаем ${\displaystyle a_{m}=0}$ и ${\displaystyle P_{m}=P_{m+1}}$.
Чтобы избежать возведения в квадрат ${\displaystyle P_{m}}$ на каждом шаге, мы запоминаем разность ${\displaystyle X_{m}=N^{2}-P_{m}^{2}}$ и обновляем её на каждой итерации, полагая ${\displaystyle X_{m}=X_{m+1}-Y_{m}}$ с ${\displaystyle Y_{m}=P_{m}-P_{m+1}=2P_{m+1}a_{m}+a_{m}^{2}}$.
Первоначально мы устанавливаем ${\displaystyle a_{n}=P_{n}=2^{n}}$ для наибольшего ${\displaystyle n}$ с ${\displaystyle (2^{n})^{2}=4^{n}\leqslant N^{2}}$.

В качестве дополнительной оптимизации сохраняем ${\displaystyle P_{m+1}2^{m+1}}$ и ${\displaystyle (2^{m})^{2}}$, два члена ${\displaystyle Y_{m}}$ в случае, когда ${\displaystyle a_{m}}$ не нуль, в отдельных переменных ${\displaystyle c_{m}}$, ${\displaystyle d_{m}}$:

${\displaystyle c_{m}=P_{m+1}2^{m+1}}$

${\displaystyle d_{m}=(2^{m})^{2}}$

${\displaystyle Y_{m}={\begin{cases}c_{m}+d_{m}&{\text{если }}a_{m}=2^{m}\\0&{\text{если }}a_{m}=0\end{cases}}}$

${\displaystyle c_{m}}$ и ${\displaystyle d_{m}}$ можно эффективно обновлять на каждом шаге:

${\displaystyle c_{m-1}=P_{m}2^{m}=(P_{m+1}+a_{m})2^{m}=P_{m+1}2^{m}+a_{m}2^{m}={\begin{cases}c_{m}/2+d_{m}&{\text{если }}a_{m}=2^{m}\\c_{m}/2&{\text{если }}a_{m}=0\end{cases}}}$

${\displaystyle d_{m-1}={\frac {d_{m}}{4}}}$

Заметим, что

${\displaystyle c_{-1}=P_{0}2^{0}=P_{0}=N}$, что является конечным результатом, возвращаемым функцией, представленной ниже.

Реализация алгоритма на языке C^[10]:

int32_t isqrt(int32_t n) 
{ assert(("входное значение должно быть неотрицательным", n > 0));
  int32_t x = n;    // ${\displaystyle X_{n+1}}$
  int32_t c = 0;    // ${\displaystyle c_{n}}$
  // ${\displaystyle d_{n}}$ начинается с наибольшей степени четырёх <= n
  int32_t d = 1 << 30; // Второй старший бит устанавливаем в 1.
                        // То же самое, что ((unsigned)INT32_MAX + 1) / 2.
  while (d > n) d >>= 2;
  while (d != 0)    // для ${\displaystyle d_{n}\dots d_{0}}$
  { if (x >= c + d) // если ${\displaystyle X_{m+1}\geqslant Y_{m}}$, то ${\displaystyle a_{m}=2^{m}}$
    { x -= c + d;       // ${\displaystyle X_{m}=X_{m+1}-Y_{m}}$
     c = (c >> 1) + d;  // ${\displaystyle c{m-1}=c_{m}/2+d_{m}(a_{=}2^{m})}$
   } else
          c >>= 1;      // ${\displaystyle c_{m-1}=c_{m}/2(a_{m}=0)}$
    d >>= 2;            // ${\displaystyle d_{m-1}=d_{m}/4}$
  }
  return c;             // ${\displaystyle c_{-1}}$
}

Можно реализовать более быстрый алгоритм как в двоичной, так и в десятичной системе счисления, если использовать таблицы для выбора, то есть реализация принципа использование больше памяти сокращает время исполнения^[11].

Экспоненциальное тождество

Карманные калькуляторы обычно реализуют хорошие программы вычисления экспоненты и натурального логарифма. Вычисление квадратного корня S тогда производится с помощью свойств логарифмов (${\displaystyle \ln x^{n}=n\ln x}$) и экспоненты (${\displaystyle e^{\ln x}=x}$):

${\displaystyle {\sqrt {x}}=e^{{\frac {1}{2}}\ln x},\,x>0.}$

Или в более общем случае:

${\displaystyle {\sqrt[{n}]{x}}=e^{{\frac {1}{n}}\ln x},\,x>0.}$

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

Такой метод вычисления квадратного корня удобен для калькуляторов, поскольку они обычно не критичны ко времени выполнения операции. Однако ресурсоемкость данного метода делает его малопригодным для использования в ЭВМ, где простые арифметические операции должны обладать минимальными задержками. Тем не менее описанный метод вычисления квадратного корня применялся в ЭВМ ZX Spectrum.

Итеративный метод с двумя переменными

Этот метод применим для поиска квадратного корня из ${\displaystyle 0<S<3\,\!}$ и лучше всего сходится для ${\displaystyle S\approx 1}$. Это, однако, не является существенным ограничением для вычислений на компьютерах, поскольку в представлениях двоичных чисел с плавающей запятой и с фиксированной запятой тривиально умножить ${\displaystyle S\,\!}$ на целую степень числа 4, с последующей коррекцией ${\displaystyle {\sqrt {S}}}$ на нужную степень 2 путём изменения экспоненты или сдвигом соответственно. Таким образом, ${\displaystyle S\,\!}$ может быть сдвинуто в пределы ${\displaystyle {\frac {1}{2}}\leqslant S<2}$. Более того, приведённый ниже метод не использует делений общего вида, а только сложение, вычитание, умножение и деление на степень двойки. Последнее из этих действий тривиально реализуется. Недостатком метода является накопление ошибки, в отличие от итеративных методов с одной переменной, таких как вавилонский.

Начальный шаг метода

${\displaystyle a_{0}=S\,\!}$

${\displaystyle c_{0}=S-1\,\!}$

Итерационные шаги

${\displaystyle a_{n+1}=a_{n}-a_{n}c_{n}/2\,\!}$

${\displaystyle c_{n+1}=c_{n}^{2}(c_{n}-3)/4\,\!}$

Тогда ${\displaystyle a_{n}\rightarrow {\sqrt {S}}}$ (при ${\displaystyle c_{n}\rightarrow 0}$).

Заметим, что сходимость ${\displaystyle c_{n}\,\!}$, а потому и ${\displaystyle a_{n}\,\!}$, квадратична.

Доказательство метода достаточно простое. Сначала перепишем итерационное определение ${\displaystyle c_{n}\,\!}$ как

${\displaystyle 1+c_{n+1}=(1+c_{n})(1-c_{n}/2)^{2}\,\!}$.

Теперь «в лоб» доказывается, что

${\displaystyle S(1+c_{n})=a_{n}^{2}}$

а потому сходимость ${\displaystyle a_{n}\,\!}$ к желаемому результату ${\displaystyle {\sqrt {S}}}$ обеспечивается сходимостью ${\displaystyle c_{n}\,\!}$ к 0, что, в свою очередь, вытекает из ${\displaystyle -1<c_{0}<2\,\!}$.

Этот метод разработали около 1950 года М. В. Уилкс, Д. Дж. Уилер и С. Гилл^[12] для использования в EDSAC, одном из первых электронных компьютеров^[13]. Позднее метод был обобщён на неквадратные корни^[14].

Итеративные методы вычисления обратного к квадратному корню числа

Далее приведены итеративные методы вычисления обратного к квадратному корню из S числа, то есть ${\displaystyle 1/{\sqrt {S}}}$. Если такое значение найдено, находим ${\displaystyle {\sqrt {S}}}$ просто умножением: ${\displaystyle {\sqrt {S}}=S\cdot (1/{\sqrt {S}})}$. Эти итерации используют только умножение и не используют деления. Потому методы быстрее, чем вавилонский метод. Однако методы нестабильны, если начальное значение не близко к обратному к корню значению, итерации расходятся. Поэтому может быть выгодным сначала сделать итерацию вавилонским методом для грубой оценки корня перед началом использования этих методов.

• Применение метода Ньютона для уравнения ${\displaystyle (1/x^{2})-S=0}$ даёт метод, который сходится квадратично и использует три умножения на каждом шаге:

  ${\displaystyle x_{n+1}={\frac {x_{n}}{2}}\cdot (3-S\cdot x_{n}^{2})=x_{n}\cdot \left({\frac {3}{2}}-{\frac {S}{2}}\cdot x_{n}^{2}\right).}$

• Другая итерация получается из метода Галлея^[англ.], которая является методом Хаусхолдера^[англ.] второго порядка. Метод сходится кубически, но использует пять умножений на итерацию:

  ${\displaystyle y_{n}=S\cdot x_{n}^{2}}$

  ${\displaystyle x_{n+1}={\frac {x_{n}}{8}}\cdot (15-y_{n}\cdot (10-3\cdot y_{n}))=x_{n}\cdot \left({\frac {15}{8}}-y_{n}\cdot \left({\frac {10}{8}}-{\frac {3}{8}}\cdot y_{n}\right)\right)}$.

• Если работаем с числами с фиксированной запятой, умножение на 3 и деление на 8 можно реализовать сдвигами и сложениями. При работе с плавающей запятой метод Галлея можно свести к четырём умножениям на итерацию путём предварительного вычисления ${\displaystyle {\sqrt {\frac {3}{8}}}S}$ и поправки всех других констант для компенсации:

  ${\displaystyle y_{n}={\sqrt {\frac {3}{8}}}S\cdot x_{n}^{2}}$,

  ${\displaystyle x_{n+1}=x_{n}\cdot \left({\frac {15}{8}}-y_{n}\cdot \left({\sqrt {\frac {25}{6}}}-y_{n}\right)\right)}$.

Алгоритм Гольдшмидта

Некоторые компьютеры используют алгоритм Гольдшмидта для одновременного вычисления ${\displaystyle {\sqrt {S}}}$ и ${\displaystyle 1/{\sqrt {S}}}$. Алгоритм Гольдшмидта находит ${\displaystyle {\sqrt {S}}}$ быстрее, чем итерация Ньютона-Рапсона, на компьютерах с операциями совмещённого умножения-сложения и имеющих либо конвейерный процессор плавающей запятой, либо два независимых процессора плавающей запятой^[15].

Первый способ записи алгоритма Гольдшмидта начинается с

${\displaystyle b_{0}=S}$

${\displaystyle Y_{0}\approx 1/{\sqrt {S}}}$ (обычно используется поиск в таблице)

${\displaystyle y_{0}=Y_{0}}$

${\displaystyle x_{0}=Sy_{0}}$

и осуществляются итерации

${\displaystyle b_{n+1}=b_{n}Y_{n}^{2}}$

${\displaystyle Y_{n+1}=(3-b_{n+1})/2}$

${\displaystyle x_{n+1}=x_{n}Y_{n+1}}$

${\displaystyle y_{n+1}=y_{n}Y_{n+1}}$

пока ${\displaystyle b_{i}}$ не окажется достаточно близко к 1 или не будет проведено фиксированное число итераций. Итерации сходятся к

${\displaystyle \lim _{n\to \infty }x_{n}={\sqrt {S}}}$,

${\displaystyle \lim _{n\to \infty }y_{n}=1/{\sqrt {S}}}$.

Заметим, что можно опустить вычисление ${\displaystyle x_{n}}$ или ${\displaystyle y_{n}}$, а если оба значения желательны, то ${\displaystyle x_{n}=Sy_{n}}$ можно использовать в конце вместо вычисления на каждой итерации.

Второй способ, использующий операции совмещённого умножения-сложения начинается с

${\displaystyle y_{0}\approx 1/{\sqrt {S}}}$ (обычно используется поиск в таблице)

${\displaystyle x_{0}=Sy_{0}}$

${\displaystyle h_{0}=y_{0}/2}$

и осуществляются итерации

${\displaystyle r_{n}=0{,}5-x_{n}h_{n}}$

${\displaystyle x_{n+1}=x_{n}+x_{n}r_{n}}$

${\displaystyle h_{n+1}=h_{n}+h_{n}r_{n}}$

пока ${\displaystyle r_{i}}$ не станет достаточно близко к 0, либо не будет осуществлено фиксированное число итераций. Значения сходятся к

${\displaystyle \lim _{n\to \infty }x_{n}={\sqrt {S}}}$

${\displaystyle \lim _{n\to \infty }2h_{n}=1/{\sqrt {S}}}$.

Ряды Тейлора

Если N является приближением к ${\displaystyle {\sqrt {S}}}$, лучшее приближения может быть найдено использованием ряда Тейлора функции квадратного корня:

${\displaystyle {\sqrt {N^{2}+d}}=N\sum _{n=0}^{\infty }{\frac {(-1)^{n}(2n)!}{(1-2n)n!^{2}4^{n}}}{\frac {d^{n}}{N^{2n}}}=N\left(1+{\frac {d}{2N^{2}}}-{\frac {d^{2}}{8N^{4}}}+{\frac {d^{3}}{16N^{6}}}-{\frac {5d^{4}}{128N^{8}}}+\cdots \right)}$

Порядок сходимости равен числу используемых членов ряда. При использовании двух членов метод эквивалентен вавилонскому методу. При использовании трёх членов каждая итерация использует почти столько же операций, сколько использует приближение Бакхшали, но сходимость слабее. Поэтому этот метод не является особенно эффективным способом вычисления. Для максимизации скорости сходимости, следует выбрать N так, чтобы ${\displaystyle {\frac {|d|}{N^{2}}}\,}$ было как можно меньше.

Разложение в цепную дробь

См. также: Решение квадратных уравнений с помощью цепных дробей

Квадратичные иррациональности (числа вида ${\displaystyle {\frac {a+{\sqrt {b}}}{c}}}$, где a, b и c целые числа), и, в частности, квадратные корни из целых чисел, имеют периодические цепные дроби^[англ.]. Иногда целью является не нахождение численного значения квадратного корня, а его разложение в цепную дробь, а следовательно его рационального приближения. Пусть S будет положительным числом, корень из которого требуется найти. Теперь пусть a будет начальным приближением, а r будет остаточным членом, тогда мы можем записать ${\displaystyle S=a^{2}+r.}$ Поскольку мы имеем ${\displaystyle S-a^{2}=({\sqrt {S}}+a)({\sqrt {S}}-a)=r}$, мы можем выразить квадратный корень из S как

${\displaystyle {\sqrt {S}}=a+{\frac {r}{a+{\sqrt {S}}}}.}$

Применяя это выражение для ${\displaystyle {\sqrt {S}}}$ к знаменателю дроби, получим

${\displaystyle {\sqrt {S}}=a+{\frac {r}{a+(a+{\frac {r}{a+{\sqrt {S}}}})}}=a+{\frac {r}{2a+{\frac {r}{a+{\sqrt {S}}}}}}.}$

Компактная запись

Числитель/знаменатель разложения для непрерывных дробей (см. слева) затруднительно записывать, а также трудно укладывается в существующую систему форматирования документов. По этой причине была разработана специальная нотация для компактного представления целой и периодической частей непрерывных дробей. Одно из таких соглашений использует лексическую «ломаную линию» для представления черты между числителем и знаменателем, что позволяет записывать дробь горизонтально, а не вертикально:

${\displaystyle {\sqrt {S}}=a+{\frac {r|}{|2a}}+{\frac {r|}{|2a}}+{\frac {r|}{|2a}}+\cdots }$

Здесь каждая горизонтальная черта (в дроби) представлена тремя чертами — двумя вертикальными и одной горизонтальной, которые отделяют ${\displaystyle r}$ от ${\displaystyle 2a}$.

Ещё более компактная нотация имеет специальный вид

${\displaystyle [a;2a,2a,2a,...]}$

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

${\displaystyle [a;{\overline {2a}}]}$

Для √2 значение ${\displaystyle a}$ равно 1, так что представлением будет

${\displaystyle [1;{\overline {2}}]}$

Следуя этим путём мы получаем обобщённую непрерывную дробь^[англ.] для квадратного корня ${\displaystyle {\sqrt {S}}=a+{\cfrac {r}{2a+{\cfrac {r}{2a+{\cfrac {r}{2a+\ddots }}}}}}}$

Первым шагом вычисления такой дроби для получения квадратного корня является подстановки для корня и выбор числа знаменателей. Например, в канонической форме ${\displaystyle r}$ равен 1 и для √2, ${\displaystyle a}$ равен 1, так что численно непрерывной дробью для 3 знаменателей будет

${\displaystyle {\sqrt {2}}\approx 1+{\cfrac {1}{2+{\cfrac {1}{2+{\cfrac {1}{2}}}}}}}$

Шаг 2. Непрерывная дробь свёртывается снизу вверх, один знаменатель за раз, чтобы получить рациональную дробь, числитель и знаменатель которой являются целыми числами. Процесс свёртывания тогда выглядит следующим образом (беря первые три знаменателя):

${\displaystyle 1+{\cfrac {1}{2+{\cfrac {1}{2+{\cfrac {1}{2}}}}}}=1+{\cfrac {1}{2+{\cfrac {1}{\frac {5}{2}}}}}}$

${\displaystyle =1+{\cfrac {1}{2+{\cfrac {2}{5}}}}=1+{\cfrac {1}{\frac {12}{5}}}}$

${\displaystyle =1+{\cfrac {5}{12}}={\frac {17}{12}}}$

Наконец (шаг 3), делим числитель на знаменатель рациональной дроби, чтобы получить приближённое значение корня:

${\displaystyle 17\div 12=1{,}42}$ округлено до трёх знаков.

Действительное значение корня √2 равно 1,41 с точностью до трёх значащих цифр. Относительная ошибка равна 0,17 %, так что рациональная дробь хороша почти до трёх знаков. Если брать больше знаменателей, получим последовательное улучшение приближения — четыре знаменателя дают дробь ${\displaystyle {\frac {41}{29}}=1{,}4137}$, что даёт почти 4 цифры точности, и т. д.

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

√S	цепная дробь	~десятичное	Подходящие дроби
√2	${\displaystyle [1;{\overline {2}}]}$	1,41421	${\displaystyle {\frac {3}{2}},{\frac {7}{5}},{\frac {17}{12}},{\frac {41}{29}},{\frac {99}{70}}}$
√3	${\displaystyle [1;{\overline {1,2}}]}$	1,73205	${\displaystyle {\frac {2}{1}},{\frac {5}{3}},{\frac {7}{4}},{\frac {19}{11}},{\frac {26}{15}},{\frac {71}{41}},{\frac {97}{56}}}$
√5	${\displaystyle [2;{\overline {4}}]}$	2,23607	${\displaystyle {\frac {9}{4}},{\frac {38}{17}},{\frac {161}{72}}}$
√6	${\displaystyle [2;{\overline {2,4}}]}$	2,44949	${\displaystyle {\frac {5}{2}},{\frac {22}{9}},{\frac {49}{20}},{\frac {218}{89}}}$
√10	${\displaystyle [3;{\overline {6}}]}$	3,16228	${\displaystyle {\frac {19}{6}},{\frac {117}{37}}}$
${\displaystyle {\sqrt {\pi }}}$	${\displaystyle [1;1,3,2,1,1,6...]}$	1,77245	${\displaystyle {\frac {2}{1}},{\frac {7}{4}},{\frac {16}{9}},{\frac {23}{13}},{\frac {39}{22}}}$
${\displaystyle {\sqrt {e}}}$	${\displaystyle [1;1,1,1,5,1,1...]}$	1,64872	${\displaystyle {\frac {2}{1}},{\frac {3}{2}},{\frac {5}{3}},{\frac {28}{17}},{\frac {33}{20}},{\frac {61}{37}}}$
${\displaystyle {\sqrt {\phi }}}$	${\displaystyle [1;3,1,2,11,3,7...]}$	1,27202	${\displaystyle {\frac {4}{3}},{\frac {5}{4}},{\frac {14}{11}}}$

Примечание: Перечислены все подходящие дроби вплоть до знаменателя 99.

В общем виде чем больше знаменатель рациональной дроби, тем лучше аппроксимация. Также можно доказать, что отсечение непрерывной дроби приводит к рациональной дроби, с лучшим приближением к корню любой дроби со знаменателем, меньшим или равным знаменателю этой дроби. Например, никакая дробь со знаменателем, не превосходящем 70, не будет так же хороша, как аппроксимация к √2 числом 99/70.

Метод последовательности Люка

Последовательность Люка первого рода ${\displaystyle U_{n}(P,Q)}$ определяется рекуррентным отношением

${\displaystyle U_{n}(P,Q)={\begin{cases}0&{\text{если }}n=0\\1&{\text{если }}n=1\\P\cdot U_{n-1}(P,Q)-Q\cdot U_{n-2}(P,Q)&{\text{в противном случае}}\end{cases}}}$

и его характеристическим многочленом является

${\displaystyle x^{2}-P\cdot x+Q=0}$

, он имеет дискриминант ${\displaystyle D=P^{2}-4Q}$ и корни

${\displaystyle {\begin{matrix}x_{1}={\dfrac {P+{\sqrt {D}}}{2}},&x_{2}={\dfrac {P-{\sqrt {D}}}{2}}\end{matrix}}}$

Всё это даёт следующее положительное значение

${\displaystyle \lim _{n\to \infty }{\dfrac {U_{n+1}}{U_{n}}}=x_{1}}$

. Так что если мы хотим получить ${\displaystyle {\sqrt {a}}}$, мы можем выбрать ${\displaystyle P=2}$ и ${\displaystyle Q=1-a}$, а затем вычислить ${\displaystyle x_{1}=1+{\sqrt {a}}}$ используя ${\displaystyle U_{n+1}}$ и ${\displaystyle U_{n}}$для больших значений ${\displaystyle n}$. Наиболее эффективный способ вычисления ${\displaystyle U_{n+1}}$ и ${\displaystyle U_{n}}$ —

${\displaystyle {\begin{bmatrix}U_{n}\\U_{n+1}\end{bmatrix}}={\begin{bmatrix}0&1\\-Q&P\end{bmatrix}}\cdot {\begin{bmatrix}U_{n-1}\\U_{n}\end{bmatrix}}={\begin{bmatrix}0&1\\-Q&P\end{bmatrix}}^{n}\cdot {\begin{bmatrix}U_{0}\\U_{1}\end{bmatrix}}}$

Итог:

${\displaystyle {\begin{bmatrix}0&1\\a-1&2\end{bmatrix}}^{n}\cdot {\begin{bmatrix}0\\1\end{bmatrix}}={\begin{bmatrix}U_{n}\\U_{n+1}\end{bmatrix}},}$

а тогда при ${\displaystyle n\to \infty }$:

${\displaystyle {\sqrt {a}}={\frac {U_{n+1}}{U_{n}}}-1}$

Аппроксимации, зависящие от представления в виде числа с плавающей запятой

Число представляется в виде числа с плавающей запятой как ${\displaystyle m\times b^{p}}$. Этот формат записи называется также экспоненциальной записью. Квадратный корень из этого числа равен ${\displaystyle {\sqrt {m}}\times b^{\tfrac {p}{2}}}$ и аналогичные формулы могут быть представлены для кубических корней и логарифмов. Это не упрощает задачу, но если требуется только аппроксимация, то ${\displaystyle b^{\tfrac {p}{2}}}$ является хорошей оценкой порядка мантиссы. Далее, понимаем, что некоторые степени p могут оказаться нечётными, тогда для ${\displaystyle 3141{,}59=3{,}14159{\times }10^{3}}$ вместо работы с дробными степенями основания умножаем на него и вычитаем единицу из степени, делая её чётной. Уточнённое представление превращается в ${\displaystyle 31{,}4159{\times }10^{2}}$, так что квадратный корень будет равен ${\displaystyle {\sqrt {31,4159}}{\times }10^{1}}$.

Если взять лишь целую часть мантиссы, она может принимать значения от 1 до 99 и это можно использовать в качестве индекса в таблице из 99 предварительно вычисленных корней для завершения оценки. Компьютер, использующий шестнадцатеричное основание может потребовать большей таблицы, но при использовании основания 2 таблица будет состоять лишь из трёх величин — возможными битами целой части уточнённого представления мантиссы могут быть 01 (если степень чётная, так что нет никакого сдвига, и заметим, что нормализованное число с плавающей точкой всегда имеет ненулевую старшую цифру), или, если степень была нечётной, 10 или 11, это два первых бита исходной мантиссы. Тогда 6,25 (= 110,01 в двоичном представлении) нормализуется к ${\displaystyle 1{,}1001\times 2^{2}}$ с чётной степенью, так что парой битов мантиссы будет 01, в то время как 0,625 (= 0,101 в двоичном представлении) нормализуется к ${\displaystyle 1{,}01\times 2^{-1}}$ с нечётной степенью, так что требуется преобразование числа к ${\displaystyle 10{,}1\times 2^{-2}}$, а тогда парой бит будет 10. Заметим, что младший бит порядка отражается в старший бит сгруппированной парами мантиссы. Чётная степень имеет нулевой младший бит и уточнённая мантисса будет начинаться с нуля, в то время как нечётная степень имеет 1 в младшем бите и уточнённая мантисса будет начинаться с 1. Таким образом, когда степень делится пополам, это эквивалентно тому, что младший бит порядка сдвигается в первый бит попарно сгруппированной мантиссы.

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

Многие компьютеры следуют стандарту IEEE для чисел с плавающей запятой^[англ.] (или достаточно близкое представление) и очень быстрое приближение для квадратного корня может быть получено в качестве стартового значения метода Ньютона. Техника данного приближения вытекает из факта, что формат плавающего числа (по основанию два) аппроксимирует логарифм по основанию 2. То есть, ${\displaystyle \log _{2}(m\times 2^{p})=p+\log _{2}(m)}$

Так что для 32-битного числа с плавающей запятой в формате IEEE (в котором степень имеет смещение^[англ.] на 127^[16]) вы можете получить приближённый логарифм путём интерпретации числа как 32-битного целого, умножения его на ${\displaystyle 2^{-23}}$ и вычета смещения 127, то есть

${\displaystyle x_{\text{int}}\cdot 2^{-23}-127\approx \log _{2}(x).}$

Например, число 1,0 в шестнадцатеричной системе имеет вид 0x3F800000, что можно представить как ${\displaystyle 1065353216=127\cdot 2^{23}}$, если рассматривать его как целое. Используя вышеприведённую формулу вы получите ${\displaystyle 1065353216\cdot 2^{-23}-127=0}$, как и ожидалось от ${\displaystyle \log _{2}(1{,}0)}$. Аналогичным образом вы получите 0,5 из 1,5 (=0x3FC00000).

Чтобы получить квадратный корень, делим логарифм на 2 и преобразуем результат обратно. Ниже программа демонстрирует идею. Заметим, что младший бит порядка намеренно переводится в мантиссу. Одним из способов обоснования шагов этой программы, в предположении что ${\displaystyle b}$ является смещением степени, а ${\displaystyle n}$ является числом запоминаемых бит в мантиссе, заключается в доказательстве

${\displaystyle (((x_{\text{int}}/2^{n}-b)/2)+b)\cdot 2^{n}=(x_{\text{int}}-2^{n})/2+((b+1)/2)\cdot 2^{n}.}$

/* Предполагаем, что плавающее число имеет формат IEEE 754 */
#include <stdint.h>
float sqrt_approx(float z)
{
	union { float f; uint32_t i; } val = {z};	/* Преобразуем тип не меняя битового представления */
	/*
	 * Для обоснования работы кода докажите, что
	 * ((((val.i / 2^m) - b) / 2) + b) * 2^m = ((val.i - 2^m) / 2) + ((b + 1) / 2) * 2^m)
	 * где
	 * b = смещение степени
	 * m = число бит в мантиссе
	 */
	val.i -= 1 << 23;	/* Вычитаем 2^m. */
	val.i >>= 1;		/* Делим на 2. */
	val.i += 1 << 29;	/* Добавляем ((b + 1) / 2) * 2^m. */

	return val.f;		/* Интерпретируем снова как плавающее */
}

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

	val.i = (1 << 29) + (val.i >> 1) - (1 << 22) + a;

где a — смещение для уменьшения ошибок аппроксимации. Например, с a = 0 результаты точны для чётных степеней двойки 2 (например, 1,0), но для других чисел результат будет несколько великоват (например, 1,5 для 2,0 вместо 1,414... с ошибкой 6 %). При a = −0x4B0D2 максимальная относительная ошибка сокращается до ±3,5 %.

Если приближение нужно использовать как начальное значение для метода Ньютона в уравнении ${\displaystyle (1/x^{2})-S=0}$, то обратная форма, показанная в следующем разделе, предпочтительнее.

Обратное значение квадратного корня

Основная статья: Быстрый обратный квадратный корень

Вариант описанной выше процедуры представлен ниже и он может быть использован для вычисления обратного к квадратному корню, то есть ${\displaystyle x^{-{1 \over 2}}}$. Этот вариант написал Грег Уолш. Приближение сдвигом даёт относительную ошибку менее 4 % и ошибка уменьшается до 0,15 % после одной итерации метода Ньютона^[17]. В компьютерной графике это очень эффективный способ нормализации вектора.

float invSqrt(float x) {
    float xhalf = 0.5f * x;
    union {
        float x;
        int i;
    } u;
    u.x = x;
    u.i = 0x5f375a86 - (u.i >> 1);
    /* Следующая строка может быть повторена произвольное число раз для увеличения точности */
    u.x = u.x * (1.5f - xhalf * u.x * u.x);
    return u.x;
}

Некоторые СБИС реализуют нахождение обратной величины к квадратному корню с помощью полиномиальной оценки с последующей итерацией Голдшмидта^[18].

Корень из отрицательного или комплексного числа

Если ${\displaystyle S<0}$, то его главный корень равен

${\displaystyle {\sqrt {S}}={\sqrt {\vert S\vert }}\,\,i\,.}$

Если ${\displaystyle S=a+bi}$, где a и b вещественные числа и ${\displaystyle b\neq 0}$, то его главный корень равен

${\displaystyle {\sqrt {S}}={\sqrt {\frac {\vert S\vert +a}{2}}}\,+\,\operatorname {sgn}(b){\sqrt {\frac {\vert S\vert -a}{2}}}\,\,i\,.}$

Это можно проверить возведением в квадрат^[19][20]. Здесь

${\displaystyle \vert S\vert ={\sqrt {a^{2}+b^{2}}}}$

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

См. также

• Алгоритм альфа max плюс бета min^[англ.]
• Целочисленный квадратный корень

Примечания

1. Кроме главного корня имеется отрицательный квадратный корень, равный по модулю главному корню, но с противоположным знаком, за исключением случая нуль, когда имеется два одинаковых корня, равных нулю.
2. Множители два и шесть используются ввиду того, что они аппроксимируют среднее геометрическое нижнего и верхнего возможных значений с заданным числом знаков: ${\displaystyle {\sqrt {{\sqrt {1}}\cdot {\sqrt {10}}}}={\sqrt[{4}]{10}}\approx 1{,}78\,}$ и ${\displaystyle {\sqrt {{\sqrt {10}}\cdot {\sqrt {100}}}}={\sqrt[{4}]{1000}}\approx 5{,}62\,}$.
3. Неокруглённая оценка имеет максимальную абсолютную ошибку 2,65 в точке 100 и максимальную относительную ошибку в 26,5% в точках y=1, 10 и 100
4. Если число находится ровно посередине между двумя квадратами, наподобие 30,5, берём большее число, которое в нашем случае 6
5. Это уравнение касательной прямой к y=x² в точке y=1.
6. Fowler, Robson, 1998, с. 376.
7. Heath, 1921, с. 323–324.
8. Bailey, Borwein, 2012, с. 646–657.
9. Bucking down to the Bakhshali manuscript  (неопр.). Simply Curious blog (5 июня 2018). Дата обращения: 21 декабря 2020. Архивировано 26 октября 2020 года.
10. Fast integer square root by Mr. Woo's abacus algorithm (archived)
11. Integer Square Root function  (неопр.). Дата обращения: 30 декабря 2021. Архивировано 30 сентября 2007 года.
12. Wilkes, Wheeler, Gill, 1951.
13. Campbell-Kelly, 2009.
14. Gower, 1958, с. 142–143, 1958.
15. Markstein, Peter (November 2004). Software Division and Square Root Using Goldschmidt's Algorithms (PDF). 6th Conference on Real Numbers and Computers. Dagstuhl, Germany. CiteSeerX 10.1.1.85.9648. Архивировано (PDF) 28 апреля 2022. Дата обращения: 30 декабря 2021.
16. К экспоненте числа добавляется 127, что позволяет интерпретировать экспоненту как число без знака.
17. Fast Inverse Square Root Архивная копия от 6 февраля 2009 на Wayback Machine by Chris Lomont
18. "High-Speed Double-Precision Computation of Reciprocal, Division, Square Root and Inverse Square Root" by José-Alejandro Piñeiro and Javier Díaz Bruguera 2002 (abstract)
19. Abramowitz, Stegun, 1964, с. 17.
20. Cooke, 2008, с. 59.

Литература

• David Fowler, Eleanor Robson. Square Root Approximations in Old Babylonian Mathematics: YBC 7289 in Context // Historia Mathematica. — 1998. — Т. 25, вып. 4. — doi:10.1006/hmat.1998.2209.
• Thomas Little Heath. A History of Greek Mathematics. — Oxford: Clarendon Press, 1921. — Т. 2. — С. 323–324.
• David Bailey, Jonathan Borwein. Ancient Indian Square Roots: An Exercise in Forensic Paleo-Mathematics // American Mathematical Monthly. — 2012. — Т. 119, вып. 8.
• Miltonn Abramowitz, Irene A. Stegun. Section 3.7.26 // Handbook of mathematical functions with formulas, graphs, and mathematical tables. — Courier Dover Publications, 1964. — С. 17. — ISBN 978-0-486-61272-0.
• J. C. Gower. A Note on an Iterative Method for Root Extraction // The Computer Journal. — 1958. — Т. 1 1, вып. 3.
• M. Campbell-Kelly. Origin of Computing // Scientific American. — 2009. — Сентябрь.
• Roger Cooke. Classical algebra: its nature, origins, and uses. — John Wiley and Sons, 2008. — ISBN 978-0-470-25952-8.
• M. V. Wilkes, D. J. Wheeler, S. Gill. The Preparation of Programs for an Electronic Digital Computer. — Addison-Wesley, 1951.

СсылкиWeisstein, Eric W. Square root algorithms&nbsp;(англ.) на сайте Wolfram MathWorld.

• Square roots by subtraction
• Integer Square Root Algorithm by Andrija Radović
• Personal Calculator Algorithms I : Square Roots (William E. Egbert), Hewlett-Packard Journal (may 1977) : page 22
• Калькулятор для обучения квадратному корню