Permutation question from linear optics

by_cycle_structure

As suggested by Timothy Chow:

kappa = 1

n
tau	cycles						freq

2
1	(1, 1)						1

4
1	(1, 1, 1, 1)				1
2	(2, 1, 1)					2

6
1	(1, 1, 1, 1, 1, 1)			1
2	(2, 1, 1, 1, 1)				6
3	(2, 2, 1, 1)				3
3	(3, 1, 1, 1)				2

8
1	(1, 1, 1, 1, 1, 1, 1, 1)	1
2	(2, 1, 1, 1, 1, 1, 1)		12
3	(2, 2, 1, 1, 1, 1)			20
4	(2, 2, 2, 1, 1)				4
3	(3, 1, 1, 1, 1, 1)			8
4	(3, 2, 1, 1, 1)				8
4	(4, 1, 1, 1, 1)				2


-----------------------

kappa = 0

n
tau	cycles					freq

2
1	(2)						1

4
1	(2, 1, 1)				4
2	(3, 1)					4

6
1	(2, 1, 1, 1, 1)			9
2	(2, 2, 1, 1)			15
3	(2, 2, 2)				1
2	(3, 1, 1, 1)			18
3	(3, 2, 1)				6
3	(4, 1, 1)				9

8
1	(2, 1, 1, 1, 1, 1, 1)	16
2	(2, 2, 1, 1, 1, 1)		80
3	(2, 2, 2, 1, 1)			40
2	(3, 1, 1, 1, 1, 1)		48
3	(3, 2, 1, 1, 1)			112
4	(3, 2, 2, 1)			16
4	(3, 3, 1, 1)			8
3	(4, 1, 1, 1, 1)			48
4	(4, 2, 1, 1)			24
4	(5, 1, 1, 1)			16

MO431422.py

#!/usr/bin/pypy3

from collections import Counter
from itertools import permutations

# For fast calculation of connected components in a graph
class UnionFind:
	def __init__(self, n):
		self.parent = list(range(n))
		self.rank = [0] * n


	def find(self, x):
		if self.parent[x] == x:
			return x

		y = self.find(self.parent[x])
		self.parent[x] = y
		return y


	def merge(self, x, y):
		rx = self.find(x)
		ry = self.find(y)
		if rx == ry:
			return

		if self.rank[rx] > self.rank[ry]:
			self.parent[ry] = rx
		elif self.rank[rx] < self.rank[ry]:
			self.parent[rx] = ry
		else:
			self.parent[ry] = rx
			self.rank[rx] += 1


# Probably far more than necessary
Catalan = [1, 1, 2, 5, 14, 42, 132, 429, 1430, 4862, 16796, 58786, 208012, 742900, 2674440, 9694845, 35357670, 129644790, 477638700, 1767263190]


def display(tau):
	return str(tuple(x + 1 for x in tau))


def cycle_structure(tau):
	rv = []
	seen = 0
	for x0 in range(len(tau)):
		if seen & (1 << x0):
			continue

		# New cycle
		x = x0
		cycle_len = 0
		while (seen & (1 << x)) == 0:
			seen |= 1 << x
			cycle_len += 1
			x = tau[x]

		rv.append(cycle_len)

	return tuple(rv)


def xi(tau):
	n = len(tau)
	uf = UnionFind(n)
	rot = tau[1:] + tau[:1]

	for i in range(0, n, 2):
		uf.merge(i, i+1)
		uf.merge(rot[i], rot[i+1])

	return len(set(uf.find(x) for x in range(n)))


def kappa(tau):
	return xi(tau) + len(cycle_structure(tau)) - len(tau)


def all_as(n):
	aa = Counter()
	for tau in permutations(range(n)):
		cycles = cycle_structure(tau)
		hash = len(cycles)
		kappa = xi(tau) - (n - hash)
		term = 1
		for ci in cycles:
			term *= -Catalan[ci - 1]

		aa[kappa] += term

	return aa


def kappa_distrib(n):
	return Counter(kappa(tau) for tau in permutations(range(n)))


if __name__ == "__main__":
	print("Distribution of kappa frequencies")
	for n in range(2, 12, 2):
		print(list(reversed(sorted(kappa_distrib(n).items()))))

	print()
	print()

	print("Values of a_kappa")
	for n in range(2, 14, 2):
		print(list(reversed(sorted(all_as(n).items()))))

results.txt

Frequency distribution of kappa
ell         kappa=1         kappa=0        kappa=-1        kappa=-2        kappa=-3        kappa=-4        kappa=-5        kappa=-6        kappa=-7        kappa=-8
1                 1               1
2                 3               8               9               4
3                12              58             156             238             192              64
4                55             408            1878            5416           10175           12032            8052            2304
5               273            2831           19290           86430          274209          614463          951060          963220          569568          147456
           A001764?        A062236?


Values of a_kappa

ell             kappa=1             kappa=0            kappa=-1            kappa=-2            kappa=-3            kappa=-4            kappa=-5            kappa=-6            kappa=-7            kappa=-8            kappa=-9           kappa=-10           kappa=-11           kappa=-12
  1                   1                  -1
  2                  -1                   4                   1                 -20
  3                   2                 -16                 -10                 376                -672               -2688
  4                  -5                  64                  70               -4096               16107               68704             -453324             -988416
  5                  14                -256                -420               34560             -221298             -987712            16515720            -5528960          -438848192          -716931072
  6                 -42                1024                2310             -250880             2300298            10592512          -329369950          1075857344         17272388792        -72987096128       -634446096960       -866834841600
  7                 132               -4096              -12012             1648640           -20120100           -94531584          4779253180        -31427929472       -356795260976       4688113960320      16300296812800    -248341830442496   -1308713750753280   -1577448898560000
                Catalan         signed 4^k?     signed A002802?

Under pypy3 the calculation for ell=6 took just under 8 minutes and for ell=7 took about 27.4 hours.