ripatti Jun 11 at 07:19

Винтик и Шпунтик, часть 2: гиперкубы, шляпы и фартуки

Medium

26 min

952

Кубы, шляпы, фартуки, винтики и (возможно) шпунтики

Это вторая часть моих наработок по решению задачи про Винтика и Шпунтика в рамках челленджа @vvvphoenix. В первой части мы выразили ответ в виде формулы включений-исключений. Хотя в подобных формулах и получается огромное число слагаемых, часто оказывается, что либо они почти все равны нулю, либо объединяются в сравнительно небольшое число групп одинаковых, либо ещё что-нибудь. И в итоге огромная формула вычисляется чуть ли не на бумаге. В этой части мы будем сворачивать и оптимизировать нашу формулу для ускорения вычислений.

Настоятельно рекомендую ознакомиться с предыдущей вводной частью для того, чтобы понимать что вообще происходит в этой части.

Представление задачи в виде гиперкубов

Чтобы лучше видеть все симметрии задачи, рассмотрим следующий способ упорядочивания вершин двудольного графа, полученного в предыдущей части.

Каждому-значному числу без звёздочек правой доли графа поставим в соответствие точку в-мерном пространстве с координатами, соответствующим цифрам в этом числе. Например, для тернарного случая, всевозможные трёхзначные числа с цифрами от 0 до 2 полностью заполнят трёхмерный куб $3 \times 3 \times 3$ .

Для чисел со звёздочками также поставим в соответствие точки в трёхмерном пространстве, ставя каждой звёздочке координату. Тогда для тернарного случая вершины левой доли графа будут расположены в трёх квадратах размера $3 \times 3$ , которые нависают около куба с трёх различных сторон.

Две вершины из разных долей будут соединены ребром тогда и только тогда, когда они находятся на одной линии, параллельной одной из осей координат.

Для кватернарного случая у нас уже будет гиперкуб $4 \times 4 \times 4 \times 4$ , прикрытый с четырёх сторон четырьмя кубами размера $4 \times 4 \times 4$ .

Как теперь считать слагаемые нашей формулы включений-исключений? Для какого-то фиксированного множестваточек внутри куба поставим в соответствие, а для всех остальных — . Теперь для каждой точки в "сателлитах" посчитаем сумму чисел внутри куба на линии, которая проходит через данную точку. В конце полученные числа перемножим. Это будет значение функции f(X) . Картинка ниже иллюстрирует процесс вычисления.

Тут 13 двоек и 4 тройки, f(X) равно их произведению

Однако для такой формулировки я не придумал эффективного метода сворачивания формулы, поэтому ты сделаем следующий финт: мы поменяем доли двудольного графа местами. Понятно, что от этого количество совершенных паросочетаний не изменится.

Тогда множествобудет задавать какие-то вершины из висящих квадратов. Опять поставим в соответствие таким вершинам, а всем остальным —. Для каждой вершины внутри куба найдём сумму по всем соответствующим вершинам вне куба. Соответствующие вне куба — это все такие, у которых все координаты кроме одной (где звёздочка!) совпадают с координатами точки внутри куба. Другими словами, мы из точки внутри куба "осматриваемся" по всемнаправлениям, параллельным осям координат, и суммируем числа по всем внешним точкам, которые увидим. Теперь f(X) — это произведение по всем числам внутри куба. Иллюстрация процесса ниже.

Значение f(X) тут то же, что и для картинки выше

Сворачивание формулы

Для простоты начнём с тернарного случая. Назовём один из "парящих" квадратов (скажем, верхний) шляпой (hat) потому что этот квадрат прикрывает наш куб $3 \times 3 \times 3$ сверху. Остальные квадраты, прикрывающие куб с боков назовём фартуками (apron). Итак, наш куб прикрыт одной шляпой и двумя фартуками. В понятия шляп и фартуков мы также добавим значения 0 и 1, приписываемые точкам внутри соответствующих квадратов, то есть по сути шляпы и фартуки — это двоичные маски длины $3 \times 3 = 9$ .

Пусть— множество извсевозможных шляп, а— множество из 2^9 всевозможных фартуков. Тогда ответ считается по формуле

$\sum_{h \in H} \sum_{a \in A} \sum_{a' \in A} s(h,a,a') F(h,a,a') = \sum_{h \in H} s(h) \sum_{a \in A} \sum_{a' \in A} s(a,a') F(h,a,a').$

Здесь функциязадаёт знакилив зависимости от чётности количества нулей во всех подаваемых на вход масках, а функциявозвращает значение f(X) , сначала преобразовав все подаваемые на вход маски в подходящее множество. Эта небольшая вольность в обозначениях нужна, чтобы выражения получались не трёхэтажными (а хотя бы двухэтажными).

Весь куб 3х3х3 и он же, разрезанный на 3 слоя

Теперь разрежем весь наш куб на слои, параллельные шляпе, а все фартуки разрежем на лоскуты $1 \times 3$ . Пусть A^* — множество из 2^3 всевозможных лоскутов. Тогда формула примет вид

$\sum_{h \in H} s(h) \sum_{a_1 \in A^*} \sum_{a'_1 \in A^*} \sum_{a_2 \in A^*} \sum_{a'_2 \in A^*} \sum_{a_3 \in A^*} \sum_{a'_3 \in A^*} s(a_1,a'_1)F(h,a_1,a'_1) \times$ $\times s(a_2,a'_2)F(h,a_2,a'_2)s(a_3,a'_3)F(h,a_3,a'_3),$

что довольно несложно сворачивается в

$\sum_{h \in H} s(h) \left( \sum_{a \in A^*} \sum_{a' \in A^*} s(a,a') F(h,a,a') \right)^3 = \sum_{h \in H} s(h) \left( \sum_{a \in A^*} s(a) \sum_{a' \in A^*} s(a') F(h,a,a') \right)^3.$

Как оно так свернулось? Ну раскройте скобку обратно и посмотрите. Тут можно сделать следующее наблюдение: как только мы зафиксировали шляпу, далее слои можно обрабатывать независимо друг от друга — именно этим мы и пользуемся.

Продолжим шинковать маски на куски. Разрежемна три части, а также разрежем на три части. Формула преобразуется в следующую:

$\sum_{h_1 \in H^*} \sum_{h_2 \in H^*} \sum_{h_3 \in H^*} s(h_1,h_2,h_3) \times$ $\times \left( \sum_{a \in A^*} s(a) \sum_{a'_1 \in A^{**}} \sum_{a'_2 \in A^{**}} \sum_{a'_3 \in A^{**}} s(a'_1) F(h_1,a,a'_1) s(a'_2) F(h_2,a,a'_2) s(a'_3) F(h_3,a,a'_3) \right)^3,$

что совершенно аналогично тому, что мы делали выше, сворачивается в

$\sum_{h_1 \in H^*} \sum_{h_2 \in H^*} \sum_{h_3 \in H^*} s(h_1,h_2,h_3) \times$ $\times \left( \sum_{a \in A^{*}} s(a) \left( \sum_{a' \in A^{**}} s(a')F(h_1,a,a') \right) \right. \times$ $\times \left. \left( \sum_{a' \in A^{**}} s(a')F(h_2,a,a') \right) \left( \sum_{a' \in A^{**}} s(a')F(h_3,a,a') \right) \right)^3.$

Ещё одного вложенного куба не получается из-за того, что каждая из самых внутренних скобок зависит от своего куска шляпы h_1 , h_2 и h_3 . Каждая из самых внутренних скобок зависит от масок h_i и, и эти значения мы можем заранее предподсчитать! Тогда итоговая формула примет вид

$\sum_{h_1 \in H^*} \sum_{h_2 \in H^*} \sum_{h_3 \in H^*} s(h_1,h_2,h_3) \left( \sum_{a \in A^*} s(a) dp(h_1,a)dp(h_2,a)dp(h_3,a) \right)^3,$

где

$dp(h,a) = \sum_{a' \in A^{**}} s(a') F(h,a,a').$

В общем случае, когда тернарная задача не $3 \times 3 \times 3$ , а $n \times m \times k$ , формула примет вид

$\sum_{(h_1, \cdots , h_m) \in H} s(h_1 \cdots,h_m) \left( \sum_{a \in A^*} s(a) \prod_{i \in [1,m]} dp(h_i,a) \right)^k.$

Теперь оценим сложность вычисления формулы для тернарной задачи (n,m,k) . Все значенияможно вычислить за $O^*(2^{2n+1})$ , после чего сама формула может быть вычислена за $O^*(2^{mn+n})$ , что покрывает все затраты на предподсчет. Нам, кстати, выгодно считать, что $n \le m \le k$ — так вычисления идут быстрее.

Решаем тернарный случай быстрее

Наши старые знакомые — задачи (3,3,3), (2,4,4) и (2,3,6)

С использованием формулы выше, напишем код, решающий тернарный случай, на языке Python (обычно я пишу на C++, но пока мы делаем прототип решения, и мне не хочется возиться с длинной арифметикой):

def calc_dp( m1, m2, n ):
	x,y = 1,1
	for i in range(0,n):
		c = ((m1>>i)&1) + ((m2>>i)&1)
		x *= c+1
		y *= c
	return x-y

def sign( x, L ):
	if (L-x.bit_count())&1:
		return -1
	return 1

def solve3d( A, B, C ):
	print( (A, B, C) )
	dp = [ [calc_dp( i, j, A ) for j in range(0,1<<A)] for i in range(0,1<<A)]
	ans = 0
	for hat in range(1<<(A*B)):
		sum = 0
		for apron in range(1<<A):
			prod = 1
			for part in range(B):
				prod *= dp[(hat>>(part*A))&((1<<A)-1)][apron]
			sum += sign(apron, A) * prod
		ans += sign(hat, A*B) * (sum ** C)
	print( ans )

solve3d( 3, 3, 3 )
solve3d( 2, 4, 4 )
solve3d( 2, 3, 6 )

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

Свёрнутая формула для кватернарного случая

Для кватернарного случая сворачивание происходит аналогичным образом, в итоге для задачи (n, m, k, t) получается формула

$\sum_{(h_1, \cdots , h_k) \in H} s(h_1 \cdots,h_k) \left( \sum_{a \in A^*} s(a) \prod_{i \in [1,k]} dp(h_i,a) \right)^t,$

где все так же

$dp(h,a) = \sum_{a' \in A^{**}} \sum_{a'' \in A^{***}} s(a',a'') F(h,a,a',a''),$

только h_i иимеют размерность $n \times m$ ,— размерностьи a'' — размерность. Следующая картинка позволяет взглянуть на происходящее с геометрической точки зрения.

У Незнайки тоже была большая синяя шляпа

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

Предподсчёт всехтеперь выполняется за $O^*(2^{2mn+m+n})$ , однако его можно ускорить до $O^*(2^{2mn+n})$ , применив тот трюк, что мы использовали два раза, еще один третий раз. В любом случае, дальше у нас идёт более ресурсоёмкое вычисление самой формулы за время $O^*(2^{nmk + nm})$ . Это значительное ускорение относительно вычисления "в лоб". Например, для кватернарного случая (4,4,4,4) количество операций с порядка $2^{256}$ уменьшается до порядка $2^{80}$ . К сожалению, это всё ещё медленно для вычислений на практике. Однако уже достаточно быстро, чтобы посчитать за разумное время ответы для расширенной версии задачи.

Вспомним про египетские дроби

Тут самое время вспомнить про три "школьные" задачи из первой части. Активности в комментариях они что-то не прибавили, ну да ладно, изложу здесь их решение.

Задача. Докажите, что при каждомуравнение

$\dfrac{1}{x_1} + \dfrac{1}{x_2} + \ldots + \dfrac{1}{x_n} = 1$

в натуральных числах имеет конечное число решений.

Доказательство. Докажем более общее утверждение: для любого $0 < a \le 1$ при каждомуравнение

$\dfrac{1}{x_1} + \dfrac{1}{x_2} + \ldots + \dfrac{1}{x_n} = a$

в натуральных числах имеет конечное число решений. Доказывать будем индукцией по.

База индукции: для n=1 для любогоочевидно, что решений либо одно, либо ни одного — в зависимости от того, являетсяегипетской дробью или нет.

Теперь предположим, что для n=k и для любого $0 < a \le 1$ утверждение верно. Рассмотрим случай n = k+1 для произвольного $0 < a \le 1$ . Пусть 1 / x_1 — самая большая дробь в сумме $1/x_1 + \cdots + 1/x_{k+1} = a$ . Тогда $x_1 \le (k+1)/a$ , поскольку иначе

$\dfrac{1}{x_1} + \cdots + \dfrac{1}{x_{k+1}} \le \dfrac{k+1}{x_1} < \dfrac{k+1}{(k+1)/a} = a.$

А так как x_1 число натуральное, оно может принимать лишь конечное число значений — не более $\lfloor (k+1)/a \rfloor$ вариантов. Переберём все эти варианты, и для каждого из них остаток суммы дробей $1/x_2 + \cdots + 1/x_{k+1} = a' = a - 1/x_1$ даёт конечное число решений согласно предположению индукции для a' > 0 (для $a' \le 0$ их тоже конечное число — ровно нисколько). Итак, утверждение доказано.

Наше исходное утверждение является следствием только что доказанного для a=1 .

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

Решение. Мой код, в принципе, в точности повторяет доказательство из задачи 1:

from fractions import Fraction

def f(n, a, last, arr):
	if n==0:
		if a==Fraction(0,1):
			print( arr, arr[0]*arr[1]*arr[2] + arr[0]*arr[1] )
		return
	if a==Fraction(0,1):
		return
	while Fraction(1,last)>a:
		last += 1
	while Fraction(n,last)>=a:
		arr.append(last)
		f(n-1, a-Fraction(1,last), last, arr)
		arr.pop()
		last += 1

f(4, Fraction(1,1), 1, [])

Единственное что — я дополнительно храню последнюю дробь в переменной last для того, чтобы генерировать дроби в решении в порядке убывания и избегать повторяющихся решений, отличающихся только порядком дробей. Ну и еще иду рекурсивно до n=0 , а не n=1 .

Задача. Найдите все решения для n=4 .

Запустим код выше и получим:

[2, 3, 7, 42] 48
[2, 3, 8, 24] 54
[2, 3, 9, 18] 60
[2, 3, 10, 15] 66
[2, 3, 12, 12] 78
[2, 4, 5, 20] 48
[2, 4, 6, 12] 56
[2, 4, 8, 8] 72
[2, 5, 5, 10] 60
[2, 6, 6, 6] 84
[3, 3, 4, 12] 45
[3, 3, 6, 6] 63
[3, 4, 4, 6] 60
[4, 4, 4, 4] 80

Кстати, код можно позапускать и с другими параметрами. Он выдает корректные ответ; во всяком случае для $n \le 6$ количество найденных решений совпадает с табличкой на OEIS (для n=7 я не дождался завершения работы программы).

Рядом с задачами я вывел значения nmk + nm , которые косвенно указывают на сложность решения задачи. Задача (4,4,4,4) одна из самых сложных, если не самая сложная (в задаче (2,6,6,6) показатель сложности больше, однако на самом деле её решить проще с той идеей, что будет описана в третьей части).

Кстати, для задачи (2,3,7,42) свёрнутая формула ускоряет решение с порядка $2^{1764}$ до порядка $2^{48}$ операций.

Зачем вообще нужны дополнительные задачи? Во-первых, некоторые из них проще исходной — приятно получить хоть какой-то прогресс в решении задачи. Во-вторых, при помощи них можно тестировать более усложнённое (и ускоренное) решение. Ну и в-третьих, вдруг на основе ответов на все эти задачи удастся вывести более простую формулу.

Код для кватернарного случая и результаты вычислений

Код на Python довольно компактный:

import sys
from timeit import default_timer as timer

def sign( x, L ):
	if (L-x.bit_count())&1:
		return -1
	return 1

def calc_dp2t( M, n, m ):
	trits = []
	for i in range(n*m):
		trits.append( M%3 )
		M //= 3
	res = 0
	for x in range(1<<n):
		for y in range(1<<m):
			prod = 1
			for i in range(n):
				for j in range(m):
					prod *= trits[i*m+j] + ((x>>i)&1) + ((y>>j)&1)
			res += sign(x, n) * sign(y, m) * prod
	return res

def solve4d( A, B, C, D ):
	start = timer()
	print( (A, B, C, D) )
	to_ternary = []
	for i in range(1<<(A*B)):
		tmp = 0
		for j in range(A*B):
			if ((i>>j)&1):
				tmp += 3 ** j
		to_ternary.append( tmp )
	dp = [ calc_dp2t( i, A, B ) for i in range(3**(A*B)) ]
	ans = 0
	subtask_bit = 0
	last_ans = 0
	for hat in range((1<<(A*B*C))-1,-1,-1):
		if subtask_bit<A*B*C and ((hat>>subtask_bit)&1)==0:
			print( "subtotal", subtask_bit, ans, file=sys.stderr )
			print( "subtotal", subtask_bit, ans, last_ans/ans, flush=True )
			last_ans = ans
			subtask_bit += 1
		sum = 0
		for apron in range(1<<(A*B)):
			prod = 1
			for part in range(C):
				h = to_ternary[(hat>>(part*A*B))&((1<<(A*B))-1)]
				prod *= dp[h+to_ternary[apron]]
			sum += sign(apron, A*B) * prod
		ans += sign(hat, A*B*C) * (sum ** D)
	print( "total", ans )
	end = timer()
	print( "time", end - start )
	print()

solve4d( 3, 3, 4, 12 )

Основное отличие от описанного выше заключается в вычислении таблицы. А проблема следующая: для задачи (4,4,4,4) эта таблица получается размера $2^{16} \times 2^{16}$ и просто не лезет в память. Однако заметим, что маскиипри вычислениив точности накладываются друг на друга и дальнейшее вычисление зависит только от суммарного количества единичных бит на каждой изпозиций. Таким образом, размерможно сократить с $2^{16} \times 2^{16}$ до $3^{16}$ , то есть задавать троичной маской. Это и сокращает потребление памяти, и ускоряет просчёт всей таблицы.

Я запустил этот код на нескольких самых "простых" задачах, а также на (4,4,4,4) . Вычисления производились на PyPy 3.11 (с ускоренным bit_count()) около суток. Результаты получились следующие:

(4,4,4,4)

(4, 4, 4, 4)
subtotal 0 1314773196405160880253287002591853900553104513882342157021949988456542744598543879344935129170374402775187456 0.0
subtotal 1 442578222306575980794016193452499497834563043715036454978110213681371679770318624250355469000358330797916160 2.9707137182507175
subtotal 2 145519359244090994821878476189914088529577478614208054102501888824029041699742587095023660476742483239239680 3.0413700596647413
subtotal 3 46512853090570596000943855309347174795369932897243079808813337517470179087102102297227506406675622794362880 3.128583812322441
subtotal 4 14354140362192317828379230568413897400274511153588249769957155544026944627372158381840267478938173162127360 3.2403788674856355
subtotal 5 4601937725722353260445973563227983189885432816105575724822420229743310277789037895560637129338133365653504 3.119151370076219
subtotal 6 1439066159980865891480356264822965800437568110158710220229209756991752507597317041103103763125063672397824 3.197863902090187
subtotal 7 436752264463226355391056142006284109952490771369913706603959975992791188626338005131947776458154605281280 3.2949254693607486
subtotal 8 127732467142035448254028626320891774836584040179595131226742505104242133750370508724639076125698022178816 3.419273691610202
subtotal 9 38413342307987398919037399189522881145559744170146108393147556873513227522007708546792694113909029732352 3.3252109675308223
subtotal 10 11241896112452093792559291249846653443789972847217375979988440920162954072258853569618900890911644844032 3.4169807231574425
subtotal 11 3184581691344181290344497454493521356653434573929489580215587915936266471388015829270695865858469658624 3.5301013451807535
subtotal 12 866542225166339187020985080909987461405932268889877857338595776048029368081683993640148219016582791168 3.675045022454477
subtotal 13 238471947764164511594884055711277089648018987887292234190239171402401242840504757365611556812605095936 3.6337281315087897
subtotal 14 63603015904651276339046961270309227371126167628303991155673766563656049819429600986070600896770736128 3.7493811318894568
subtotal 15 16340710422285505275778717982342111281369973806811662503664613525367279096577024222454326070146498560 3.8923042059364388
subtotal 16 4008891591913203507073631466741124143704490028170337313821633509089390985326153456067893883530379264 4.07611681374678
subtotal 17 1188655294239612282392602020390129007960329931370239245384837429685360827198683198631102922188914688 3.3726275492490094
subtotal 18 343021954846841514372203437964426352545951755554012716147489678452065368741347093609766422629580800 3.4652455256700523
subtotal 19 95842451931612500254931899842840993489234854968238417493385094259826936003581479506934809151668224 3.579018983066102
subtotal 20 25735114315955958662798960154939468113658101841026532366077713106047043559796787186762361776111616 3.7241898658359363
subtotal 21 7183641338420254735717505639152002471681348576729082062273955422756811264189312179573989817974784 3.5824609141212127

(3,3,4,12)

(3, 3, 4, 12)
subtotal 0 38128920463430842687180769101716401194920564337918864447732456633464924874048164343977054301871727968256000000000000 0.0
subtotal 1 12695330519656578429973042615239234440382708577665467096194713884185092173734919572512086162959075490019360745979904 3.0033814719825247
subtotal 2 4119366959621555640671105274104128221920058314399545429553778858563541519217708061850796719554954918047796548337664 3.081864433078546
subtotal 3 1300534236476014958645969701565751138723457960401767141860805665515407852203653114319337394554712855705445805850624 3.1674421511451896
subtotal 4 403810502416207576046779643585609393545936687195251488329753925623311209100143309857231473581537458105874198298624 3.2206548088626827
subtotal 5 121803884236538714466031760089841035372503730127676244746826016191911174972256715176993658026011428032375406198784 3.3152514383861704
subtotal 6 35626045790860489543985620154939862264152808775690077183527997032004532455217077229583192094998595633414798311424 3.4189560343457006
subtotal 7 10202872650190086143517137649128368003577668848541371195912188251939584003822575963887192342436432756687757115392 3.491766192945352
subtotal 8 2827505901549377348766126151973706483130476788141878026185344772436759300770833526619085017142470244701139107840 3.608435492424351
subtotal 9 756588719290948825761231700071041344178644220603488238287060399804814047746368471760656178654701190710770532352 3.73717692248865
subtotal 10 211226686550172062308384752863770351583399618316878402567435312038645251135755374587537738995482462998251962368 3.581880356349947
subtotal 11 56989046231884990239604143240841063743678752253476148412611613042034593187335328535674168406307263917592674304 3.7064436153344875
subtotal 12 14825297049095414416370126185570244182733692369522412728630514860882656273431085020463035254793966679043866624 3.844040766479094
subtotal 13 3758427351050574465224427031116962841434567911035813837563352562014730896513005848588638804346908536828067840 3.944547988921849
subtotal 14 915902064409235235602622638109136332933557054164306342191802454220527366323118325939001026770002370573107200 4.103525362697804
subtotal 15 213954612377425913258863965186927980256014515334206127479406775348101284349560497483284458077490488803328000 4.280824116067857
subtotal 16 48360804100252451394360554049204587106236258796538073954585334606543853767865978681785403686483761226055680 4.4241326495294775
subtotal 17 10434359713779055913857120405584963047192598530202882342564395189798186506456099691726832265017887545098240 4.6347648947150795
subtotal 18 2141420838952314883166136427023434677194155614128079960931547760309559049782646899248314505428855265689600 4.8726338718568005
subtotal 19 451687142172187710027002606587586510146532962456055982187769897518432517048695796539591851059711343329280 4.740938226964149
subtotal 20 90433037103913851693848579356169982472089378189050100816717831729237637974649178855394218510971935129600 4.994713841725428
subtotal 21 17112396347831338317891580280617398723197619579142710216080382068208166152247640091437766899004258385920 5.2846506863063905
subtotal 22 3089154956318234637180138943067041140411441743631178992415694973387975846105010125195194138985045688320 5.539507272961957
subtotal 23 522870455958882265448961769405095705854776826369694976926081932755572847665515739662290493618572492800 5.908069429268272
subtotal 24 82479293311030171951747788128428924422006756823217876017881180430804934782464643407023842586253066240 6.339414839396514
subtotal 25 12204241522718380533547558202017821582068676614319827009590180213816116113528949604878268221515366400 6.758248200635305
subtotal 26 1661729675167003157777926790520862082855835569371796160433867211199281785463322171354071319969792000 7.344300162113829
subtotal 27 206257941891295809539926207265773655024671250168531941534168524877454663830914721301349699000729600 8.05656092526505
subtotal 28 25121447846382979104630940764407150121853816533832083675565146782569688975337883368554440635187200 8.210432103776737
subtotal 29 2738762026468572563763085918483014105378510412257662267000683378355546448966475512251433772646400 9.172555922565932
subtotal 30 262860522724664498296479203316620887287709055621293794995058225663557861060970585202516348108800 10.419069391174087
subtotal 31 22080240343214095353292253901607974523434093574673412921505880036812201854353761609771175116800 11.904785393581516

(2,3,7,42)

(2, 3, 7, 42)
subtotal 0 251613614336283736633564420658829243587868063698968184853079328038721672756912461453486605274141773427332128521987531156493577594931801489082231083176513761059156126951280533914430261812947239346405542748749824000000000000000000000000000000000000000000000000000000000000000000000000000000000000 0.0
subtotal 1 83274969815433379152526994636748945338151136413070736460626237490086619232061411279259698985750283254515458450536605243610354380025405695555442266873406836694628562112473960309090375379210303146455016073888626929704597893492931809209376020702433378304000000000000000000000000000000000000000000 3.021479502111475
subtotal 2 27124482125483099612113947720214097463778546564307164310899219763333390510233650234684542337766875826826528823124430485141189612005476241070229567137099026603997522304169075126407671994000042627007533459648788792402052790066593898772797653857402355712000000000000000000000000000000000000000000 3.0701035850265184
subtotal 3 8690533313808429080297585594356865033348428648802925898102853606945999039279871301604118871122774220734971587106093324384503246357279024701581530811578628392653029002896869071594727216731933481497189692554686759567531464800591351722342440525574963200000000000000000000000000000000000000000000 3.121152770035975
subtotal 4 2737336640454619625780580366624904925399588020920473994024478844595746221505274468150896828897342256895620695457021771014397243277254612143987613863970833382293085705112278686956381309379869854024955027548347831729930133807834358551941326402958131200000000000000000000000000000000000000000000 3.174813497679663
subtotal 5 847133834883668454842878615202605883470044923548931522851621203347916849820443895360532942137537021758433377927004955806873618229766658716268371488534630100138813671194186351889739236476164339358675187989448838471874454234538217217378661820897689600000000000000000000000000000000000000000000 3.2312918310369705
subtotal 6 257423632088931244676019686727514973909038981403731333730688528303427583146309148686549359433237749611906954396629149383762064732655715431700568139590627099670578766274418726680075762811068725493321683184937731083845465135971018017582372310797516800000000000000000000000000000000000000000000 3.290816107322315
subtotal 7 76759470210880337075831226826819576876825607578126573437265677758188902867777896541764619797058424021780957154352967528506244914629283360992312713070550167912757242724324474335060541848221712432563663330906806062555306390393479807396931719240089600000000000000000000000000000000000000000000 3.3536400314086916
subtotal 8 22443985690604910618346764516465916129577437478187576769639973111381714327187005924456154723836228078010153738258126089131331408653656325032267857695673814712897146355622896499767102056387045778024828442807028310716268752503725193550009575512473600000000000000000000000000000000000000000000 3.4200462996646794
subtotal 9 6430294996036501772175898907266616063655521632666661863157465605108568672535351754954170062269241372188736499847240394054429253881087684801760909023073603190246573566136754039650413345206059678229890540769839633407772287078351783628346880530841600000000000000000000000000000000000000000000 3.4903508632868183
subtotal 10 1803775877515573089106153809290739310393452356192401153781861649037350663429498796545761456746607195827763386243735200998563156189710695795283606136662320890475870976680702575753352318543149983783301655002740935348242763389974192569539703812915200000000000000000000000000000000000000000000 3.5649079667775885
subtotal 11 494983094008962139600724725587901225609662007395670362145294743874587386500771050525669424348142790801137685374122187275032274782195785202408392006504589256582612897137374301530332056491093639956042796174045817677520444816079942759740505010995200000000000000000000000000000000000000000000 3.6441161311317716
subtotal 12 132759289737102900419901463558624291087843689441800793905789016614520407143482682995204504601630467555594954210782200796807709553077366376376401974579626937702923205405970771504289312083940784985555210243831504989108751163556910786093864137523200000000000000000000000000000000000000000000 3.7284252950520775
subtotal 13 34768800716530213655616107912209120795102808412819553387062914345846738833037358149592145753995749752110056926099602895692168983946318473595575919878987865941425219684255137302511988520903440102991585002519643986558790622735497357254422901555200000000000000000000000000000000000000000000 3.81834538440047
subtotal 14 8882152437527763775186282589458867482055676785903003850356204213996927768620694912018929847894789801757985379933424418113272190101312730135368249195598577318850538235255564800404481660186192887121747098795560489558388371661935034190475244339200000000000000000000000000000000000000000000 3.914456654631304
subtotal 15 2210908360140887553957191451818550781692413697947491349505096267380032819948680148599150769966077054974408358625025529822806587723408535990944783447335561585203047873917051766876404633700728797053678406103371606956467923354414001072950476800000000000000000000000000000000000000000000000 4.0174222494512435
subtotal 16 535587804906189454536838790103645709167518972367494778941890829955581397735877384839145278632477374587245341520944040975043093934087832914742466598737502759610718798576371991250840006327200812987400792372281005859480367136726439821974700032000000000000000000000000000000000000000000000 4.128003550282736
subtotal 17 126107330789801647479029320750155638412037848457414460393345579764511317758638040945513322014903338379930693065031397241081606300206062399371099037213192357409158987710107654481353318627903799857971399691429588856075709867767731570058723328000000000000000000000000000000000000000000000 4.247079067900648
subtotal 18 28820133198547531467623065582531885440061708513412024167618757424373209622832817606796799040287786427630987817200702488504203833113922089715879577141152668490064325669667913248214654574663513775860085900213743457485560431890477786817626112000000000000000000000000000000000000000000000 4.375667868049866
subtotal 19 6383254894227453216780274211401977528713672096753564626214781683845284274809805081985876929356265807947523926103053550905036074535465549507799423983860230948876008016384712915137824695212033990458550932379529383063473098117573989367808000000000000000000000000000000000000000000000000 4.514958853454269
subtotal 20 1367933837270383198223755808157255474363044434959396928077131852434192276132891115648774687721418759116790676060440896394962902996458508227920319646906521450880895805158897813494548290542626467575062855128628492828017108028648161961574400000000000000000000000000000000000000000000000 4.666347684596204
subtotal 21 283129138469959664573362436579587945017880106458271894325876526025900530864316567287378887435532324206527326611733894540812999544132467026726310015470151089466541679336071465783757184253716261708823200936025016297703967589064547825090560000000000000000000000000000000000000000000000 4.831483769783457
subtotal 22 56486524232540074867571372428688133996710078280265315465566750949730239912238540153450804394568916942847890243732065214966843367562027271451965479750399826113851428277297519649914420258890291158040620572444923163028969973681395681198080000000000000000000000000000000000000000000000 5.012330680931826
subtotal 23 10839353657082662694508492144087618368838056386297295771429067090343014958715537384012002902303632381212824282475082324045124454150893002096228180484803661019375366827751189948468708388399160545958964113753350188122003332190093443072000000000000000000000000000000000000000000000000 5.211244694062601
subtotal 24 1995801458487002548857889270218218466095217976846018674389376660541958880164015503260779900223334992692067809455854416299762691678527009468791434745181173354826155609406302688382367599406768557085369080232208544553463028257696776192000000000000000000000000000000000000000000000000 5.431078132040183
subtotal 25 351663423210900961405078403162113526292465236721756270052846646245012168625188453518996869002210166134744019497193192221411504034250908343576966077550967720276889753170196963073956036288682469423477934847721281626948765842459852800000000000000000000000000000000000000000000000000 5.675317154863936
subtotal 26 59120304004939298458685885818820995915418564863798103278165276330530180402674555555567233022987337479208899930654475915165094019858101375034093227069107994315745117941533291412800973454111654897704633883034114307535762352047652864000000000000000000000000000000000000000000000000 5.948268181799618
subtotal 27 9451212506497010637177982936191282993807411810139025861682295722243119754544256664615927973797651789547277813041478475408589396645690494851559475022885606222888728861046191003872371679594127077445956505182455079179082207099617280000000000000000000000000000000000000000000000000 6.255314221778259
subtotal 28 1431291944643660196275157528070521876606353404628938622242360419020542391943044944997397161664973011328336530021166821588866353189088796250431377238469167015727732806319683316190537770291200676495894325421723624753558568105410560000000000000000000000000000000000000000000000000 6.6032737359183695
subtotal 29 204443603808864164920443549815653972549217810815962628917157788835228724272134962695381354968802234387111681610539276909170871584672754002208111539901663608430579675489761112141250981919515843570149547728822512850690667100241920000000000000000000000000000000000000000000000000 7.000913298230575
subtotal 30 27406422459469679619187164092177362638781316855079734713831154376833409898879595207043015500639097664143764465968791592364414077972171719256304231718303953823537031271164544491125223167973962295000093369092911783599706241761280000000000000000000000000000000000000000000000000 7.4596968689075736
subtotal 31 3427983988994699044731715468592124806239917529876889716362612263324334289208825854997131828700598626735077159356949349380762469901049751125530785712490028329571284360082690137426924564278574369015560420114410210745650207784960000000000000000000000000000000000000000000000000 7.994909704204007
subtotal 32 397337071805914744447593811671107120356556126514444911047676265300468523343503982740742731686642832034365896408843783557632642532660532519761336650781412252594584496831858473291579684096457034778870470583338518873978867548160000000000000000000000000000000000000000000000000 8.62739530800476
subtotal 33 42331441470476921958894636810097795670168162157545819343452161694516190952067233370017621404798253245801377788896687337052948023680228120902626436942620688768086585972769450454059734379911419992116900349791681845214429839360000000000000000000000000000000000000000000000000 9.38633455425864
subtotal 34 4104318791261777514531938994860414276817188579723096836712993974638642000413845857618555930430244773534138410901476114018666244740157652408827121687803818949418985948670698474080812851727131216012143422122330232722932367360000000000000000000000000000000000000000000000000 10.313877557611235

(2,4,5,20)

(2, 4, 5, 20)
subtotal 0 4291877245884182528448370656483600820900737231405265560604980944525705952268697555390350510374637634572786006078047474756998906471644749338464852377600000000000000000000 0.0
subtotal 1 1423991526400441356739117303412985904339282029947659842596738954346663303609754421816958450182345014856260655373126500148003736695890229641808445440000000000000000000000 3.013976674940734
subtotal 2 463398325871857777591184354943570392429592050863447808824270789739260160069952615252563692644368866605133234215118219628304684680135709945152038502400000000000000000000 3.072931961334305
subtotal 3 147805894195757223675346785621794410715411384522698242825717107792088936323253920059172730445998270887437661359264152563642048664381190248495631564800000000000000000000 3.1351816407140256
subtotal 4 46174961112758753533709815563733794647390680595905221816183458623371304845125845963897565885705529738497423148773575041594493004195965648098872524800000000000000000000 3.200996614481533
subtotal 5 13901048730675749758812444550894339016612874325372307021657047546279806924191151875904864832190753020996245557956522644439063074680853173467348992000000000000000000000 3.321689032775164
subtotal 6 4090333093262662669210568392543139441661738042826191628170125779338281685812545586245387557526979265845782301336817989052127363160399399245735526400000000000000000000 3.398512642790064
subtotal 7 1175345764562915631062258852349716273368430243739028656338644069137905649446098700587290889361052117917736558845995780343619906772170793038957772800000000000000000000 3.480110463310143
subtotal 8 329513009727171222224219916447326602271712980477432901968812859992063639465708274533998267948048180494384839968684243201481341590552523294375936000000000000000000000 3.5669176325877796
subtotal 9 94911895851721305688572261484725600271515329370982186575537063091827706239202407226221772961110889397820517492910029893224046073296883497507886660877168441997393920 3.471777765791992
subtotal 10 26661968036085441433225494202404051688556918310639006422930245365107414227785652225050622988597344843299927068806639116073966625457011558312920369727299020222627840 3.5598233304932143
subtotal 11 7296698801915019296118445797571159872277575023732313547196030099978912098505062121216924241722120332938882054988365215610870612148737920384654680358055476818083840 3.6539767859251646
subtotal 12 1943259542737674394280666626364320323505260847431255390952825356868698276052731455848989337010041321889458942221052963496325058796156929074471686124009361275617280 3.7548760942325754
subtotal 13 495014683996298828132828332908200704095345330078539330815277725805715968543794314067006304824962440560430555747282680072205055403798817213758434616840953255690240 3.925660400716929
subtotal 14 122329398087004490615142036515274818743968537656508147742780814729002722681373443692242337983638477416802431163678785579927594062863803410888675481424018454609920 4.046571729587265
subtotal 15 29285536696163668545990412236918650674997161597915588403491044936051763545303639212896640567824058401761658185965435790186145661308434395413095640225266983239680 4.177126728328982
subtotal 16 6781434976931396804064293724251489644734490204293339633827627113472979692910712435711804235222660587492910685567597319163175923582343440796891598578962493603840 4.31848669135443
subtotal 17 1601629893989855100956579284006208313421012890318948033152882249287143241402500020212207454695195227476213975848337951157840215288603224527580883344695912038400 4.23408366837985
subtotal 18 365281845564043035365532185879119290676831706706309444280525579765751224621732850440447522944661787618830549596994953901731208456524497228806411590711430348800 4.384641376076954
subtotal 19 80295824071222765643493019977801258170663749512987359737858830709263296880022090217379396385587221297470958423878007241161156155346504096012685537517240320000 4.549201029931972
subtotal 20 16976676004323944280141480210105449300333319834961073365001378916225786868169180891767255638606905420350043997450505220310729684521419457650476427179668275200 4.729773016270764
subtotal 21 3387016620657016657524028779347763366587812582228443744863222755511661665150243381770432256125272240494377244848362581134804708275330834102702640621997260800 5.0122800994790495
subtotal 22 646408465298788285109296696960384431008965385800442175395725369756101721837278280054364254876808161091560552004135855723107169100329207169971214798657945600 5.2397466965279325
subtotal 23 117678776729972609905230258251927262214389962222952641950866752257607518955094651498678003707021923757854140903856321687837599834734117755193091547987968000 5.492991032546559
subtotal 24 20371732133399513811905243996686086738031826926825889456069554110113739579947224345074557197766290142610588866731856285393873933516458293989912793317376000 5.7765719654755285
subtotal 25 3539592697379126257078005738562342944973530938989649222192447415860185556779145420599548898044060146784351190667760713775427592591139200647171682598912000 5.755388790490969
subtotal 26 581638210681991328788534306898951524203808172648721484180723719387620035759234284231763363020044490486845864451699486885921413453473918264240385294336000 6.08555736602798
subtotal 27 89990024401541116688379015661960143677701756077122151543769200530780931749297959536531937413714778723121597716020642441542515410283203821577512681472000 6.463363184420145
subtotal 28 13042375728192574454784314273892028860525368036252678071510798072298675801598112212443052096348547631880092316509497142410049719391208461555717898240000 6.899818428556499
subtotal 29 1729049160985170503498686662861922253815212740394483943860447285605398698274787987258965639576422342628745675315944486255429119805493812047032352768000 7.543091325848331
subtotal 30 211811292517897513921836330470760312491314563557132099657661728161892181116114681746348022706903370637415495077267759378300100996038457707533959168000 8.16315853810802
subtotal 31 23777238090026927635128829521615052717210589612697530159965256638714135341539620207568863336582406236570466564144198176588531138096529425179869184000 8.90815374417852

Конечно же, ни одна задача до конца не досчиталась. Самая "простая" задача (3,3,4,12) считалась бы около месяца в один поток. Однако я разбил эту задачу на 32 независимые подзадачи и таки посчитал её полностью в несколько потоков за неделю.

Итак, количество совершенных паросочетаний для задачи (3,3,4,12) равно

1625717055687087185033963109305255897277923694223134108704431016700065145660728200396800

В этом числецифр. Разложение на простые множители следующее (раскладывал тут):

$\scriptsize{1625717055687087185033963109305255897277923694223134108704431016700065145660728200396800 =}$ $=2^{83} \cdot 3^6 \cdot 5^2 \cdot 7 \cdot 11 \cdot 467 \cdot 32470531 \cdot {}$ ${} \cdot 7899 339661 786140 369486 826344 309653 421924 810349$

В последнем простом сомножителецифр.

Прогнозируем ответы для посчитанных задач

Это он ещё экспонентой не экстраполировал

Частично посчитанные суммы позволяют примерно прикинуть порядок итогового ответа (как это расписано в первой части). Я просто загрузил результаты вычислений в Excel:

Во втором столбце каждой из таблиц показаны результаты частичного суммирования первых 2^i слагаемых (это которые скобки в степени), гдев идеале идёт отдо nmk . Однако до nmk у нас ни одно вычисление не дошло, поэтому мы его прогнозируем. В третьем столбце записано отношение частичной суммы для i-1 к частичной сумме для. Другими словами — во сколько раз уменьшился ответ в текущей строке по отношению к предыдущей.

В самом низу каждой из таблиц прогнозное значение ответа для всей задачи. Оно рассчитывается следующим образом: мы считаем среднее геометрическое по третьему столбцу (оно записывается в правый нижний угол таблицы), чтобы определить во сколько раз в среднем уменьшается ответ при переходе от одной строки к следующей. Затем мы делим значение в первой строке (для i=0 ) на это среднее в степени nmk . Формально:

$S_{nmk} = S_0 / \left( \sqrt[\Large{z}]{x_1 x_2 \cdots x_z} \right)^{nmk}.$

Например, для задачи (3,3,4,12) оценка получилась около $8.3 \times 10^{90}$ , хотя верный ответ равен примерно $1.6 \times 10^{87}$ . Действительно, данные оценки несколько завышены, поскольку от строчки к строчке ответ уменьшается с ускорением, как видно из следующего графика:

Ось по вертикали в логарифмическом масштабе

В "скоростях" уменьшения ответов заметны изломы в точках, кратных, а еще в задаче (4,4,4,4) видны дополнительные изломы в точках, кратных. В остальном кривые растут сверхэкспоненциально.

Для задачия бы оценил среднюю скорость (на промежутке от 1 до 64) какили, а неиз таблицы. Это даст вместо оценки в таблице порядка $10^{73}$ что-то между $10^{58}$ и $10^{63}$ , что недалеко от оценки @vvvphoenix в $4 \times 10^{56}$ в его статье.

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

Заключение

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

В следующей части мы дополнительно ускорим программу, используя следующее наблюдение. Заметим, что если в какой-нибудь шляпе поменять местами два слоя, то соответствующее ей слагаемое в формуле не изменится. Значит, можно разбить все шляпы на так называемые классы эквивалентности, и для каждого класса посчитать соответствующее слагаемое только один раз. Мы погрузимся в теорию групп, раскраски, и вот в это вот всё. Но подробности уже в следующей части.

Hubs: