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

[B@R5uk]

Andrey A, в связи с вашим предложением применять цепные дроби для нахождения корня меня мучает вопрос: а какова средняя вычислительная сложность нахождения первого нетривиального решения уравнения Пелля для чисел в диапазоне $0...2^N-1$ ? А так же какова сложность проверки того факта, что эти самые нетривиальные решения вообще есть?

[Andrey A]

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

... что эти самые нетривиальные решения вообще есть?

Тривиальное решение одно: $1^2-m\cdot 0^2=1$, а наличие бесконечной последовательности нетривиальных решений для любого $m\neq c^2$ - факт установленный. Доказательство следует из самой классической процедуры разложения квадратного радикала и конечного числа кв. вычетов $z_i=p_i^2-mq_i^2<m/2, где $\frac{p_i}{q_i}$ - $i$-ая подходящая дробь. Цепная дробь периодична, последовательность вычетов - тоже, а единица заложена уже в "нулевом" тривиальном решении. Вычислительная сложность - об этом как раз Вас хотелось поспрашивать. Я, простите, не математик. Знаю, что числа вида $(ab)^2\pm b$ и $(ab)^2\pm 2b$ имеют короткий период в 2 - 4 знака, но для общего случая статистики не встречал. Предполагаю, что процедура нахождения 1-го решения более затратная, чем хотелось бы, иначе использовалась бы мат. пакетами (хотя - кто знает, тема-то не популярная). Но в любом случае это разовая акция, а далее выигрыш может быть ощутимым, поскольку для заданной точности любой другой метод задействует бОльшие знаменатели, будь это степени десятки или что-то еще. Ну и предложенный метод разложения. Новая хау. Если это требует разъяснений - я готов, там такая тройная рекуррентная последовательность образовывается. Сложение и умножение относительно малых целых чисел, больше ничего. Excel сразу отвечает, возможно иллюзия. Получив последовательность $a_1,a_2,...,a_k$, собираем искомую простую дробь: $u_0=a_k;u_1=a_{k-1}a_k+1;...;u_i=a_{k-i}u_{i-1}+u_{i-2};...$ (два последних члена).
Для грубых приближений смысла нет огород городить, а для точных - вроде бы перспективно.

[Andrey A]

Andrey A в сообщении #1020357 писал(а):

Вычисления ведутся до $z_k=1$, и простая дробь восстанавливается от последнего знака.

P.S. Лучше так: Вычисления ведутся до $z_k=\pm 1$, и простая дробь восстанавливается от последнего знака.
В сущности ведь используется тождество $(p^2-mq^2)^2=(p^2+mq^2)^2-m(2pq)^2$. Если $p^2-mq^2=z$ или $mq^2=p^2-z$, делая соответствующие подстановки получаем $(2p^2-z)^2-m(2pq)^2=z^2$. Для $z_k=-1$, а это не редкий случай, получаем первую дробь $\frac{P_1=2p_k^2+1}{Q_1=2p_kq_k}$, и далее по известной схеме. При больших $P_1,Q_1$ возведение в квадрат (назовем это "схема 2") может оказаться "слишком" быстрой процедурой. Тогда можно использовать схему из первого моего поста (схема 1). Перейти от 1-ой ко 2-ой можно всегда, но не обратно.

[whitefox]

slavav в сообщении #1020130 писал(а):

Я сравнил метод Ньютона/Герона и ваш. Ваш лучше - он делает двоичный логарифм итераций, а Ньютон/Герон больше. Не смотря на квадратичную сходимость последнего.

Вы уверены в правильности такой оценки?

Возьмём, например, 1024-битное число. По Вашей оценке, приведённый ТС алгоритм, выполнит порядка 1024 итераций.

В качестве первого приближения в алгоритме Герона возьмём произвольное 512-битное число. Пусть это приближение имеет только один верный бит. Квадратичная скорость сходимости означает, что после каждой итерации число верных битов удваивается, следовательно потребуется порядка 9 итераций чтобы получить 512 верных битов.

[slavav]

whitefox в сообщении #1022716 писал(а):

slavav в сообщении #1020130 писал(а):

Я сравнил метод Ньютона/Герона и ваш. Ваш лучше - он делает двоичный логарифм итераций, а Ньютон/Герон больше. Не смотря на квадратичную сходимость последнего.

Вы уверены в правильности такой оценки?

Я просто проверил. Для чисел в районе $10^9$ алгоритм делает 18-19 итераций. Начальное значение выбиралось как $\frac{a}{2}$.

whitefox в сообщении #1022716 писал(а):

Возьмём, например, 1024-битное число. По Вашей оценке, приведённый алгоритм выполнит порядка 1024 итераций.

В качестве первого приближения в алгоритме Герона возьмём произвольное 512-битное число. Пусть это приближение имеет только один верный бит. Квадратичная скорость сходимости означает, что после каждой итерации число верных битов удваивается, следовательно потребуется порядка 8 итераций чтобы получить 512 верных битов.

Я не умею считать биты. Давайте оценим метод точно. Итерация:
$x_{i+1} = \frac{x_i + \frac{a}{x_i}}{2}$.

Обозначим $C_i = \frac{x_i}{\sqrt{a}}$. Получим:
$C_{i+1} = \frac{C_i + \frac{1}{C_i}}{2}$.

Ещё одно обозначение: относительная погрешность $\varepsilon_i = C_i - 1$. Получим:
$\varepsilon_{i+1} = \frac{{\varepsilon_i}^2}{2(1 + \varepsilon_i)}$.

Обычно, когда говорят "число верных битов удваивается" то имеют ввиду: $\varepsilon_{i+1} \leqslant {\varepsilon_i}^2$. Последнее утверждение гарантирует высокую скорость сходимости, если $\varepsilon_i$ заметно меньше единицы. А если $\varepsilon_i$ заметно больше единицы, то Герон уточняет только один бит. Это я и имел ввиду, говоря, что он медленно стартует. Он медленно стартует с неточного приближения.

[whitefox]

slavav в сообщении #1022761 писал(а):

Я просто проверил. Для чисел в районе $10^9$ алгоритм делает 18-19 итераций. Начальное значение выбиралось как $\frac{a}{2}$.

Ну что ж, давайте проверим.

Возьмём число 220491584. Это число было выбрано случайно, его битовый размер равен 28 бит.

В качестве первого приближения возьмем любое 14-битное число, например 16383 $(11 1111 1111 1111_2),$ то есть потребуется не более $\left(\lceil\log_2 14\rceil+1\right)=5$ итераций.

Первая итерация даст число 14920.
Вторая — 14849.
Третья — 14848.
Четвёртая — 14848, останавливаемся.

Итого потребовалось 4 итерации. Как у Вас получилось 18 - 19?

[slavav]

whitefox в сообщении #1022786 писал(а):

slavav в сообщении #1022761 писал(а):

Начальное значение выбиралось как $\frac{a}{2}$.

[whitefox]

А зачем? Ведь квадратный корень содержит вдвое меньше цифр.
Неэффективное применение алгоритма ещё не говорит о его неэффективности.

[slavav]

Пытаюсь оправдаться:
На языке высокого уровня трудно найти число бит в числе. Этому посвящена глава 5.3 в Алгоритмических трюках младшего Уоррена.
Там же в главе 11.1 описаны алгоритмы поиска корня на целых числах. Половина главы относится к быстрому получению первого приближения (после этого достаточно 5 делений для 32-битовых целых).

[SSerge]

slavav в сообщении #1022813 писал(а):

На языке высокого уровня трудно найти число бит в числе.

Выделить старший ненулевой бит несложно, сложнее быстро посчитать в какой позиции он оказался.
У процессоров с архитектурой ARM, впрочем, есть очень полезная команда CLZ:

Цитата:

Count Leading Zeros returns the number of binary zero bits before the first binary one bit in a value.

[Dmitriy40]

У процессоров x86 начиная с 80386 есть такая же команда BSR. Не самая быстрая, но быстрее вычислений на С. Да и вообще почти на всех процессорных архитектурах есть. Вот в микроконтроллерах редкость.

[Alexu007]

Квадратный корень.

1. Считываем входное число при osn = 100. В дальнейшем считаем что osn = 10.
2. Подбираем цифру, квадрат которой не превосходит старший разряд входного числа. Это есть первая цифра ответа.
3. Найдем разницу между старшим разрядом входного числа и квадратом найденной цифры. Ответ припишем слева ко второму
разряду входного числа. Получим некое число A.
4. Удвоим имеющуюся часть ответа, получим число a. Подберем наибольшую цифру x, такую что произведение чисел ax
(a - десятки, x - единицы) и x не превосходило А. Число x - вторая цифра ответа.
5. Число A - ax * x припишем слева к третьему разряду входного числа. Получим некое число B. Повторяем итерации
4-5 до тех пор пока не закончатся разряды во входном числе.

Более нагляднее рассмотрим пример на основе числа 138384.

1) Храним число как 84-83-13
2) 3*3 <=13. Текущий ответ: 3
3) 13-3*3 = 4. Припишем 4 к 83. Получаем А = 483.
4) a = 3*2 = 6. 6x * x <=483. x = 7. Текущий ответ: 37.
5)483-67*7 = 14. Припишем 14 к 84. Получаем B = 1484.
6). b = 37*2 = 74. 74x*x <=1484. x = 2. Текущий ответ = 372. Он же и окончательный.

http://cppalgo.blogspot.co.uk/2010/09/blog-post_11.html

[B@R5uk]

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

Однако, в процессе реализации столкнулся с проблемой останова. Дело в том, что если пытаться завершить процедуру при совпадении двух последних приближений, то есть шанс, что она никогда не закончится: приближения будут прыгать туда-сюда на единичку. Такое происходит, если пытаться вычислить корень из чисел вида $N^2-1$. Результатом должно быть число N, но когда это число подставляется в рекуррентную формулу, то получается вот что:

$${{x}_{k-1}}=N$$

$${{x}_{k}}=\left\lfloor \frac{1}{2}\left( {{x}_{k-1}}+\left\lfloor \frac{{{N}^{2}}-1}{{{x}_{k-1}}} \right\rfloor \right) \right\rfloor =\left\lfloor \frac{1}{2}\left( N+\left\lfloor \frac{{{N}^{2}}-1}{N} \right\rfloor \right) \right\rfloor =\left\lfloor \frac{1}{2}\left( N+N-1 \right) \right\rfloor =N-1$$

$${{x}_{k+1}}=\left\lfloor \frac{1}{2}\left( {{x}_{k}}+\left\lfloor \frac{{{N}^{2}}-1}{{{x}_{k}}} \right\rfloor \right) \right\rfloor =\left\lfloor \frac{1}{2}\left( N-1+\left\lfloor \frac{{{N}^{2}}-1}{N-1} \right\rfloor \right) \right\rfloor =\left\lfloor \frac{1}{2}\left( N-1+N+1 \right) \right\rfloor =N$$

Вот, думаю теперь, как побороть эту проблему.

[Dmitriy40]

B@R5uk
Сохраняйте значение предыдущей итерации и проверяйте условие $|x_{k}-x_{k-1}| \le 1$. Где-то видел что можно проверять смену знака производной (например что $x_k<x_{k-1}$), но это кажется зависит от выбора начального приближения.
Какое из двух значений принять за корень легко проверить обратным умножением, если получили больше исходного, значит другое.

Кстати не вижу особой проблемы в выборе начального приближения для не слишком больших чисел на ЯВУ, уж выполнить 15-20 сравнений всяко быстрее даже одного деления (на x86/x64).

-- 16.01.2020, 01:38 --

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

[B@R5uk]

На stackoverflow нашёл вариант решения:

Используется синтаксис Java

public static BigInteger sqrt(BigInteger x) {
    BigInteger div = BigInteger.ZERO.setBit(x.bitLength()/2);
    BigInteger div2 = div;
    // Loop until we hit the same value twice in a row, or wind up alternating.
    for(;;) {
        BigInteger y = div.add(x.divide(div)).shiftRight(1);
        if (y.equals(div) || y.equals(div2))
            return y;
        div2 = div;
        div = y;
    }
}
 

Я, правда, для обычных чисел хотел этот алгоритм приспособить. И я решил проверять условие $0\le A-B^2\le 2B$, где A — число, из которого вычисляется корень, а B — очередное приближение. Таким образом не надо делать две лишних итерации, но требуется дополнительное умножение на каждом шаге. В конце концов (long)Math.sqrt(argValue) оказалось быстрее.

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

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

Можно составить таблицу размером где-нибудь 64-256 чисел (хотя первая четверть не понадобится особо) и из неё выбирать первые несколько бит начального приближения. Делается это с помощью сдвигов как-то так:

Используется синтаксис Java

argBits = 64 - Long.numberOfLeadingZeros(argValue);
shiftSize = (argBits - refBits + 1) / 2;
tmpValue = refRoots[(int) (argValue >> (2 * shiftSize))] << shiftSize;
 

По идее начальный запас точности позволяет сократить количество итераций. Хотя установить в единицу первый бит результата — это уже хороший запас точности.