Квадратный корень на целочисленной арифметике

[Andrey A]

B@R5uk а что, схема $\dfrac{P_{n+1}=2P_n^2-1}{Q_{n+1}=2P_nQ_n}$ себя не оправдала? Смотрю свой пост пятилетней давности, там много лишнего. Получить первое решение Пелля – процедура более затратная чем хотелось бы, это да, но в любом случае это разовая процедура. А дальше кроме умножения никаких чудес. Само название темы предполагает хорошую точность, иначе не стоило бы разговора. Но тогда разовые затраты не существенны. Сравните (только с учетом точности):

$\sqrt{3} \approx \frac{2}{1},\frac{7}{4},\frac{97}{56},\frac{18817}{10864},$ $\frac{708158977}{408855776},\frac{1002978273411373057}{579069776145402304},\frac{2011930833870518011412817828051050497}{1161588808526051807570761628582646656},$ $\frac{8095731360557835890888779535060256832479295062749579257164654370487894017}{4674072680304961790168962360144614650442718636276775741658113370728376064},...$

Всего $7$ итераций. Подходящие дроби – это дроби с оптимальными знаменателями. Иными словами, десятичная дробь той же точности будет значительно длиннее, делайте выводы.

[B@R5uk]

Andrey A, вы совершенно не поняли постановку задачи в этой теме. Вы хотите найти рациональную аппроксимацию для корня как можно точнее. Я же хотел найти целую часть корня — величину $\lfloor\sqrt{A}\rfloor$ — как можно быстрее. Ваша задача тоже интересна и тоже имеет право на жизнь, но я пока не вижу прямой связи между вашей и моей задачами.

[Dmitriy40]

B@R5uk в сообщении #1435473 писал(а):

Я же хотел найти целую часть корня — величину $\lfloor\sqrt{A}\rfloor$ — как можно быстрее.

Пока величина $A<2^{64}$ (т.е. влезает в uint64) на архитектуре x86/x64 практичнее всего пользоваться FPU командой FSQRT, она может брать на вход 64-х битное число без ошибок округления и выдавать корень из него за 10-20 тактов, останется лишь округлить вниз и записать в uint32. Думаю встроенные функции в любом ЯВУ так и делают. Но даже если они и SSE пользуются, с 53-х битовой точностью аргумента, всё равно погрешность результата кажется не превышает 1 и его можно быстро проверить обратным умножением.
Быстрее этого получить корень можно лишь в каких-то исключительных случаях, ИМХО.

[Andrey A]

Ах целую? Ну да, это я из пушки по воробьям. То есть операция нахождения среднего геометрического и, тем более, логарифмы не определены. Но деление с остатком определено? Делим $A$ на некое $m$, получаем неполное частное, берем среднее арифметическое и далее по кругу. Разве не так? Впрочем, я тут не копенгаген.

[Dmitriy40]

Andrey A
Дык Вы алгоритм Герона и описали словами. С ним проблема что при операциях в целых он для некоторых аргументов не останавливается, чередуются результаты $N$ и $N+1$. В этом и был последний вопрос.

(Оффтоп)

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

B@R5uk
Вот тут в примере 4 приведён алгоритм без дефекта с чередованием корней. Но сильно подозреваю что это из-за выбора начального приближения, которое ну совсем не оптимально. А в примере 7 там же показан и несложный выбор более оптимального начального приближения для уменьшения количества итераций.
Просто информация Вам на всякий случай, быстрее встроенной функции это быть не должно.

[Andrey A]

Dmitriy40 в сообщении #1435485 писал(а):

... в этом и был последний вопрос.

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

p.s. хотя нет, для чисел близких к квадрату оно всё равно случается. Но если округлять строго вниз, то в случае кольца $n \to n+1$ целая часть $A$ есть $n$. Такое бывает когда $A$ – квадрат без единицы, и в цикл попадает деление нацело.

[Dmitriy40]

Andrey A
Простите, но именно об этом и был первый пост, целиком. Продолжение видимо бессмысленно:

B@R5uk в сообщении #1435418 писал(а):

В конце концов (long)Math.sqrt(argValue) оказалось быстрее.

[Andrey A]

Dmitriy40
Всё нормально. Просто сунулся не в своё дело, бывает такое.

[Andrey A]

Тут надо бы рассматривать не процедуру нахождения целой части, а процедуру нахождения величины $\sqrt{A}$ с точностью до $k$ десятичных знаков. Тогда искомое находится естественным путем и, думается, тем быстрее чем больше $k$, хотя разница не существенна. Процедура такая:
1) Задаем некоторое $b_1>0$ и вычисляем $(b_1+A/b_1)/2$, округляя до $k$ знаков.
2) Назначаем предыдущий результат $=b_2$ и повторяем процедуру по кругу пока совпадут $k$ десятичных знаков. Для целой части достаточно совпадения первого знака.

Если было уже, не страшно. Повторение – мать учения.

[B@R5uk]

В своих бесконечный поисках откопал на Википедии совсем не то, что искал, а именно — красивый алгоритм извлечения квадратного корня с остатком из длинных чисел. В вики как всегда алгоритм перепечатали из источника с ошибкой, и без всяких дополнительных разъяснений. Хорошо хоть ссылку на оригинальную статью привели. Автор алгоритма, судя по статье, некий Поль Зиммерманн, и почему-то он назвал свой алгоритм "Квадратный корень Карацубы". Но да ладно.

Процедура алгоритма представляет собой следующее.
Вход: одно число n.
Выход: два числа s и r, такие что $n=s^2+r,\,\,\,0\le r<s$.
Шаги алгоритма:
1) если число n достаточно мало (влезает в 32 бита, например), то вычисляем величины $s=\left\lfloor\sqrt{n}\right\rfloor,\,\,\,r=n-s^2$ стандартными средствами и возвращаем результат.
2) представляем исходное число в виде $n=a_3b^3+a_2b^2+a_1b+a_0,\,\,\,0\le a_k<b,\,\,\,b/4\le a_3$
3) рекурсивно вызывая процедуру находим значения s' и r' для числа $n'=a_3b+a_2$ (первая половина исходного числа)
4) находим частное q и остаток u от деления числа $r'b+a_1$ на число $2s'$
5) находим целую часть корня и остаток: $s=s'b+q,\,\,\,r=ub+a_0-q^2$
6) если остаток получился отрицательным ($r<0$), то делаем коррекцию: $r\leftarrow r+2s-1,\,\,\,s\leftarrow s-1$
7) возвращаем результат.

Замечания. Очевидно, что в качестве b при вычислении на компьютере удобно брать степени двойки (или десятки, если работа идёт с числами в десятичном представлении). Не смотря на то, что в оригинале процедура расчёта рекурсивная, её можно свести к обычному циклу. Этому способствует то, что в теле процедуры есть только один самовызов. Так что остаётся только правильно подобрать последовательность чисел b, чтобы идти от малых корней к большим, а не на оборот, как в рекурсии. Автор, видимо, вдохновлялся алгоритмом Карацубы с его "разделяй и властвуй". Но конкретно в случае вычисления корня "властвовать" нужно только над одной половинкой, так что нет необходимости городить рекурсию.

Алгоритм очень сильно перекликается с методом Ньютона расчёта квадратного корня. Например, он так же удваивает число правильных разрядов с каждой "итерацией" (если идти из "дна" рекурсии "наружу"). А так же они оба используют деление больших чисел, что в обоих случаях является самой трудоёмкой операцией и определяет сложность алгоритма. Однако, алгоритм Зиммерманна на каждой итерации выполняет деление не полноразрядных чисел, как метод Ньютона, а чисел экспоненциально убывающей длины (если идти снаружи в глубь рекурсии). Это даёт значительно более высокую производительность. Кроме того, за счёт правильного разбиения и за счёт подсчёта на каждом шаге ровно стольки бит, сколько требуется, алгоритм не имеет никаких проблем с зацикливанием или "переработкой", как метод Ньютона (когда, например, для окончания расчёта не хватило 1-2 верных бит, а следующая итерация "удваивает точность", давая очередные 200 или 200'000 бит).

В качестве иллюстрации можно посчитать целочисленный корень из числа $98765432101112131415161718192021$.
Сначала разбиение:
$98765432\,\,10111213\,\,14151617\,\,18192021$ ($b=10^8$)

$9876\,\,5432\,\,1011\,\,1213$ ($b=10^4$)

$98\,\,76\,\,54\,\,32$ ($b=100$)

$9\,\,8\,\,7\,\,6$ ($b=10$)

Корень из $98$ равен $9$, остаток $17$

Делим $17'7$ на $18=2\cdot 9$, получаем $9$, остаток $15$
Корень из $98'7'6$ равен $9'9$, остаток $15'6-9^2=75\ge 0$

Делим $75'54$ на $198=2\cdot 99$, получаем $38$, остаток $30$
Корень из $9876'54'32$ равен $99'38$, остаток $30'32-38^2=1588\ge 0$

Делим $1588'1011$ на $19876=2\cdot 9938$, получаем $799$, остаток $87$
Корень из $98765432'1011'1213$ равен $9938'0799$, остаток $87'1213-799^2=232812\ge 0$

Делим $232812'14151617$ на $198761598=2\cdot 99380799$, получаем $117131$, остаток $69416279$
Корень из $9876543210111213'14151617'18192021$ равен $99380799'00117131$, остаток $69416279'18192021-117131^2=6941614198520860\ge 0$

Окончательный результат: $98765432101112131415161718192021=9938079900117131^2+6941614198520860$. Чудесным образом нигде не потребовалась коррекция.