Можно ограничиться нечетными числами: для четного числа, если количество нулевых младших разрядов нечетно, то корня нет, а если четно, то можно сдвинуть число вправо, посчитать корень от нечетного, и сдвинуть обратно влево на половину от первоначального количества нулевых бит.
Для нечетного N и 2k, k > 3, если N ? 1 mod 8, то есть 4 разных корня по модулю 2k, а иначе корней нет. Нам нужен наименьший из этих четырех корней x. При этом другие три корня это 2k - x, 2k-1 + x и 2k - 2k-1 - x
Хочется что-то подобное вычислению обратного по модулю 2k — удваивая количество верных бит за итерацию.
Пусть у нас уже есть корень x0 из N по модулю 2k: N - x02 = 2ka
И мы хотим найти x1 = x0 + 2k-1y, такое чтобы в N - x12 было больше младших нулевых бит.
N - (x0 + 2k-1y)2 = 2ka - 2kx0 * y - 22k-2y2
Поделим на 2k: a - x0 * y - 2k-2y2
И приравняем к 0 по модулю 2k-2: y = a * x0-1 mod 2k-2
Получилии x1 = x0 + 2k-1a * (x0-1 mod 2k-2)
И окончательно x1 = x0 + (N - x02)/2 * (x0-1 mod 2k-2)
Из k бит на следующей итерации получится 2(k-1) бит. Параллельно считаем на каждой итерации обратное к корню.
Тестовый код:
uint8_t sqr16(uint16_t n) {
if (n % 8 != 1) return 0;
uint16_t sqr = (n + 1) / 2; //4 bit
uint16_t inv = 2 - sqr;
sqr += inv * (n-sqr*sqr)/2; //6 bit
inv *= 2 - sqr * inv;
sqr += inv * (n-sqr*sqr)/2; //10 bit
//inv *= 2 - sqr * inv;
if (sqr & 256)
sqr = 0u - sqr;
sqr = (uint8_t)sqr; // lowest root
if (n == sqr*sqr) return sqr;
return 0;
}
Добавив пару итераций, получим корень из uint_64
Комментарии (13)
daiver19
01.10.2019 01:43Двоичный поиск, конечно, чуть медленее, зато может искать floor/ceil квадратного корня и не зависит от типа данных.
Smog1on1the1water
01.10.2019 13:17Зачем изобретать велосипед, когда все уже придумано до нас?
Генри Уоррен, «Алгоритмические трюки для программистов», на основе алгоритма, описанного в 1945 году!!! Джоном фон Нейманном:
int isqrt(unsigned x) { unsigned m, y, b; m = 0x40000000; y = 0; while (m != 0) { b = y | m; y >>= 1; if (x >= b) { x -= b; y |= m; } m >>= 2; } return y; }
Группу операторов в if можно даже упростить, используя знаковый сдвиг вправо на 31 бит, избежав одной проверки и условного перехода:
int t = (int) (x | ~(x - b)) >> 31; // -1, если x >= b, иначе 0 x -= (b & t); y |= (m & t);
Там же в книге есть и другие варианты поиска целочисленных корней (в том числе с разбором варианта идеи, предложенной топикстартером), в том числе и кубических. И много другого вкусного.
Код выше приведен как есть в книге, его нужно будет причесать в плане разрядности в зависимости от компилятора и операционной системы.
4eyes
01.10.2019 16:35-1Заморачивался оптимизацией получения квадратного корня, когда участвовал в AI Cup. Пробовал фишки «из Quake», приблизительный обратный корень SSE и целочисленные приближения. На современном процессоре (у меня haswell) оно совершенно не стоит того. Деталей не помню, но ускорение дает только приблизительный обратный корень на SSE, и это ускорение несущественно — процентов 10-20 ценой потенциальных багов на потере точности и округлениях.
sasha_alesin
01.10.2019 17:53Понятно, что можно реализовать метод Ньютона в целых числах, но он требует деления на каждом шаге. А нельзя ли по другому?
А то есть операция взятия остатка не использует внутри себя операцию деления? Физически же остаток считается делением в столбик, или нет?
ZiggiPop
Забавно, но именно сегодня задался решением той же проблемы — является ли целое число квадратом (простенькая задачка в codewars). В результате прикинул, что самым дешевым способом может оказаться анализ числа по признакам, приведенным в википедии:
wundte
А можете дать ссылку на ваше решение. Недавно сам делал эту задачу.
paluke Автор
Ну это необходимые, но не достаточные условия
WhiteBlackGoose
Вот именно. Зная их, можно лишь отсеять точно-не-квадраты.
Ketovdk
я конечно понимаю, что мой метод работает медленнее, зато надежный как швейцарские часы и простой:
взять sqrt(x), округлить, возвести в квадрат, проверить, что совпадает с x.
Ну и в статье есть не очевидные утверждения:
Вроде как это даже правда и даже кем-то доказано, но догадаться до этого самостоятельно сложно. И не понятно, что делать с числами, у которых корня нет
WhiteBlackGoose
Взять sqrt очень тяжко, как и делать округление. paluke, сделайте пожалуйста тест производительности, для round(sqrt(x)) и вашей функции. Интересно было б посмотреть
Ketovdk
ну я указал, что это медленнее и зависит от того, насколько критична производительность в этом случае, я скорее к тому, что преждевременная оптимизацией заниматься не стоит. Тем более, что в приведенном способе все-равно придется это сделать, если обнаружится, что корня по модулю нет и в худшем случае производительность от этой оптимизации станет только хуже на несколько проверок. Если нужно прям очень производительно, то нужно просто преподсчитать все квадраты и искать бинарным поиском за O(log(sqrt(n))
paluke Автор
На х86 при включении оптимизации sqrt будет инлайниться компилятором в sqrtsd (sse) или fsqrt (x87), и это на самом деле быстрее.
Этот метод может быть быстрее для BigInteger чисел. Ну или где-нибудь в embedded, где нет сопроцессора (правда непонятно, зачем в embedded может понадобиться корень).
Протестировал. sse чуть быстрее, х87 чуть медленее