Ага, правда, иногда случайности не случайны. Известна история, когда два студента Принстона обнаружили, что комбинация их фамилий звучит по-английски очень неприлично, и решили, что должны обязательно написать совместную статью. Теперь в математике есть машина Кокса-Цукера (Cox-Zucker Machine)
Это, конечно, баловство, но я уместил решение в восемь строчек на питоне! Тут смысла не очень много, это просто формальное преобразование ранее полученных формул с целью уменьшения количества кода.
def cnt(N, b):
f1, f2 = [1]*N, [1]*N
for k in range(1, N//(b-1) + 2):
for n in reversed(range(N)):
f1[n] += sum(f1[n - min(b-1, n, n+1-(k-1)*(b-1)) : n])
for t in range(1, min(2*b-1, n+1, n+2-(k-1)*(b-1))):
f2[n] += (b - abs(b-1-t))*f2[n-t]
return [f2[i] * 4 - f1[i] * 4 + 1 for i in range(N)]
Спасибо, интересная задача! Я немного подумал над ней, но в комментариях хабра мало места, чтобы подробно написать. Примерно так получилось (рассуждения, в принципе, повторяют ваши, но только в более привычных обозначениях): Во-первых, возьмем основание системы счисления произвольное b вместо 10, мне так формулы легче воспринимаются. Для положительных x обозначим сумму цифр Теперь сначала рассмотрим одномерный случай: посчитаем количество таких x, что для всех , будем называть такие числа достижимыми (из нуля). Для этого введем две функции: - количество достижимых чисел в отрезке и , равная единице, если сумма цифр последнего числа в отрезке не превосходит n, иначе нулю. То есть, она говорит, можем ли мы продолжать смотреть следующие отрезки, или надо остановиться. Тогда f и q можно выразить такими формулами:
Теперь точно так же рассмотрим двумерный случай. Удивительно, что условие на то, надо ли нам рассматривать следующий квадрат или нет, сводится к той же самой функции q (не зависящей от размерности!). Я строго не смог доказать это, но там должны быть примерно ваши рассуждения. Для двумерной f формула почти такая же:
Или можно ее сократить до одной суммы, если взять diag(t) равной количеству пар x, y в квадрате bxb, в сумме не превосходящих t:
Теперь, имея одномерную и двумерную функции, можно получить ответ по принципу inclusion-exclusion. Что касается кода, то проверить правильность результата можно обычным поиском (в глубину или ширину, неважно): код на питоне
Вычислить сразу все решения до заданного n можно по этой рекуррентной формуле: код на питоне, вроде бы результаты совпадают с вашими c b=10, N=900.
Сложность алгоритма O(n^2 * b), так как надо итерировать f до значений k=O(n).
В математике вероятность определяется как нормированная мера на сигма-алгебре множеств, в этом контексте событие определяется как измеримое множество. Поэтому явное перечисление множеств (событий) и задание меры на них является фундемантальным и математически корректным способом ввести вероятность.
Спасибо, интересно. Я всегда считал, что естественное обобщение - произведение n-1 вектора в n-менром пространстве как вектор длиной площадь n-1-мерного параллелограмма и перпендикулярный им, и его можно посчитать как определитель матрицы, в которой строки - координаты векторов и базис. Но, согласно википедии, существует много вариантов обобщений, в зависимости от того, какие свойства векторного произведения мы хотим сохранить.
А точно есть такие правила? В Википедии (особенно в Computer Science) много ссылок на форумы и даже личные страницы, когда, например, нет или мало научных статей на эту тему - часто встречается, когда люди делают исследования и не считают нужным их публиковать в авторитетных источниках. Например, в этой статье первые четыре ссылки ведут на форумы и личные страницы (хотя, насколько я знаю, основные результаты по этой теме опубликованы)
Это стандартная задача на поиск в глубину, таких много на литкоде. Самый простой (и медленный) вариант на питоне:
from sys import setrecursionlimit
from textwrap import wrap
setrecursionlimit(100000)
words = open('words.txt').read().split()
def no_repeats(word): return len(set(word)) == len(word)
words = [word for word in words if len(word) == 5 and no_repeats(word)]
def dfs(result, index):
if len(result) >= 25:
print(wrap(result, 5))
return
if index >= len(words):
return
if no_repeats(result+words[index]):
dfs(result+words[index], index+1)
dfs(result, index+1)
dfs("", 0)
Спасибо, интересная статья! Насчет сортировки циклических сдвигов - а вы не пробовали стандартный prefix doubling algorithm - он, насколько я помню, довольно просто пишется, и работает за O(N*log N), сильно быстрее наивного, хотя и медленнее линейных. Но и классический алгоритм Укконена по построению суффиксного дерева на питоне занимает всего 50 строк, хотя самому его написать (и, главное, отладить), наверное, сложно.
не уверен насчет оптимизатора js, но после оптимизатора с++ не будет ни того, ни другого: деление превратится в умножение на reciprocal (throughput = 1 инструкция/такт), а if - в conditional move (2 инструкции/такт). Мне кажется, ни то, ни то не будет узким местом, скорее уж обращение к памяти
я пытался строить все варианты с половиной кубиков (18 для n=8, 22 или 23 для n=9). В вашем случае границы, как на рис.4, сложность в том, что на любой позиции над средней линией можно поставить квадратик размера 1 (допустим, он бы еще не был использован), и тоже получится правильная граница. То есть, "жадное" заполнение сверху вниз и слева направо не работает, надо перебирать больше вариантов, по крайней мере, близко к границе. Или я чего-то не учитываю?
Я немного поэкспериментировал с задачей, но комментарии хабра слишком коротки, чтобы подробно описать результат, поэтому вот очень примерное описание решения задачи для n=8 (n=9 я пока даже не трогал):
Неоптимимизированное решение: стандартный DFS (depth first search with backtracking), состояние представленно как три 16-байтных массива (текущие высоты, длины, и количество оставшихся квадратов). После заполнения текущего нулевого уровня, все заполненные строки схловываются, как в тетрисе, так что заполняется всегда нулевой уровень, что упрощает код.
Время работы: 24 мин, количество шагов DFS: 64G, код
Симметрия: рассматриваются только такие состояния, в которых самый первый квадрат не меньше второго и не меньше последнего в первом ряду, остальные состояния достраиваются с помощью симметрии.
Эвристика: после добавления квадрата, рассматриваются два предыдущих: если средний квадрат в этой тройке ниже, чем крайние, то проверяется, что у нас еще остались достаточно маленьуие квадраты, чтобы в будущем заполнить этот промежуток.
Низкоуровневые оптимизации перехода состояний иcпользуют 128-битные __uint128_t и регистры SSE, а также ускоренное нахождение размеров квадратов, которые можно использовать на текущем шаге, через битовые маски
Еще пробовал, но не получилось: 5. Еще более низкоуровневые оптимизации: все операции в 128-бит или SSE, уменьшение размера состояния до 32 байт, уменьшение дорогостоящих операций обновления состояния: ускорение получается меньше 10%, а код очень сильно усложняется. 6. Я пробовал делать DFS в половину глубины, а потом сшивать попарно. К сожалению, идея не работает: обратный DFS не может сгененировать некоторые варианты (например, когда в первой половине есть "яма" длиной 2 и высотой 4)
Спасибо, интересная задача, крутое исследование! Вы, наверно, можете в oeis.org завести соответствующую последовательность (будет что-то типа 1,0,0,0,0,0,0,2332,216285) Я тоже хочу попробовать ускорить ее на CPU, как будет больше времени)) Тут возможны низкоуровневые оптимизации для представления текущего заполнения, интересно, насколько они будут полезны. Еще можно попробовать такую оптимизацию: при добавлении нового квадрата, если он получается выше предыдущего, проверять, что существуют свободные квадраты, которые можно поставить на предыдущий квадрат (если пред-предыдущий квадрат выше) - но не факт, что это сильно сократит количество вариантов для перебора.
Я написал простую реализацию на c++ (однопоточную и без учета симметрий), она решила n=8 за 25 минут.
то добавляем к задаче ещё 1 квадратик (всевозможные его варианты) и тем самым подразбиваем задачу примерно на 9 подзадач
Если я правильно понял, то вы перешли от dfs (depth-first search), не очень подходящего для GPU, к bfs (breadth-first search)?
216'285 вариантов заполнения ..., или в 8 раз больше, если включить в ответ все повороты и симметрии
А всех вариантов действительно будет в 8 раз больше? Это верно, только каждый вариант разбиения несимметричен ни в какой симметрии, так как иначе ему будет соответствовать меньше, чем 8 вариантов. Я проверил для n=8 и действительно получил 18656 = 2332*8 вариантов. Интересно, можно ли доказать, что не существует симметричных разбиений?
Можно получить явную формулу для решения. Как уже отмечали в комментариях, есть четыре класса A={0}, B={4,6}, C={1,3,7,9}, D={2,8}, связанных рекуррентным соотношением
Если взять квадрат преобразования, то соответствующий граф разбивается на две части:
Это, кстати, следует из того, что граф переходов получается двудольный, то есть, рёбра есть только между двумя множествами вершин {0,1,3,7,9} и {2,4,6,8} Получившиеся матрицы 2x2 легко диагонализуются (с собственными числами ), и можно получить формулу, аналогичную Фибоначчи, для каждого класса. Например, количество путей из цифры 1:
Так как первое слагаемое меньше единицы, то можно формулу представить в таком виде:
Естесственно, как и для Фибоначчи, эти формулы быстро теряют точность даже для небольших n, так что быстрое возведение соответствующих матриц в степень n по-прежнему самый быстрый способ.
Некий китаец посчитал количество решений до n=59 (oeis.org/A071983), и у меня нет идей, как он в принципе мог это сделать. Там похоже, что количество решений растет как 4^n, поэтому любое перечисление решений будет неэффективно. Тем не менее, простой DFS с оптимизациями (представление ребер графа и подмножеств вершин битовыми масками) на практике оказался почти самым быстрым: 4 часа для n=46 (во всех алгоритмах я считаю однопоточное исполнение). До этого его скорость растет примерно как 3^n, но понятно, что на больших n она выростет как минимум до 4^n.
DP оказалось неэффективно: для n=39 у меня уже не хватает памяти (растет экспоненциально), а решение для n=38 занимает 4.5 минуты.
Самое эффективное решение, которое я нашел, - комбинация DFS и DP: для нечетного n ищем с помощью DFS количество цепочек F(a, S), начинающихся в a, для всех подмножнств S мощности (n+1)/2. Эти позволяет рассчитать количество гамильтоновых путей, имеющих a ровно посередине. Ищет n=47 за 6 часов, до этого растет примерно как (1.7)^n. К сожалению, память растет тоже экспоненциально (но медленнее, чем стандартный DP).
Самый лучший теоретический алгоритм основан на принципе inclusion/exclusion, он имеет сложность (2^n)*(n^3) при полиномиальной памяти, но он требует возведения матриц nxn в степень n, поэтому он на практике слишком медленный (реализация на питоне не доходит даже до n=30). Рано или поздно он обгонит по скорости все прочие алгоритмы, но, скорее всего, время работы будет превышать время существование вселенной. Для этого алгоритма, правда, есть огромное поле оптимизаций (полное распараллеливание на 2^n умножений матриц, которые можно делать на GPU) - но по моим подсчетам, все равно будет слишком медленно, если только не найти эвристик, которые позволят не считать большую часть подмножеств.
А почему у Вас получилась такая оценка? всего вариантодля еще дают , и для построения каждого надо примерно операций (примерное количество ребер из ). То есть, получается И кажется, можно обойтись без - если взять - количество гамильтоновых путей в , начинающихся в .
Ага, правда, иногда случайности не случайны. Известна история, когда два студента Принстона обнаружили, что комбинация их фамилий звучит по-английски очень неприлично, и решили, что должны обязательно написать совместную статью. Теперь в математике есть машина Кокса-Цукера (Cox-Zucker Machine)
Это, конечно, баловство, но я уместил решение в восемь строчек на питоне! Тут смысла не очень много, это просто формальное преобразование ранее полученных формул с целью уменьшения количества кода.
Спасибо, интересно было вспомнить (хотя у мы только симплекс-метод проходили, когда-то вручную его программировал, без всяких библиотек).
это вот этот курс? оставлю ссылку, может, кто-то захочет посмотреть (а вдруг?!!)
Спасибо, интересная задача!
для всех 
, будем называть такие числа достижимыми (из нуля).
- количество достижимых чисел в отрезке ![[0, b^k]](https://habrastorage.org/getpro/habr/formulas/0/00/004/004304a805b694745c153506aacadd11.svg)
и 
, равная единице, если сумма цифр последнего числа в отрезке ![[0, b^k]](https://habrastorage.org/getpro/habr/formulas/0/00/004/004304a805b694745c153506aacadd11.svg)
не превосходит n, иначе нулю. То есть, она говорит, можем ли мы продолжать смотреть следующие отрезки, или надо остановиться.
Я немного подумал над ней, но в комментариях хабра мало места, чтобы подробно написать. Примерно так получилось (рассуждения, в принципе, повторяют ваши, но только в более привычных обозначениях):
Во-первых, возьмем основание системы счисления произвольное b вместо 10, мне так формулы легче воспринимаются. Для положительных x обозначим сумму цифр
Теперь сначала рассмотрим одномерный случай: посчитаем количество таких x, что
Для этого введем две функции:
Тогда f и q можно выразить такими формулами:
Теперь точно так же рассмотрим двумерный случай. Удивительно, что условие на то, надо ли нам рассматривать следующий квадрат или нет, сводится к той же самой функции q (не зависящей от размерности!). Я строго не смог доказать это, но там должны быть примерно ваши рассуждения. Для двумерной f формула почти такая же:
Или можно ее сократить до одной суммы, если взять diag(t) равной количеству пар x, y в квадрате bxb, в сумме не превосходящих t:
Теперь, имея одномерную и двумерную функции, можно получить ответ по принципу inclusion-exclusion.
Что касается кода, то проверить правильность результата можно обычным поиском (в глубину или ширину, неважно): код на питоне
Вычислить сразу все решения до заданного n можно по этой рекуррентной формуле: код на питоне, вроде бы результаты совпадают с вашими c b=10, N=900.
Сложность алгоритма O(n^2 * b), так как надо итерировать f до значений k=O(n).
упс, хотел написать длинный комментарий, но случайно опубликовал его. Напишу его заново и опубликую
В математике вероятность определяется как нормированная мера на сигма-алгебре множеств, в этом контексте событие определяется как измеримое множество.
Поэтому явное перечисление множеств (событий) и задание меры на них является фундемантальным и математически корректным способом ввести вероятность.
Спасибо, интересно. Я всегда считал, что естественное обобщение - произведение n-1 вектора в n-менром пространстве как вектор длиной площадь n-1-мерного параллелограмма и перпендикулярный им, и его можно посчитать как определитель матрицы, в которой строки - координаты векторов и базис.
Но, согласно википедии, существует много вариантов обобщений, в зависимости от того, какие свойства векторного произведения мы хотим сохранить.
А точно есть такие правила? В Википедии (особенно в Computer Science) много ссылок на форумы и даже личные страницы, когда, например, нет или мало научных статей на эту тему - часто встречается, когда люди делают исследования и не считают нужным их публиковать в авторитетных источниках.
Например, в этой статье первые четыре ссылки ведут на форумы и личные страницы (хотя, насколько я знаю, основные результаты по этой теме опубликованы)
Это стандартная задача на поиск в глубину, таких много на литкоде. Самый простой (и медленный) вариант на питоне:
Спасибо, интересная статья!
Насчет сортировки циклических сдвигов - а вы не пробовали стандартный prefix doubling algorithm - он, насколько я помню, довольно просто пишется, и работает за O(N*log N), сильно быстрее наивного, хотя и медленнее линейных. Но и классический алгоритм Укконена по построению суффиксного дерева на питоне занимает всего 50 строк, хотя самому его написать (и, главное, отладить), наверное, сложно.
как это не меняется? (5, 1) -> (6, 4)
не уверен насчет оптимизатора js, но после оптимизатора с++ не будет ни того, ни другого: деление превратится в умножение на reciprocal (throughput = 1 инструкция/такт), а if - в conditional move (2 инструкции/такт). Мне кажется, ни то, ни то не будет узким местом, скорее уж обращение к памяти
я пытался строить все варианты с половиной кубиков (18 для n=8, 22 или 23 для n=9).
В вашем случае границы, как на рис.4, сложность в том, что на любой позиции над средней линией можно поставить квадратик размера 1 (допустим, он бы еще не был использован), и тоже получится правильная граница. То есть, "жадное" заполнение сверху вниз и слева направо не работает, надо перебирать больше вариантов, по крайней мере, близко к границе. Или я чего-то не учитываю?
А сколько у вас времени работало для n=8?
Я немного поэкспериментировал с задачей, но комментарии хабра слишком коротки, чтобы подробно описать результат, поэтому вот очень примерное описание решения задачи для n=8 (n=9 я пока даже не трогал):
Неоптимимизированное решение: стандартный DFS (depth first search with backtracking), состояние представленно как три 16-байтных массива (текущие высоты, длины, и количество оставшихся квадратов). После заполнения текущего нулевого уровня, все заполненные строки схловываются, как в тетрисе, так что заполняется всегда нулевой уровень, что упрощает код.
Время работы: 24 мин, количество шагов DFS: 64G, код
Симметрия: рассматриваются только такие состояния, в которых самый первый квадрат не меньше второго и не меньше последнего в первом ряду, остальные состояния достраиваются с помощью симметрии.
Время: 14 мин, шагов DFS: 35G, код
Эвристика: после добавления квадрата, рассматриваются два предыдущих: если средний квадрат в этой тройке ниже, чем крайние, то проверяется, что у нас еще остались достаточно маленьуие квадраты, чтобы в будущем заполнить этот промежуток.
Время: 12 мин, шагов DFS: 23G, код
Низкоуровневые оптимизации перехода состояний иcпользуют 128-битные __uint128_t и регистры SSE, а также ускоренное нахождение размеров квадратов, которые можно использовать на текущем шаге, через битовые маски
Время: 9.5 мин, шагов DFS: 23G, код
Еще пробовал, но не получилось:
5. Еще более низкоуровневые оптимизации: все операции в 128-бит или SSE, уменьшение размера состояния до 32 байт, уменьшение дорогостоящих операций обновления состояния: ускорение получается меньше 10%, а код очень сильно усложняется.
6. Я пробовал делать DFS в половину глубины, а потом сшивать попарно. К сожалению, идея не работает: обратный DFS не может сгененировать некоторые варианты (например, когда в первой половине есть "яма" длиной 2 и высотой 4)
Спасибо, интересная задача, крутое исследование! Вы, наверно, можете в oeis.org завести соответствующую последовательность (будет что-то типа 1,0,0,0,0,0,0,2332,216285)
Я тоже хочу попробовать ускорить ее на CPU, как будет больше времени)) Тут возможны низкоуровневые оптимизации для представления текущего заполнения, интересно, насколько они будут полезны.
Еще можно попробовать такую оптимизацию: при добавлении нового квадрата, если он получается выше предыдущего, проверять, что существуют свободные квадраты, которые можно поставить на предыдущий квадрат (если пред-предыдущий квадрат выше) - но не факт, что это сильно сократит количество вариантов для перебора.
Я написал простую реализацию на c++ (однопоточную и без учета симметрий), она решила n=8 за 25 минут.
Если я правильно понял, то вы перешли от dfs (depth-first search), не очень подходящего для GPU, к bfs (breadth-first search)?
А всех вариантов действительно будет в 8 раз больше? Это верно, только каждый вариант разбиения несимметричен ни в какой симметрии, так как иначе ему будет соответствовать меньше, чем 8 вариантов. Я проверил для n=8 и действительно получил 18656 = 2332*8 вариантов. Интересно, можно ли доказать, что не существует симметричных разбиений?
Кажется, что это не верно для 2S, а ваша оценка точная. Например, для набора (6, 22, 33, 66):
. Но 
больше 2S, но его нельзя представить с номиналом 66.
Можно получить явную формулу для решения.
Как уже отмечали в комментариях, есть четыре класса A={0}, B={4,6}, C={1,3,7,9}, D={2,8}, связанных рекуррентным соотношением
Если взять квадрат преобразования, то соответствующий граф разбивается на две части:
Это, кстати, следует из того, что граф переходов получается двудольный, то есть, рёбра есть только между двумя множествами вершин {0,1,3,7,9} и {2,4,6,8}
), и можно получить формулу, аналогичную Фибоначчи, для каждого класса. Например, количество путей из цифры 1:
Получившиеся матрицы 2x2 легко диагонализуются (с собственными числами
Так как первое слагаемое меньше единицы, то можно формулу представить в таком виде:
Естесственно, как и для Фибоначчи, эти формулы быстро теряют точность даже для небольших n, так что быстрое возведение соответствующих матриц в степень n по-прежнему самый быстрый способ.
Некий китаец посчитал количество решений до n=59 (oeis.org/A071983), и у меня нет идей, как он в принципе мог это сделать.
Там похоже, что количество решений растет как 4^n, поэтому любое перечисление решений будет неэффективно.
Тем не менее, простой DFS с оптимизациями (представление ребер графа и подмножеств вершин битовыми масками) на практике оказался почти самым быстрым: 4 часа для n=46 (во всех алгоритмах я считаю однопоточное исполнение). До этого его скорость растет примерно как 3^n, но понятно, что на больших n она выростет как минимум до 4^n.
DP оказалось неэффективно: для n=39 у меня уже не хватает памяти (растет экспоненциально), а решение для n=38 занимает 4.5 минуты.
Самое эффективное решение, которое я нашел, - комбинация DFS и DP: для нечетного n ищем с помощью DFS количество цепочек F(a, S), начинающихся в a, для всех подмножнств S мощности (n+1)/2. Эти позволяет рассчитать количество гамильтоновых путей, имеющих a ровно посередине. Ищет n=47 за 6 часов, до этого растет примерно как (1.7)^n. К сожалению, память растет тоже экспоненциально (но медленнее, чем стандартный DP).
Самый лучший теоретический алгоритм основан на принципе inclusion/exclusion, он имеет сложность (2^n)*(n^3) при полиномиальной памяти, но он требует возведения матриц nxn в степень n, поэтому он на практике слишком медленный (реализация на питоне не доходит даже до n=30). Рано или поздно он обгонит по скорости все прочие алгоритмы, но, скорее всего, время работы будет превышать время существование вселенной.
Для этого алгоритма, правда, есть огромное поле оптимизаций (полное распараллеливание на 2^n умножений матриц, которые можно делать на GPU) - но по моим подсчетам, все равно будет слишком медленно, если только не найти эвристик, которые позволят не считать большую часть подмножеств.
А почему у Вас получилась такая оценка? всего
вариантодля 
еще 
дают 
, и для построения каждого 
надо примерно 
операций (примерное количество ребер из 
). То есть, получается 
- если взять 
- количество гамильтоновых путей в 
, начинающихся в 
.
И кажется, можно обойтись без