«The decimal number 585 is 1001001001 in binary. It is palindromic in both bases. Find n-th palindromic number». Или, по-русски: «Десятичное число 585 в двоичной системе счисления выглядит как 1001001001. Оно является палиндромом в обеих системах счисления. Найдите n-й подобный палиндром».
Но само существование значительно более быстрого, с принципиально другой вычислительной сложностью, алгоритма не давало мне покоя, и в конце концов я вернулся к его разбору.
В конце концов, алгоритм оказался не таким уж и сложным, зато, на мой взгляд, очень красивым.
Оригинальное объяснение и реализация используют префиксное дерево, но на мой взгляд, чуть более глубокое понимание математики алгоритма позволяет обойтись более простыми и быстрыми структурами. Думаю, что лучше всего разобрать алгоритм на примере.
Мы будем искать двоичные палиндромы среди десятичных палиндромов шириной 8. Исходных палиндромов ровно 9000: от 10000001 до 99999999.
Теперь возьмём 2 множества чисел:
- А[90] = { 10000001, 11000011, 12000021, ..., 98000089, 99000099 }
- B[100] = { 00000000, 00011000, 00022000, ..., 00988900, 00999900 }
Если описывать эти множества формально, то множество А является подмножеством десятичных палиндромов шириной 8, у которых средние 4 разряда — нули, а множество B состоит из 0, множества десятичных палиндромов шириной 2, умноженных на 103 и множества десятичных палиндромов шириной 4, умноженных на 102.
А если неформально, то множество А — это все возможные «края», а множество B — все возможные «серединки» десятичных палиндромов шириной 8. Очевидно, что множество всех десятичных палиндромов шириной 8 — это множество всех попарных сумм элементов множеств A и B.
for each (a in A) {
for each (b in B) {
palindrome = a + b;
}
}Далее, для краткости я буду использовать «a» вместо «элемент множества A», и «b» вместо «элемент множества B».
Теперь переведём элементы обоих множеств в двоичную систему счисления:
A
10000001 100110001001011010000001
11000011 101001111101100011001011
12000021 101101110001101100010101
…
98000089 101110101110101110011011001
99000099 101111001101001111100100011
B
00000000 000000
00011000 10101011111000
00022000 101010111110000
…
00988900 11110001011011100100
00999900 11110100000111011100
Я выделил по 6 верхних и нижних разрядов у всех a, и 6 нижних разрядов у всех b. Ширина 6 выбрана не случайно — это максимальная степень 2, не превышающая ширину «краёв» A = 102.
Теперь возьмём каждый a, и посмотрим, что общего у всех палиндромов, образованных от него прибавлением b. А общее у них, то, что все они находятся в интервале между самим a и следующим за ним элементом A.
Например, для a = 10000001, все образованные от него десятичные палиндромы { 10000001, 10011001, 10022001, ..., 10988901, 10999901 } находятся в интервале [10000001, 11000011).
Это значит, что все десятичные палиндромы, образованные от a = 10000001 могут начинаться только со следующих 6 двоичных цифр:
100110 — первые двоичные цифры a = 10000001
100111
101000
101001 — первые двоичные цифры a = 11000011
Следовательно, чтобы быть двоичными палиндромами, все эти десятичные палиндромы могут заканчиваться только на обратную перестановку начальных двоичных цифр:
100110 -> 011001
100111 -> 111001
101000 -> 000101
101001 -> 100101
А учитывая, что сам а = 10000001 заканчивается на двоичные цифры 000001, то из всех возможных b нас интересуют только те, которые заканчиваются на двоичные цифры:
011001 — 000001 = 011000
111001 — 000001 = 111000
000101 — 000001 = 000100
100101 — 000001 = 100100
Только для таких b нужно проверять, является ли 10000001 + b двоичным палиндромом.
Используя данный подход, алгоритм для нахождения палиндромов в основаниях BASE0, BASE1 в интервале [1, N] может быть описан следующим образом:
Для каждой ширины W из [1, количество разрядов N в BASE0]
Сгенерировать множества A и B, используя BASE0. Ширина «краёв» A WA = O(W/4), WA ? WB
Перевести элементы A и B в BASE1
Рассортировать элементы B по корзинам по последним цифрам в BASE1
Для каждого a из множества A
Для каждого X из множества возможных начальных цифр a + b
Для каждого b, заканчивающегося на (X — конечные цифры a)
Проверить, является ли a + b палиндромом в BASE1
К сожалению, я уже разучился строго доказывать вычислительную сложность алгоритмов, приведу лишь несколько соображений:
- Множества A и B содержат O(N1/4) элементов.
- Для каждого a в среднем существует не более BASE0 вариантов X.
(Мы изначально выбираем ширину интересующих нас начальных цифр в BASE1 так, чтобы получившееся из них число было не больше BASE0WA, а max(A) больше min(A) в BASE0 раз.) - Для каждого X в среднем проверяется не более, чем BASE1 различных b.
(X распределяется равномерно и почти не коррелирует с конечными цифрами a в BASE1. Значит каждая корзина с b выбирается равномерно, а всего таких корзин не более чем в BASE1 раз меньше, чем самих b.) - Проверка на палиндром все равно занимает O(log(N)).
В общем, я предполагаю, что вычислительная сложность алгоритма O(N1/4 * log(N) * BASE0 * BASE1).
Первая же реализация этого алгоритма вполне предсказуемо была на пару порядков быстрее, а потратив еще немного времени на оптимизацию я довел время вычисления до чуть более 0.01 секунды, более чем в 1000 раз быстрее предыдыщего алгоритма. Этот результат наконец-то меня удовлетворил, но вполне естественно возникло желание проверить его на числах, бо?льших, чем 260. К этому времени я уже обнаружил интересующую меня последовательность в «Энциклопедии целочисленных последовательностей». В списке двойных палиндромов было 122 члена, самый большой из которых состоял из 131 бита. Это было впечатляюще, но и косвенно указывало на то, что никто пока не придумал алгоритм логарифмической сложности.
Вызов был серьёзный — можно было не сомневаться, что люди, получившие такой результат, вложили немало сил в разработку алгоритма и к тому же наверняка были готовы потратить дни и недели машинного времени на последующий поиск. Поэтому, я решил хотя бы приблизительно оценить, время, которое необходимо моему алгоритму на повторение:
2131 / 260 = 271
271 * 1/4 < 218 = 262144
Учитывая 0.01 секунды, необходимых для предела в 260, выходило, что мой алгоритм должен был справиться с этой задачей примерно за 1 час! Я чувствовал какой-то подвох, но вызов был уже принят.
Продолжение следует.
Комментарии (31)
Mrrl
17.12.2015 15:58Я рассмотрел слегка другой алгоритм — рекурсивный.
Когда мы выбрали K старших десятичных цифр палиндрома, мы почти всегда знаем K+1 его старших бит (обычно, почти 3*K, но для наших целей это неважно). Следовательно, знаем K+1 младший бит, а поскольку K младших цифр нам известно, то K+1-я цифра может принимать 5 разных значений — быть либо чётной, либо нечётной (в зависимости от ситуации). Так что число кандидатов увеличилось в 5 раз.
В редких случаях, когда старшие K+1 бит неизвестны, нам придётся поставлять все 10 кандидатов на K+1-ю цифру, но половина из них сразу отсеется.
Какой-то из пунктов 2 или 3 в Ваших рассуждениях вызывает сомнения — в нём (не знаю пока, в каком) должно быть (BASE0/BASE1)^(W/4) вариантов (или (BASE1/BASE0)^(W/4) — из текста непонятно, кто из них 2, а кто 10).Mrrl
17.12.2015 16:04За 20 минут алгоритм нашёл числа до 10^25 (т.е. 80 первых членов последовательности). На поиск чисел до 10^40 ушло бы около 5 лет.
ilyanik
17.12.2015 16:16Я брал BASE0 > BASE1, но не думаю, что это принципиально.
80 первых членов я нахожу за 3.5 секунды.
ilyanik
17.12.2015 16:49Это не строгое доказательство, а попытка объяснить реальные результаты :)
Например, для ширины в 24 десятичных знака, и следовательно 900000 элементов А, проверено 3958921 палиндромов, в среднем по 4.4 различных b на каждый a.Mrrl
17.12.2015 17:35проверено 3958921 палиндромов
Действительно палиндромов? Или значений X? Палиндромов должно быть в 64 раза больше — ведь все числа из множества B в двоичной записи заканчиваются на 6 нулей. Или я чего-то совсем не понял.
Какая там получается ширина края в двоичной записи — 19 или 20?ilyanik
17.12.2015 17:59Ширина края — 19. Т.е. 1000000 палиндромов-серединок распределяются по 2^19 корзин (в среднем по 1.9 в каждой корзине).
Палиндромы-края изменяются от 100000хххххххххххх000001 до 999999хххххххххххх9999999, и для каждого палиндрома-края проверяется всего несколько корзин.Mrrl
17.12.2015 18:59К сожалению, заполнено только 2^13 корзин, по 122 в каждой…
ilyanik
17.12.2015 19:04Эээ..., а почему 2^13 корзин? :)
Это же мой алгоритм, у меня 2^19? Откуда вы вообще взяли число 13?Mrrl
17.12.2015 19:32Но вы же кладёте в них числа из множества B, которые в десятичной системе оканчиваются на 6 нулей, не так ли?
ilyanik
17.12.2015 20:16А, из-за того, что 10 делится на 2?
Да, распределение по корзинам получается очень неравномерное. Но из-за равномерного выбора каждой из 2^19 корзин, в среднем выходит 1.9: на каждую корзину с ~122 элементами 63 пустых.
ilyanik
17.12.2015 18:05Например, первого палиндрома края я проверяю 5 корзин:
100000_000000000000_000001 = 1010100101101000000_1011000111111000010100101011110110100000000000000000000001
100001_000000000000_100001 = 1010100101101000100_0010101000100101111111111010011101111001011000011010100001
ilyanik
17.12.2015 17:58Ширина края — 19. Т.е. 1000000 палиндромов-серединок распределяются по 2^19 корзин (в среднем по 1.9 в каждой корзине).
Палиндромы-края изменяются от 100000хххххххххххх000001 до 999999хххххххххххх9999999, и для каждого палиндрома-края проверяется всего несколько корзин.
Mrrl
17.12.2015 23:45Ещё они пропустили, по меньшей мере, два 39-значных числа (оба начинаются на 12...)…
Mrrl
17.12.2015 23:58Пропущенное 40-значное число: 1017421766189445102992015449816671247101
То есть, их последнее число — не 122-е, а 126-е.
Mrrl
18.12.2015 00:08Ну, и ещё одно 40-значное число: 7114907950920173924554293710290597094117. На всё ушло 49 минут на одном ядре.
ilyanik
18.12.2015 18:37Круто! А как же 5 лет?
Mrrl
18.12.2015 20:45+1Чуток соптимизировал.
Вот ещё:
128 14327425216354951264146215945361252472341 129 31636759764024794204540249742046795763613 130 38090421176450233778487733205467112409083 131 56545858306667087923432978076660385854565 132 98801466348600079992129997000684366410889 133 128795669673344381770077183443376966597821 134 139035351443367699760067996763344153530931 135 170341815153453197154451791354351518143071 136 309612431907274418544445814472709134216903 137 595943598626320807905509708023626895349595 138 905357630732463833436634338364237036753509 139 1816060344791285708869688075821974430606181 140 3381578420704603001613161003064070248751833 141 7342779513978827245484845427288793159772437 142 9101547767757547725021205277457577677451019
44-значных программа не нашла.
Думаю, на этом остановлюсь (до очередного прогрессивного способа) — на следующий разряд ушли бы сутки счёта.ilyanik
18.12.2015 21:57Мне кажется, вам надо принять эстафету «Вперед, на поиски палиндромов» :)
Потому что ваш алгоритм уже быстрее O(N1/4)Mrrl
18.12.2015 22:09До 40-значных включительно — O(N7/34), потом — O(N0.35).
Но алгоритм отличается от вашего только шириной краёв в десятичной записи: я беру не 1/4, а 5/17 от длины числа (пока хватает памяти на серединки). Тогда в непустые корзины действительно попадает по одному числу. Ну, и края фильтруются по младшим битам.ilyanik
18.12.2015 23:09Ну, с шириной краев и середины я тоже «играл», тем более, что потом всё равно кончается память. Но дальше вы по краям тоже ведь не просто проходите…
Mrrl
18.12.2015 23:45Память тратится только на середину. 10^8 чисел удержать можно (если хранить их как 128-битные). Края строятся динамически, их запоминать не нужно.
ilyanik
19.12.2015 11:31Но краев-то 5*10^13 (для ширины 44). Только на их генерацию должно было уйти много часов, даже при скорости генерации 10^9 в секунду.
Mrrl
19.12.2015 14:05С учётом предварительной фильтрации — всего лишь 5^14. Ушло 3-4 часа.
Если бы я не поленился и написал 192-битные числа, получилось бы ещё быстрее.ilyanik
19.12.2015 14:20Так вы профильтровали 5*10^13 краев, или сразу сгенерировали 5^14?
Mrrl
19.12.2015 15:55Можно сказать, что сразу — фильтрация шла для каждой очередной цифры.
ilyanik
20.12.2015 11:32Здорово!
Но, если я не ошибаюсь, ваше ускорение не применимо в случае взаимно простых оснований.
ilyanik
21.12.2015 17:40А вы как-то оптимизировали/сортировали структуру, содержащую серединки?
Я только что посмотрел — на «случайный» доступ к этой огромной структуре тратится 70-80% всего времени, в разы больше чем на всю математику.Mrrl
21.12.2015 18:10У меня основное время уходит на разворот длинных чисел в двоичной записи.
Структура, содержащая серединки, устроена просто:
struct{
UInt128 val;
int next;
}
где next — индекс в массиве следующего элемента в списке. next=0 — ячейка пуста, next=-1 — этот элемент последний. Первые 2^X элементов — корзины (начала списков), дальше — свободные места для продолжения списков. Когда массив оказался длиннее 2ГБ, пришлось переделать его на массив массивов.
Если захочу считать дальше, переделаю его на файл и буду хранить серединки, упорядоченные начиная с младшего бита. Понятие корзины при этом исчезнет, всё будет делаться сопоставлениями двух отсортированных списков. Но это потом.
Mrrl
Предсказываю, что для 131 бита потребуется перебор 5^20 чисел (примерно 10^14).
ilyanik
Почему 520? Только первая цифра десятичного палиндрома нечетная, на остальные ограничений нет.