Алгоритм извлечения квадратного корня. Построение, доказательство и реализация.

Пусть $a$ -- натуральное число, являющееся квадратом другого натурального числа $b$:
$$
a,b\in \mathbb{N},\ a=b^2.
$$

Рассмотрим случай, когда $a$ известно, а $b$ -- нет.
Для нахождения $b$ мы применим алгоритм, описанный здесь:
http://spacemath.xyz/algoritm-izvlecheniya-kvadratnogo-kornya/

Посмотрим, как устроен этот алгоритм.

Прежде всего, определим количество цифр в числе $b$. Пусть $a$ содержит $k$ цифр,
и пусть
$$
\begin{equation*}
p =
\begin{cases}
\frac{k}{2}, &\text{если $k$ чётно,}
\\
\frac{k+1}{2}, &\text{если $k$ нечётно.}
\end{cases}
\end{equation*}
$$
Покажем, что $p$ есть количество цифр в $b$.
Действительно, наименьшее число, содержащее $(p+1)$ цифру, есть $10^p$.
Его квадрат,
$$
(10^p)^2 = 10^{2p},
$$
содержит (2p+1) цифру. В зависимости от чётности числа $k$, $a$ содержит либо $2p$
цифр, либо $(2p-1)$ цифру. Это значит, что $a\lt (10^p)^2$.
Следовательно, количество цифр в $b$ меньше или равно $p$.
Наибольшее число, содержащее $(p-1)$ цифру, есть $10^{(p-1)}-1$, и его квадрат содержит лишь $(2p-2)$ цифры, поэтому
$$
(10^{(p-1)}-1)^2 \lt a.
$$
Следовательно, количество цифр в $b$ больше или равно $p$.

$p$ представляет собой количество граней числа $a$.

Теперь, когда нам известно $p$ --количество цифр числа $b$, мы видим, что
$$
10^{(p-1)}\le b \le (10^p -1).
$$

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

Мы последовательно найдём все цифры числа $b$, начиная со старшего разряда $(10^{(p-1)})$.
В качестве первого приближения $b_1$ числа $b$ возьмём $m_{p-1} \cdot 10^{(p-1)}$), где
$$
m_{p-1} = \max_{i\in \{ 1,\ 2,\ \ldots,\ 9 \},\ (i\cdot 10^{(p-1)})^2 \le a}\ i.
$$

$m_{p-1}$ будет точным значением старшего разряда числа $b$. В самом деле,
в старшем разряде $b$ не может находиться какое-нибудь $n$, для которого $((n+1) \cdot 10^{(p-1)})^2 \le a$,
так как
$$
n \cdot 10^{(p-1)} \lt n \cdot 10^{(p-1)} + (10^{(p-1)}-1) \lt (n+1) \cdot 10^{(p-1)},
$$
и квадрат числа $(n \cdot 10^{(p-1)} + (10^{(p-1)}-1))$ был бы строго меньше $a$
(несмотря на то, что мы все оставшиеся разряды заполнили девятками).

Найдём значение для следующего разряда $m_{p-2}$ числа $b$ ($10^{(p-2)}$):
$$
m_{p-2} = \max_{i\in \{ 0,\ 1,\ \ldots,\ 9 \},\ ((m_{p-1}\cdot 10 +i)\cdot 10^{(p-2)})^2 \le a}\ i.
$$
В разряде $10^{(p-2)}$ так же не может быть значение меньше этого максимума.
Если $n_1\lt m_{p-2}$, то
$$
((m_{p-1}\cdot 10 + n_1)\cdot 10^{(p-2)} + 10^{(p-2)}-1)^2 \lt ((m_{p-1}\cdot 10 + m_{p-2})\cdot 10^{(p-2)})^2 \le a.
$$
Получается, что квадрат такого числа будет меньше $a$, даже если мы заполним все разряды младше $10^{(p-2)}$ девятками.

Вообще, если мы имеем приближение $b_1$ числа $b$, такое, что значения первых $j$ разрядов $b_1$ совпадают c $b$, а остальные -- равны нулю, то мы можем найти следующий ($(j+1)$-ый, если считать слева направо, начиная с 1) разряд числа $b$ так:
$$
m_{p-(j+1)} = \max_{i\in \{ 0,\ 1,\ \ldots,\ 9 \}, \ (b_1 + i \cdot 10^{p-(j+1)})^2 \le a}\ i.
$$

Теперь усложним немного этот алгоритм.

Будем использовать те же обозначения, что и раньше.
Пусть
$$
t_1 = m_{p-1},
$$
и для $k\gt 1$
$$
t_k = t_{k-1}\cdot 10 + m_{p-k}.
$$
$b_i$ -- приближения для $b$:
$$
b_i = t_i\cdot 10^{p-i}.
$$
Грубо говоря, первые $i$ цифр (считая слева направо) числа $b_i$ совпадают с первыми $i$ цифрами $b$, остальные цифры -- нули.
Пусть для $i\in \{2,\ \ldots, \ p \}$
$$
s_i = t_{i-1}\cdot 10 + x.
$$
В предыдущем алгоритме, чтобы найти очередной $t_i,\ i\gt 1$, мы искали
$$
m_{p-i} = \max_{x=\{0,\ \ldots,\ 9\}, \ (s_i\cdot 10^{(p-i)})^2 \le a^2} x.
$$
Теперь мы на первом шаге вычислим разность
$$
r_1=a-b_1^2,
$$
на втором шаге вычислим
$$
r_2 = r_1 - (b_2-b_1),
$$
и в общем виде на шаге $i, \ i=2,\ \ldots, p\ $ будет вычислена
разность
$$
r_i = r_{i-1} - (b_{i}-b_{i-1}),
$$
а $m_{p-i}$ выбираться будет так:
$$
m_{p-i} = \max_{x=\{1,\ \ldots,\ 9\}, \ (s_i-b_{i-1})^2 \le r_{i-1}} x.
$$
На каждом шаге $(s_i\cdot 10^{(p-i)}-b_{i-1})^2$ кратно $10^{2(p-i)}$. Поэтому последние $(p-i)$ граней не влияют на вычисление, и мы "сносим" их постепенно.

Реализация

<?php

echo 12323948**2 === 151879694306704;
echo "\n";

$a = '151879694306704';
$digitsNumber = strlen($a);
// Количество граней:
$fringesNumber = intdiv($digitsNumber, 2) + $digitsNumber % 2;

$result = '';

for ($i = 1; $i <= $fringesNumber; $i++) {

    $zeroes = str_repeat('0', $fringesNumber - $i);
    for ($digit = 1; $digit <= 10; $digit++) {

        $approximation = strval($result . $digit) . $zeroes;

        if (bccomp(bcmul($approximation, $approximation), $a) > 0) {
            $resultDigit = $digit - 1;
            break;
        }
    }

    $result .= strval($resultDigit);
}

echo "sqrt($a)=$result\n";

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

<?php

echo 1518796943**2 === 2306744154066145249;
echo "\n";

$a = '2306744154066145249';
// Количество цифр
$digitsNumber = strlen($a);
// Количество граней
$fringesNumber = intdiv($digitsNumber, 2) + $digitsNumber % 2;

// Массив граней
$fringes = [];
for ($i = 1; $i <= $fringesNumber; $i++) {
    $fringes[$fringesNumber - $i + 1] = ($digitsNumber - 2 * $i >= 0) ? substr($a, $digitsNumber - 2 * $i, 2) : substr($a, 0, 1);
}

$result = '';
$reminder = '';

/**
 * Возвращает (20 * $result + $digit) * $digit
 * 
 * @param string $result
 * @param int $digit
 * @return string
 */
function sumCalculation(string $result, int $digit): string
{
    return bcmul(bcmul(empty($result) ? '0' : $result, '2') . strval($digit), strval($digit));
}

for ($i = 1; $i <= $fringesNumber; $i++) {

    $reminder .= $fringes[$i];

    for ($digit = 1; $digit <= 10; $digit++) {

        // Как только мы превысили остаток, выбираем предыдущее число
        if (bccomp(sumCalculation($result, $digit), $reminder) > 0) {
            $reminder = bcsub($reminder, sumCalculation($result, $digit-1));
            $resultDigit = $digit - 1;
            break;
        }

    }

    $result .= strval($resultDigit);
}

echo "sqrt($a)=$result\n";