MO419698

MO419698.py

"""
Identify and count (by size) the subsets with odd beta_n using MacMahon's inclusion-exclusion formula for beta_n
as presented in theorem 2.3 of https://arxiv.org/pdf/1710.11033.pdf
"""

def bit_subsets_of(bitmask):
	# Nice bit hack from Thomas Beuman
	current = bitmask
	yield current
	while current:
		current = bitmask & (current - 1)
		yield current


def bitwise_decomposition(mask):
	for part1 in bit_subsets_of(mask):
		if part1 == mask:
			yield (part1,)
		elif part1 != 0:
			for tail in bitwise_decomposition(mask ^ part1):
				yield (part1,) + tail


def contributing_subsets(n):
	for perm in bitwise_decomposition(n):
		# What we're permuting are the differences, so we want to return the partial sums
		undelta = [0]
		for d in perm:
			undelta.append(undelta[-1] + d)
		yield tuple(undelta[1:-1])


def fq(n):
	# Indexing by subsets of 1 to n-1, so width = 1 << (n-1).
	# But using bytearray we get 8 bits per entry, so
	all_bits = (1 << (n-1)) - 1
	bitset = bytearray(((1 << (n-1)) + 7) >> 3)

	for ss in contributing_subsets(n):
		ss_mask = sum(1 << (i-1) for i in ss)
		# Every superset of ss gets toggled
		for addendum in bit_subsets_of(all_bits ^ ss_mask):
			superset = ss_mask | addendum
			bitset[superset >> 3] ^= 1 << (superset & 7)


	# To get good cache locality it's necessary to process ss in order
	hamming = [0] * 256
	for i in range(1, 256):
		hamming[i] = (i & 1) + hamming[i >> 1]

	rv = [0] * n
	for ss in range(all_bits + 1):
		if (bitset[ss >> 3] >> (ss & 7)) & 1:
			# Compute the Hamming weight of ss
			tmp = ss
			weight = 0
			while tmp:
				weight += hamming[tmp & 255]
				tmp >>= 8

			rv[weight] += 1

	return rv


for n in range(1, 32):
	print(n, fq(n))

results.txt

1 [1]
2 [1, 1]
3 [1, 0, 1]
4 [1, 3, 3, 1]
5 [1, 2, 2, 2, 1]
6 [1, 3, 4, 4, 3, 1]
7 [1, 0, 9, 12, 9, 0, 1]
8 [1, 7, 21, 35, 35, 21, 7, 1]
9 [1, 6, 16, 26, 30, 26, 16, 6, 1]
10 [1, 7, 22, 42, 56, 56, 42, 22, 7, 1]
11 [1, 4, 15, 52, 112, 144, 112, 52, 15, 4, 1]
12 [1, 9, 37, 93, 162, 210, 210, 162, 93, 37, 9, 1]
13 [1, 6, 24, 86, 231, 420, 512, 420, 231, 86, 24, 6, 1]
14 [1, 7, 30, 110, 317, 651, 932, 932, 651, 317, 110, 30, 7, 1]
15 [1, 0, 55, 192, 541, 976, 1291, 1312, 1291, 976, 541, 192, 55, 0, 1]
16 [1, 15, 105, 455, 1365, 3003, 5005, 6435, 6435, 5005, 3003, 1365, 455, 105, 15, 1]
17 [1, 14, 92, 378, 1092, 2366, 4004, 5434, 6006, 5434, 4004, 2366, 1092, 378, 92, 14, 1]
18 [1, 15, 106, 470, 1470, 3458, 6370, 9438, 11440, 11440, 9438, 6370, 3458, 1470, 470, 106, 15, 1]
19 [1, 12, 75, 340, 1242, 3672, 8614, 15852, 22836, 25784, 22836, 15852, 8614, 3672, 1242, 340, 75, 12, 1]
20 [1, 17, 137, 697, 2516, 6868, 14756, 25636, 36686, 43758, 43758, 36686, 25636, 14756, 6868, 2516, 697, 137, 17, 1]
21 [1, 14, 100, 502, 1997, 6496, 17200, 36752, 63154, 87308, 97240, 87308, 63154, 36752, 17200, 6496, 1997, 502, 100, 14, 1]
22 [1, 15, 114, 602, 2499, 8493, 23696, 53952, 99906, 150462, 184548, 184548, 150462, 99906, 53952, 23696, 8493, 2499, 602, 114, 15, 1]
23 [1, 8, 83, 688, 3687, 13816, 38877, 85792, 152754, 224768, 279734, 300128, 279734, 224768, 152754, 85792, 38877, 13816, 3687, 688, 83, 8, 1]
24 [1, 21, 211, 1351, 6195, 21679, 60249, 136629, 257754, 410210, 556206, 646646, 646646, 556206, 410210, 257754, 136629, 60249, 21679, 6195, 1351, 211, 21, 1]
25 [1, 18, 162, 990, 4662, 17910, 57274, 153018, 341343, 635732, 989604, 1289484, 1408212, 1289484, 989604, 635732, 341343, 153018, 57274, 17910, 4662, 990, 162, 18, 1]
26 [1, 19, 180, 1152, 5652, 22572, 75184, 210292, 494361, 977075, 1625336, 2279088, 2697696, 2697696, 2279088, 1625336, 977075, 494361, 210292, 75184, 22572, 5652, 1152, 180, 19, 1]
27 [1, 12, 121, 1072, 6970, 33032, 119098, 339632, 788135, 1519860, 2477375, 3464480, 4210476, 4488176, 4210476, 3464480, 2477375, 1519860, 788135, 339632, 119098, 33032, 6970, 1072, 121, 12, 1]
28 [1, 21, 219, 1531, 8136, 35028, 125980, 383232, 990129, 2176089, 4073847, 6506835, 8881208, 10372176, 10372176, 8881208, 6506835, 4073847, 2176089, 990129, 383232, 125980, 35028, 8136, 1531, 219, 21, 1]
29 [1, 14, 146, 1326, 9235, 48044, 192132, 610860, 1586497, 3435762, 6305230, 9939090, 13616811, 16373608, 17397304, 16373608, 13616811, 9939090, 6305230, 3435762, 1586497, 610860, 192132, 48044, 9235, 1326, 146, 14, 1]
30 [1, 15, 160, 1472, 10561, 57279, 240176, 802992, 2197357, 5022259, 9740992, 16244320, 23555901, 29990419, 33770912, 33770912, 29990419, 23555901, 16244320, 9740992, 5022259, 2197357, 802992, 240176, 57279, 10561, 1472, 160, 15, 1]
31 [1, 0, 285, 1980, 14215, 71760, 291715, 971640, 2742385, 6685520, 14204349, 26215940, 42026655, 58640480, 71505675, 76363152, 71505675, 58640480, 42026655, 26215940, 14204349, 6685520, 2742385, 971640, 291715, 71760, 14215, 1980, 285, 0, 1]