User:DroneBetter/dimensional INT enumerator

Notice

This page is not finished, much remains to be better explained, and the implementation of apg's algorithm still takes hyperexponential time with respect to dimensions (and, to be optimised, will need to generate actions as conjugacy classes instead). Many other decisions were made arbitrarily without elegance in mind and forgotten about, the article will be modified as the program itself is reduced. apg also sent me an email with some more information during the article's composition.

I modified the program to move some more common patterns I identified into functions, while reading over parts I had otherwise dismissed as complete in the process of making the article, the first steps need to be updated.

Edit 2023-10-04: The section in my page, miscellaneous curiosities#hypercube colouring sequences, contains a version that deduplicates conjugacy classes, and also some other results regarding sequences enumerating functions under group products.

Article

So you've learned about the Pólya enumeration theorem and decided to make a program for finding numbers of isotropic non-totalistic transitions in higher dimensions? Excellent!

Idea

For standard INT rules, each transition corresponds with a nonequivalent 3*3 binary matrix (of which there are 102). In the table in the INT article, the convention is to consider the neighbourhood excluding the central cell (that is being determined), so we will divide our final result by 2.

We will implement three algorithms for the generalisation of this task (of counting transitions) to higher neighbourhood ranges and dimensions, computing it in different ways (increasing in complicatedness and decreasing in time complexity).

The first algorithm will, for a given side length ('size') and number of dimensions ('dims'), iterate over all 2^sizedims binary hypercubes. For each one, it will enact each group action upon it, and determine the lexicographic minimum of the states (lists or integers) that they each return. If the identity action returns this value, a counter is incremented, otherwise it's considered a duplicate. Upon finishing, the counter is the number of nonequivalent states.

This takes exponential time with respect to the neighbourhood size to compute so can't go very far, and would be an uninteresting use of our time (certainly not worth writing an article about), but the enumeration theorem allows us to compute it more efficiently.

As seen in the example in the enumeration theorem article, we can compute closed-form equations manually for two dimensions, using our human visual intuition and knowledge of computing things like areas of triangles comprised of square cells (by summations of linear sequences). However, in three dimensions, it becomes more complicated (and error-prone) to think about, and in any higher ones, intuition mainly stops being useful, and it eventually becomes an entirely analytic task to do by considering cycles under each action separately.

Our second algorithm will have even less intuition of higher-dimensional geometry than us, but will classify and actions by equivalence, and use the fact that for each action classification, the number of instances of each cycle (corresponding with the terms within the exponent for its term) are polynomials, of the same degree as the number of dimensions of the shape. For the odd and even parities, for each cycle type, it can compute one more term in the sequence of occurrences than the degree, and recursively find the coefficient of the greatest exponent in its polynomial expansion (using the iterated difference) and subtract it, to get the corresponding polynomial for each parity, and sum these inside powers of 2, multiply by classification occurrence count, sum and divide by the total number of actions.

Both of these can be modified to consider only a subset of the cells of the hypercube, for instance to compute the equations for nonequivalent lattice diamonds and octahedra (and more generally orthoplexes) instead, but have hyperexponential slowdown for higher dimensions.

Our third algorithm will be a very efficient one specifically for hypercubes, as explained to DroneBetter by Adam P. Goucher.

Similarly to the eightfold reducer, we will have a bitwise mode (controlled by the bitwise variable, determining whether states are lists or integers). However, while for the reducer this was to maintain compatibility with multiple programs' different conventions, for this it is because we will later want to have lists of elements in each cell that are not necessarily booleans. (The lists won't be nested by dimension, the position in the kth dimension of a cell at index i is ⌊isize^k⌋%size.) Note that there is no need for bitwise inflation here (like the reducer had for compatibility with our bitwise SWAR Life program), because the representations are only to be used internally.

Prerequisites

It is in Python, you must understand the functions from the bitwise SWAR Life article, and we will use some new ones from itertools library here. (They don't do anything that we couldn't already, and don't necessarily speed it up, but allow ideas to be expressed more concisely.)

Symbol	Name	Explanation	Example
	starmap	Equivalent to map, but unpacking the function's parameters before inputting them.	starmap(f,l)=map(λ i: f(*i),l)
	product	The Cartesian product, iterates over all combinations of an element from its first parameter, one from its second, etc.	product(a,b)=((x,y) for x in a for y in b)
	combinations	Returns tuples of all length-k subsets of an input iterable, ordered according to the original input.	combinations((0,1,2,3),3)=((0,1,2),(0,1,3),(0,2,3),(1,2,3))
	chain	Like zip, but concatenates elements together.	chain(((0,1,2),(0,1,3),(0,2,3),(1,2,3)))=(0,1,2,0,1,3,0,2,3,1,2,3)
	pairwise	zips a length-n list's first n-1 elements to its last n-1	pairwise(l)=zip(l[:-1],l[1:])
	zip_longest	zip, but stopping when all input iterables are exhausted instead of any, with a fillvalue parameter for those that end early	zip_longest((0,1,2),(3,4),fillvalue=0)=((0,3),(1,4),(2,0))

You should also know a small amount about higher-dimensional geometry terms, like that a k-facet is a generaliisation of a face to k dimensions.

First steps

Now we will define some other "elementary" functions here (ie. that are applicable beyond our specific purposes).

ORsum=(λ l: reduce(int.__or__,l,0))
dot=(λ a,b: sum(map(__mul__,a,b)))
moddiv=(λ a,b: divmod(a,b)[::-1])
fact=(λ n: reduce(int.__mul__,range(1,n+1)) if n else 1)
choose=(λ n,*k: fact(n)//(reduce(int.__mul__,map(fact,k))*fact(n-sum(k))) if all(map(λ k: 0<=k,k+(n-sum(k),))) else 0)

ORsum is the same as was explained previously, dot is the dot product, moddiv is divmod but little-endian (returning (a%b,a//b)). fact is factorial, choose is the multinomial function , the factorial of the first input divided by (the product of the factorials of all thereafter and the first minus their sum), ie. (64,2) = 6!4!*2!*(6-4-2)!.

redumulate=(λ f,l,i=None: accumulate(l,f,initial=i))
bind=(λ *f: λ *a: reduce(λ a,f: (λ a,i,f: (f(a) if i else f(*a)))(a,*f),enumerate(f),a))
funcxp=(λ f,l: λ i: reduce(λ x,i: f(x),range(l),i))
expumulate=(λ f,l: λ i: accumulate(range(l),λ x,i: f(x),initial=i))

redumulate is accumulate with reduce parameter ordering, bind composes the inputted functions (ie. bind(f,g)(x)=g(f(x))), funcxp composes one f with itself l times (functional exponentiation), expumulate is equivalent but returns intermediate values like accumulate does.

decompose=(λ n,l=None: (n>>i&1 for i in range(n.bit_length() if l==None else l)))
recompose=(λ i: reduce(int.__or__,(k<<j for j,k in enumerate(i))))
mapdec=(λ dims: expumulate(λ l: (λ i: (0,)*i+(1,)+l[i+1:])(l.index(0)),2**dims-1)((0,)*dims))

decompose converts an integer to a list of booleans corresponding with its (little-endian) binary representation, recompose reverses this, mapdec(dims) is equivalent to mapping decompose over range(2^dims) (but increments a binary counter instead of repeatedly converting, so is slightly more efficient).

lexbin

See also miscellaneous curiosities' explanation of this (this explanation of my initial process remains only for posterity)

We will be encoding hypercube facets as pairs of endpoints later on, and will determine their dimension by the number of bits changed in between. To do this, it will be useful to have functions that return integers of a specific length, the dimension of the hypercube, and containing specific numbers of 1s in their binary representations (weight), the dimension of the facet, so that we only need to bitwise XOR an integer by all of them to generate all facets with the specified dimension containing it.

lexinc=(λ n: (n^n+(u:=n&-n))//u>>2|n+u)
perms=(λ n,m: expumulate(lexinc,choose(n,m)-1)((1<<m)-1))
bindex=(λ i: reduce(λ m,i: (λ s,m,c,i,b: ((s+c,m,c*i//(i-m+1)) if b else (s,m+1,c*i//m)) if m else (s,1-b,c))(*m,*i),enumerate(decompose(i),1),(0,0,1))[0])
binget=(λ n,m,i: recompose(reduce(λ m,n: (λ o,m,i,n: (λ f: ((0,)+o,m,i) if i<f else ((1,)+o,m-1,i-f))(choose(n,m)))(*m,n),range(n-1,-1,-1),((),m,i))[0]))

lexinc(n) increments an integer to the next lexicographic permutation of its 1 bits, and is equivalent to A057168(n). perms returns an iterable of all permutations of m 1s within length n, beginning with (1<<m)-1 (m 1-bits) and recursively applying lexinc, bindex(n) returns the index of an integer within the sorted infinite sequence of all those with the same weight (ie. the number of integers lower than it with the same weight), and is equivalent to A068076(n). binget(n,m,i) returns the ith element of perms(n,m).

For debugging purposes later, we would like to print states (boolean lists) to the terminal. For 3D ones, we will print them with slices arranged horizontally, so that we can print successive states on new rows.

def printState(state):
    print(''.join('\n'+'\n'.join(''.join('o' if s else ' ' for s in map(((λ r: s>>r&1) if bitwise else state.__getitem__),chain(*(range(((i*size+k)*size+j)*size,((i*size+k)*size+j+1)*size) for k in range(size))))) for j in range(size)) for i in range(size**(dims-3))))

This is similar to the previous programs' print functions, except with the order of precedence rearranged (so that the representation of an index in base width has the last digit mapped to the x position within its square, the second-last to the y and the third-last to the square it's in).

Squares are said to have eightfold symmetry, because there are eight group actions (including the identity) under which the square is invariant, and all obliquely-displaced points can be mapped to eight. In the eightfold reducer, we explained this to be because you can independently reflect in the x and y axes and transpose them (and that applying the three booleans for these in order has 2³ possibilities, each one comprising a unique group action). The number of actions for an n-dimensional hypercube, however, is 2ⁿ (from the combinations of reflections) times n! (from the number of permutations of n axes).

actions=(λ dims: fact(dims)*2**dims)

Initialisation

Now we will have another function, not for setting the board width and height but the size (length in all axes) and number of dimensions, generating all actions as both a permutation matrix and a corresponding function upon the board representation.

convert=(λ m,i,a: (λ d,n,i: ('('*bool(i)+str(size)+'+~(')*n+a+("//"+str(size**d))*bool(d)+('%'+str(size))*(d<dims-1)+')'*n+(')'*n+'*'+str(size**i))*bool(i))(*m[i],i))
def setBoard(s,d,o=False,b=False,l=True):
    global size,dims,orthoplex,boardCells,bitwise,unravel,matrices,exps,matrixPeriods,λs
    size=s
    dims=d
    orthoplex=o
    boardCells=size**dims
    if l:
        bitwise=b
        unravel=True
        matrindex=(λ p: sum((e-sum(e>d for d,m in p[:i]))*fact(len(p)+~i)*2**dims+2**i*n for i,(e,n) in enumerate(p)))
        matrices=sorted((λ matrices: [tuple((j,n>>i&1) for i,j in enumerate(m)) for m in matrices for n in range(2**dims)])([(λ t: tuple(reduce(λ e,n: e+(e>=n),t[i-1::-1],e)%dims for i,e in enumerate(t)))(tuple((a[1]%(dims-i) for i,a in enumerate(redumulate(λ m,i: moddiv(m[1],i),range(dims,1,-1),(0,m)))))) for m in range(fact(dims))]),key=matrindex)
        compose=(λ a,b: tuple(map(λ i: (λ j,k: (j,k^i[1]))(*b[i[0]]),a)))
        exps=tuple(map(λ m: lap(matrindex,expumulate(λ i: compose(m,i),lcm(4,fact(dims)))(matrices[0])),matrices))
        matrixPeriods=tap(λ e: e[1:].index(0)+1,exps)
        λs=tuple(eval("λ s: "+( ( '|'.join((λ l: map((λ e: ''.join("(s&("+'|'.join(starmap(λ j,ind,d: '1<<'+str(ind),e[1]))+'))')+(λ d: ">>"+str(d) if d>0 else '' if d==0 else "<<"+str(-d))(e[1][0][2])),l))(map(λ x: (x[0],tuple(x[1])),groupby(sorted(((λ j,ind: (j,ind,ind-j))(j,eval('+'.join(map(λ i: convert(m,i,str(j)),range(dims))))) for j in range(boardCells)),key=(λ j: j[2])),key=(λ j: j[2])))))
                                 if unravel else
                                  "ORsum("*(not unravel)+('|' if unravel else ',').join((λ j,ind: '(s&1<<'+str(ind)+')'+(">>"+str(ind-j) if ind>j else '' if ind==j else "<<"+str(j-ind)))(j,eval('+'.join(map(λ i: convert(m,i,str(j)),range(dims))))) for j in range(boardCells))+")"*(not unravel))
                               if bitwise else
                                ( '['+','.join('s['+str(eval('+'.join(map(λ i: convert(m,i,str(j)),range(dims)))))+']' for j in range(boardCells))+']'
                                 if unravel else
                                  '[s'+'['+'+'.join(map(λ i: convert(m,i,'i'),range(dims)))+"] for i in range("+str(boardCells)+")]"))) for m in matrices)

The matrices list is equivalent to signed permutation matrices , each element represents an axis and is a tuple, containing the index of the axis that the position in it is mapped to and the sign (a boolean of whether it's flipped). convert intakes a matrix, dimension index and variable name (as a string) and returns a string with the coefficient of the variable (bitwise-NOT'ed and added to the board width if the sign is true). exps is a list of indexes of exponentiations of the matrices within the list, for our convenience, such that funcxp(λs[i],a)=λs[exps[i][a]] (because all compositions of actions (including one with itself) can be described by other actions in the group), and λs is a list of λ functions for enacting these, that are evaluated from constructed strings. In most programming, you may be inclined to keep functions general and have them iterate over a number of cells that depends on global variables or their parameters, to avoid having to manually change them for different instances of the same problem (and avoid restating lines over and over again), however we can unravel them to arbitrarily long expressions, because they are never seen by humans and this way, they don't have to invoke internal iterators (and can be redefined later when it changes, which is justified by them being called exponentially many times with respect to the number of times they're redefined). The unravel mode in the bitwise mode makes them group (OR) together bits that are shifted by the same amount, and AND the integer by this before a single shift in the ORsum, reducing the length slightly. (A more efficient way would involve something like the logarithmic-time reversal in the bitwise simulator's manifold mode, but would be more difficult to avoid overlapping in.)

Currently, on the test case of 3*3*3 binary cubes in the trivial algorithm, the non-bitwise unravelled implementation with unravel is unexpectedly the fastest, checking and counting 41000/second, whereas the bitwise unravelled one counts only 17000.

Now, for our future purposes of cycle enumeration, we would like a function that assigns each cell a unique value, or 0 if the key function called upon its index is false. If there is no key, it provides one for the orthoplex mode, that counts the sum of the absolute displacements from the centre in each axis and only accepts those for which it's less than the threshold.

permutationState=(λ key=None: lap(λ i: (i+1)*(key(i) if key else (sum(abs(x[0]-(size-1)/2) for x in tuple(expumulate(λ a: moddiv(a[1],size),dims)((0,i)))[1:])<=(size+1-size%2)/2) if orthoplex else True),range(boardCells)))

Algorithms

Trivial schoolbook one

This one (the checking and counting of states' lexicographic minimalities) is impractical for most purposes but will be used to test our others against.

We would like a function that acts like mapdec (counting in binary), but skips over cells of indexes that are False in a test state (for instance, the permutationState from previously).

def shapedec(test=None):
    nextJ=(λ j: next(filter(test.__getitem__,range(j,boardCells))))
    if test==None:
        test=permutationState()
    state=[False]*boardCells
    itsOver=False
    for i in range(2**boardCells-1):
        #printState(state)
        yield(state)
        j=nextJ(0)
        state[j]^=True
        while j<len(state):
            if test[j]:
                if state[j] or not j<boardCells-1:
                    break
                try:
                    j=nextJ(j+1)
                except:
                    itsOver=True
                    break
                state[j]^=True
            else:
                j+=1
        if itsOver or j==len(state):
            break
    if not itsOver:
        yield(state)
generateStates=(λ: range(2**boardCells) if bitwise else shapedec() if orthoplex else mapdec(boardCells))

Now, to determine the number of nonequivalent states.

def nonequivalents(size,dims,orthoplex=False):
    setBoard(size,dims,orthoplex)
    distincts=0
    interval=2**max(16,boardCells-10)
    t=time()
    for i,state in enumerate(generateStates()):
        if i and not i%interval:
            print(distincts*2**boardCells/i,distincts,i/2**boardCells,interval/(time()-t))
            print('\n'.join(map(strate,(state,min(map((λ f: f(state)),λs))))))
            t=time()
        distincts+=(state==min(map((λ f: f(state)),λs)))

This could all be done in a closed form as well, but it returns up to using intermediate results inline to show that it's still running (and can take a long time for some tasks). Anyway, now you can determine the states of an n-hypercube of edge length k with

nonequivalents(s,dims,orthoplex)

(in O(n!*2^n+kn) time).

Polynomial fitter

We will make a function that, for each action, will use its period and find the number of cells that oscillate at subperiods (factors of this) by seeing whether enacting the exponent of the matrix causes their index to contain the same value, account for this method counting cells within subperiods that oscillate at their factors (by subtracting all those preceding it that are factors).

def cycleLengths(size,dims,orthoplex=False):
    setBoard(size,dims,orthoplex)
    test=permutationState()
    factors=tuple(tuple(filter(λ n: not m%n,range(1,m//2+1)))+(m,) for m in matrixPeriods)
    subperiods=tuple(tuple(sum(map(((λ a,b: a==b!=0) if orthoplex else int.__eq__),test,λs[exps[i][a]](test))) for a in f) for i,f in enumerate(factors))
    return((λ t: lap(λ t: reduce(λ r,t: r+((t[0],t[1]-sum(i[1] for i in r if not t[0]%i[0])),),t,()),t))(tuple(map(λ a,b: tuple(zip(a,b)),factors,subperiods))))

The factorisation will take quadratic time, but we don't mind, it is a negligible proportion of our program's total cost.

To obtain the number of cycles within an action, we must divide each cycle period's occurrences by its length (which will always return an integer, so we can use floordivs). From the cycleLengths function's output, we can determine the number of states by applying the Pólya enumeration theorem directly, for specific cases.

specific=(λ cycles: sum(2**sum(i[1]//i[0] for i in c) for c in cycles)//len(cycles))

Such that specific(cycleLengths(size,dims,orthoplex)) will be equivalent to nonequivalents(size,dims,orthoplex), but a great deal more efficient (still having the factorial times exponential number of actions and maximum period to be checking, but multiplied by the period and the number of cells instead of 2^{number of cells}).

Now, to find minimal fitting polynomials from empirical data, we must consider parities separately, as you saw in the explanation of the theorem (where even widths had no orthogonally-displaced cells to be considering). Parity mod 2 suffices for all dimensions of hypercube, but orthoplexes are more complicated. For the 2D (diamond) case, the same argument applies, but the diagonal edges' centres alternate with respect to width, between being protruding and indented cells, but only changing on even-to-odd increments, because even ones' quadrants are the same as preceding odd ones' but with the boundary cells copied into them. Meanwhile, if you create some extended von Neumann neighbourhoods (odd-width regular octahedra) out of cubes (in an image-editing program or something like Minecraft) and colour the cells in a triagonal from the centre, you will notice that they oscillate with period 3 (between protruding cells and two different orientations of vertices between cells). For the odd and even parities, we need to double this. More generally, an n-dimensional one's n-agonal will have n distinct phases, so we must take the lcm of all numbers up to it. Thus, we must consider 2*A003418(dims), to account for our two cases.

parities=2*lcm(*range(1,1+dims*orthoplex))

Now we will define a function, generateEquation, that every line shown hereafter in this section is to be included in. It will return a string with a λ function's definition, containing an equation for each parity, that determines the number of cells.

Note that while specific() divides each cycle by its length (i[0]) immediately, we don't do so in equation-finding, because though it leaves the values themselves integer, it doesn't leave the coefficients of each term of the polynomial integer also, so by dividing later and computing the entire polynomial in the numerator (which will be divisible by the denominator) we avoid incurring floating-point errors.

typeCounts=tap(λ i: lap(λ g: (len(tuple(g[1])),g[0]),groupby(sorted(cycleLengths(i,dims,orthoplex)))),range(parities,parities*(dims+2)))

We begin with the size being equal to parities, so that each parity is on its second term (in case some actions in its first have 0 occurrences of some periods and forgo them and aren't classified in the same way as the later ones). For each size i, it finds the cycle lengths of each action, sorts and groups them (with lexicographic precedence), then maps the conversion (from groupby objects containing the key result (which is the value, being that no key is provided so it uses the identity function) and iterator with the run of instances that matched it, to the length of the run and the key result) to each element of the groupby. Remember, each element of typeCounts is a tuple for its corresponding size, containing classifications of actions, each one containing its occurrence count and a tuple of periods, each one being its value and occurrence count. For instance, for two dimensions, the first element (for 2*2 boards) is

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

Meaning "there are 2 actions for which all 4 cells oscillate with period 4, 3 for which 4 oscillate with period 2, 2 for which 2 have period 1 and 2 have period 2, and 1 for which 4 have period 4" (corresponding with the 90º rotations, two orthogonal reflections and one eversion/180º rotation (that it considers equivalent), two diagonal reflections and one identity). While for even widths they're grouped into the classifications, with numbers of occurrences 2,3,2,1, for odd widths they're 2,1,4,1 (the 3 is split into 1 and 2 because reflections are not equivalent to eversions (because cells along the line of symmetry are not mapped to each other), but the 2 for orthogonal reflections are combined with those of diagonal reflections (because the same number of cells are of period 1 and 2 respectively are the same for both)).

Next, we encode this as a list of parities, each one containing a tuple of action classifications, each of which is its number of instances and a tuple of the cycle counts across each instance of it.

togethers=tuple(tuple((typeCounts[i][j][0],tuple(t[j][1] for t in typeCounts[i::parities])) for j in range(len(typeCounts[i]))) for i in range(parities))

So for instance (in the 2D case), it would convert

(((2,((1,0),(2,0),(4, 4))),(3,((1,0),(2, 4))),(2,((1,2),(2, 2))),(1,((1, 4),))),
 ((2,((1,1),(2,0),(4, 8))),(1,((1,1),(2, 8))),(4,((1,3),(2, 6))),(1,((1, 9),))),
 ((2,((1,0),(2,0),(4,16))),(3,((1,0),(2,16))),(2,((1,4),(2,12))),(1,((1,16),))),
 ((2,((1,1),(2,0),(4,24))),(1,((1,1),(2,24))),(4,((1,5),(2,20))),(1,((1,25),))),
 ((2,((1,0),(2,0),(4,36))),(3,((1,0),(2,36))),(2,((1,6),(2,30))),(1,((1,36),))),
 ((2,((1,1),(2,0),(4,48))),(1,((1,1),(2,48))),(4,((1,7),(2,42))),(1,((1,49),))))

to

(((2,(((1,0),(2,0),(4,4)),((1, 0),(2, 0),(4,16)),((1, 0),(2, 0),(4,36))))
  (3,(((1,0),(2,4)      ),((1, 0),(2,16)       ),((1, 0),(2,36)       )))
  (2,(((1,2),(2,2)      ),((1, 4),(2,12)       ),((1, 6),(2,30)       )))
  (1,(((1,4),           ),((1,16),             ),((1,36),             ))))
 ((2,(((1,1),(2,0),(4,8)),((1, 1),(2, 0),(4,24)),((1, 1),(2, 0),(4,48))))
  (1,(((1,1),(2,8)      ),((1, 1),(2,24)       ),((1, 1),(2,48)       )))
  (4,(((1,3),(2,6)      ),((1, 5),(2,20)       ),((1, 7),(2,42)       )))
  (1,(((1,9),           ),((1,25),             ),((1,49),             )))))

Next, we will convert these to polynomials.

terms=tuple(tuple((f,)+tuple((i[0][0],(λ a: tap((Fraction if orthoplex else int),reduce(λ a,b: lap(__add__,a,b),(tap(λ n: n*j,(λ a,b,n: (choose(n,k)*a**k*b**(n-k) for k in range(dims+1)))(Fraction(1,parities),-1-Fraction(b,parities),i)) for i,j in enumerate(a)))))(reduce(λ i,d: (λ a,i,d: (λ c: ((c,)+a,tuple(k-c*j**d for j,k in enumerate(i))))(funcxp(λ i: tap(int.__sub__,i[1:],i[:-1]),d)(i)[0]//fact(d)))(*i,d),range(dims,-1,-1),((),lap(λ n: n[1],i)))[0])) for i in transpose(s)) for f,s in t) for b,t in enumerate(togethers))

This looks complicated but I will explain each operation enacted, in order. Within each parity, we transpose the second element of each action (with the counts of cycle lengths over each width), ie. from (((1,0),(2,0),(4,4)),((1, 0),(2, 0),(4,16)),((1, 0),(2, 0),(4,36)))) to (((1,0),(1,0),(1,0)),((2,0),(2,0),(2,0)),((4,4),(4,16),(4,36))), then for each element of this transposed list, we include the first element (the first size)'s first element (the period), and the fitting polynomial's coefficients. To do this, first we enact a reduce statement.

reduce(λ i,d: (λ a,i,d: (λ c: ((c,)+a,tuple(k-c*j**d for j,k in enumerate(i[:-1]))))(funcxp(λ i: tap(int.__sub__,i[1:],i[:-1]),d)(i)[0]//fact(d)))(*i,d),range(dims,-1,-1),((),lap(λ n: n[1],i)))[0]

It begins with the tuple ((),lap(λ n: n[1],i)), that is mutated in-place with each iteration of the reduce, the first element will be the polynomial (with the element at index i corresponding with the coefficient of the nⁱ term), the second is the second element (the number of occurrences) of the cycle length (because the first each contain the same value, the cycle length). It recursively finds the dth iterated difference (for d beginning as dims and then decrementing to 1), in each iteration it unpacks the mutating tuple and passes it with the difference order into

λ a,i,d: (λ c: ((c,)+a,tuple(k-c*j**d for j,k in enumerate(i[:-1]))))(funcxp(λ i: tap(int.__sub__,i[1:],i[:-1]),d)(i)[0]//fact(d))

, which takes the first element of the difference list (which, for each recursively-applied iteration of the difference function, decrements the length by 1, and will return a list of length 1), divides by the factorial of the difference order to get the polynomial coefficient (because the dth difference between two consecutive dth powers is d!), then prepends this to the polynomial coefficients tuple (which is generated in reverse order) and computed and subtracted from each element of the sequence (so it doesn't affect the computation of the previous exponent's coefficient). After the reduce has finished, the first element (the accumulating list) is taken, the second is discarded (because all differences down to the 0th have been found and it is empty).

However, this only returns a polynomial that returns the number of occurrences of the cycle type in the action with respect to the index within this parity's list (that extrapolates it), we would like one that returns it with respect to the board size itself, so we substitute the equation from one in terms of k to n, where n=parities*(k+1)+b (the index of this parity), so k = nparities - 1 - bparities. For each polynomial term with the power p, we convert it with the binomial expansion of ((nparities) + (-1 - bparities))^p (in terms of n), then we sum these polynomials together, achieved with

λ a: tap((Fraction if orthoplex else int),reduce(λ a,b: lap(__add__,a,b),(tap(λ n: n*j,(λ a,b,n: (choose(n,k)*a**k*b**(n-k) for k in range(dims+1)))(Fraction(1,parities),-1-Fraction(b,parities),i)) for i,j in enumerate(a))))

. (We convert it to Python's fraction object to avoid precision loss, then back again if in the hypercube case (for which all terms are integers).)

As another demonstration, this entire statement converts our togethers list (as was given in the 2D square example above) to

(((2,(1,( 0, 0, 0)),(2,(0, 0, 0)), (4,( 0, 0, 1)))
  (3,(1,( 0, 0, 0)),(2,(0, 0, 1)))
  (2,(1,( 0, 1, 0)),(2,(0,-1, 1)))
  (1,(1,( 0, 0, 1))))
 ((2,(1,( 1, 0, 0)), (2, ( 0, 0, 0)), (4, (-1, 0, 1)))
  (1,(1,( 1, 0, 0)), (2, (-1, 0, 1)))
  (4,(1,( 0, 1, 0)), (2, ( 0,-1, 1)))
  (1,(1,( 0, 0, 1)))))

, each element of the outermost list represents a parity, each within this an action classification, comprised of its occurrences followed by each cycle type's length and polynomial coefficients. We then only need to divide the polynomials by cycle lengths, sum them over each classification, exponentiate 2 to them, multiply by the occurrences, and sum these, which we make a closed-form equation for with.

return('λ n: ( '+''.join('+'.join(map(λ t: (str(t[0])+'*')*(t[0]!=1)+'2**('+(λ x: x if x else '0')('+'.join((λ x: x[0] if len(x)==1 else '('+'+'.join(x)+')')(tuple(str(n)*(n!=1 or not i)+('*'*(n!=1)+'n'+('**'+str(i))*(i!=1))*bool(i) for i,n in enumerate(i[1]) if n))+('/'+str(i[0]))*(i[0]!=1) for i in t[1:] if any(i[1])))+')',t))+('\n'+' '*11+'if n%'+str(parities)+('=='+str(parities+~i))*(parities+~i!=1)+' else\n '+' '*11)*bool(parities+~i) for i,t in enumerate(terms[::-1]))+')//'+str(actions(dims)))

We can print this to the terminal, or evaluate it for a function of board size that returns nonequivalent states. The process of finding it is inefficient (notably requiring the board size to be changed (and the action λs redefined) for each board width up to one proportional to the number of dimensions d (requiring a number of operations proportional to a sum of dth powers up to d, hence O(d^d+1) (by sum-integral convergence)), though the actions upon these results afterwards are very cheap), but the equation is the most efficient "algorithm" to find the numbers for the specific case.

apg's genius masterpiece

The previous equation-generating algorithm does the polynomial-fitting on the values it computes empirically, to circumvent a lack of the maths knowledge necessary to generate it more efficiently. However, when DroneBetter told apg about having made the program work for octahedra by considering 12 parities, he replied with an explanation of a more efficient method.

By the way, you could accelerate your calculations by counting the number of orbits, under the hyperoctahedral group, of each possible shape. I think that you get 2ⁿ distinct shapes of orbit, corresponding to the subsets of the set of vertices of the Wythoff simplex. For n = 3 dimensions, the eight different shapes of orbit can be seen by looking at the dual of an omnitruncated cube, the Disdyakis dodecahedron .

and representatives for each shape of orbit are:

• 000 (1 point): the centre of the polyhedron
• 001 (6 points): the 8-valent vertices
• 010 (12 points): the 4-valent vertices
• 100 (8 points): the 6-valent vertices
• 011 (24 points): the centres of the 8-4 edges
• 101 (24 points): the centres of the 8-6 edges
• 110 (24 points): the centres of the 4-6 edges
• 111 (48 points): the centres of the faces of the polyhedron

Then, for a given arrangement of cells (say, the 2197 cells in the range-6 3D Moore neighbourhood), you can count the numbers of each shape of orbit:

• 000: 1 = (60)
• 001: 6 = (61)
• 010: 6 = (61)
• 100: 6 = (61)
• 011: 15 = (62)
• 101: 15 = (62)
• 110: 15 = (62)
• 111: 20 = (63)

and check that 1 * 1+ 6 * 6 + 6 * 12 + 6 * 8 + 15 * 24 + 15 * 24 + 15 * 14 + 20 * 48 = 2197. (This generalises to the range-r Moore neighbourhood in n dimensions: for a given n-bit shape code S, the multiplicity is (rpopcount(S)).)

Then, for a given group element, you can easily compute the cycle type of the whole arrangement of cells by summing, over all 2ⁿ shapes of orbit, the product of the multiplicity of that shape with the cycle type of the element acting on that shape.

Things get more complicated if the symmetry group under consideration isn't a Coxeter group -- but unless you're going to be looking at (for example) neighbourhoods of cells in the Leech lattice, you're unlikely to run into this issue.

We will implement this, everything hereafter will be within a disdyakis function (that will intake only a dims parameter).

Very inefficient way

At first, I thought it would be a good idea to generate all of this disdyakis dodecahedron (or more generally polytope)'s vertices, then compare them pairwise (checking whether they lay in the same simplex (with each vertex the centre to the cube vertex) to determine edges, then triowise from edges for faces (checking whether each pair of edges shares at least one vertex), and so on, however this is extremely inefficient and I realised that it could be solved equivalently by only considering a single such simplex ().

Some of the bitwise functions we had until now forgone will become applicable now.

weight=(λ d: tap(λ w: tuple(expumulate(lexinc,choose(d,w)-1)((1<<w)-1)),range(d+1)))

This returns a list of lists of each of the (dw) weight-w length-d numbers in lexicographic order, for w from 0 to d (inclusive).

weave=(λ a,b: funcxp(λ x: (λ a,b,j: (λ j: (a|(b&1)<<j,b>>1,j+1))(next(filter(λ i: ~a>>i&1,range(j,j+a.bit_length()+1)))))(*x),b.bit_length())((a,b,0))[0])

This "weaves" b's bits into a's 0-bits.

agonals=tuple(tuple((x,tuple(tuple(weave(x,k) for k in j) for j in weight(dims-i))) for x in w) for i,w in enumerate(weight(dims)))

When using the convention of encoding hypercube vertices' coordinates in each axis as bits at corresponding indexes, this contains, for each vertex, a tuple of cells that are k edges away in the positive direction in each axis, for all possible k (dims minus x's weight).

facets=tap(λ f: tuple((a[0],d) for i in agonals for a in i if len(a[1])>f for d in a[1][f]),range(dims+1))

facets is a list of lists of k-facets with respect to k, constructed from this.

Now we will show the different ways of generating the representatives of the facets.

vertices=tap(λ n: tarmap(λ c,i: construce(λ x,j,d: ((x+((-1)**(i>>j&1),),j+1) if d else (x+(0,),j)),decompose(c,dims),((),0))[0],product(perms(dims,n),range(1<<n))),revange(dims+1))

For each n in the range from 0 to dims inclusive (reversed so we begin with cube vertices, then edge centres, then face centres and so on), for each combination of c (a permutation of n 1s (dimensions with displacement) in length dims) and i (an integer of length n), we construct a tuple, equivalent to decomposing our n but with its 1s' signs flipped according to j's bits (incrementing j (its index in i) only when d (the current bit of x iterated over) is 1). Meanwhile, unsigned vertices are much easier (single canonical representations of each hypercube facet (disdyakis vertice), with the convention of being sorted in reverse and all values being 1 or 0).

unsignedVertices=tuple(((1,)*n+(0,)*(dims-n),) for n in range(dims,0,-1))

Now, for generating coordinates corresponding with our disdyakis polytope's facets' centres instead of its vertices.

disdyakisFacets=tuple(redumulate(λ x,i: tarmap(λ x: tuple(tap(( (λ s: tap(sum,transpose(s[1])))
                                                               if i-2 else
                                                                (λ s: s[1])),s)
                                                          for s in product(*map(λ x: map(λ y: (x[0],y),x[1]),x))
                                                          if all( any(map(λ x: tuple.__eq__(*x),tuple(product(a,b))))
                                                                 if i-2 else
                                                                  all(map(λ a,b: a<=b,pairwise(map(λ i: i[1],sorted(enumerate(b),key=(λ i: a[i[0]]))))))
                                                                 for j,(k,a) in enumerate(s)
                                                                 for l,b in s[:j])),
                                               combinations(enumerate(x),i)),
                                 range(2,dims+1),
                                 (unsignedVertices if unsignedMode else vertices)[:-1]))

Our idea here is to form lines from pairs of vertices, then faces from trios of lines, then cells from quartets of faces, and so on. Beginning with the disdyakis vertices, we recursively enact our redumulate statement's λ function, with x being our result from the previous iteration and i being the dimension of the facet we're considering (which, for a disdyakis polytope, is a simplex). The function generates all combinations (ie. subsets, without replacement) of i types (separated by their bitwise representation in apg's encoding scheme) of these (i-1)-facets, and creates a unique point (the sum, representing the average but integer) for all combinations of vertices across types that lie in the same simplex (determined by whether, for each pair of vertices, sorting one's axial displacements according to the other's values makes it nondecreasing) or (in subsequent iterations) that all pairs of (i-1)-facets share one (i-2)-facet. It thus begins with vertices being coordinates and all higher facets being tuples of preceding facets.

This is inefficient, however, because each iteration the maximum axial displacement doubles, so we would have to store them in a hypercube with exponential size with respect to dimensions, and we iterate over many times more displacements (exponentially with respect to dimensions) than is necessary for our lack of symmetry, to the power of the number of dimensions (as it accumulates over successive iterations), so it will become infeasible very quickly. However, in making it, I realised it to be an extremely overcomplicated way. We can instead (less interestingly) look at combinations of i disdyakis vertices each iteration, and ensure they lie in the same simplex by beginning with our unsignedVertices list (that describes its boundaries, so all averages will), and looking at combinations of vertices instead of (i-1)-facets causes the edge length of the hypercube containing it to have to grow only linearly.

disdyakisFacets=tap(λ i: tarmap(λ e,x: tuple(tap(sum,transpose(map(λ s: s[1],s))) for s in product(*map(λ x: map(λ y: (x[0],y),x[1]),x))),enumerate(combinations(enumerate(unsignedVertices),i))) if i else (((0,)*dims,),),range(dims+1))

(This generates single canonical representations of each integer in the scheme, such that displacements in each axis are nonnegative, the axes are sorted in reverse by displacement and the differences between consecutive axes' displacements are at most 1.)

And we can add signs and permutations to each one by enacting each of our permutation matrices upon them.

duplicatedFacets=tap(λ f: tap(λ f: tuple(set(tap(λ f: enact(*f),product(matrices,f)))),f),disdyakisFacets)

Remember that apg's algorithm requires us to know the number of instances of each facet, ie. the six 8-valent vertices (cube faces), twelve 4-valent ones (cube edges) and eight 6-valent ones (cube vertices), and the 24 lines joining elements of pairs of this, and the 48 faces, so we would like to generate ((1,),(6,12,8),(24,24,24),(48,)).

actionCounts=tap(λ f: tuple(reversed(tuple(decompose(map(λ f: lap(λ f: choose(dims,*tap(λ x: len(tuple(x[1])),groupby(sorted(f))))*2**sum(map(bool,f)),f),f))))),disdyakisFacets)

For each facet weight, for each facet representative (of which we have a single one, with a list of axial displacements), we take

λ f: choose(dims,*tap(λ x: len(tuple(x[1])),groupby(sorted(f))))*2**sum(map(bool,f))

, which is (the total number of dimensions*((the number of axes with the displacement) for each displacement)) (using the choose function), times 2 to the power of the number of dimensions in which it is displaced by a nonzero amount.

Now we make a function to find the period of each matrix over each facet.

def period(matrix,facet):
    f=enact(matrix,facet)
    i=1
    while f!=facet:
        f=enact(matrix,f)
        i+=1
    return(i)

We combine our matrices into terms of a classification by the numbers of cells of each period in the edge-length-(2*dims+1) hypercube (that contains cells of each symmetry classification). (This is the part, I think, of which apg explained to me an optimisation that I haven't yet enacted.)

RLEmatrices=tap(λ g: (λ g: (len(g),*g[0]))(tuple(g[1])),groupby(sorted(zip(cycleLengths(3,dims,orthoplex),enumerate(matrices))),λ x: x[0]))

After doing this (again with groupby), it groups them by their tuple of periods' values and numbers of cells oscillating with them under each matrix of the equivalence class, then converts them to tuples containing the number of occurrences, the period tuple and the first matrix of the class (as its index within matrices and explicit value).

Now, we generate, for each classification in RLEmatrices, a tuple of (its occurrences, its oscillatory period and, (for each disdyakis facet's weight in apg's encoding, (for each facet of the weight, the sum of (the numbers of occurrences of cells of each period divided by the period, over all cells representing the facet)))).

cycleTypes=tap(λ m: (m[0],m[1][-1][0],tap(λ f: tap(λ f: sum(map(λ g: g[0]//tuple(g[1])[0],groupby(sorted(map(λ f: period(m[2][1],f),f))))),f),duplicatedFacets)),RLEmatrices)

Now we will make a list of, for each weight in apg's encoding, for each facet of the weight (represented as a cell of a corresponding displacement under symmetry), a state with 1 cells at all representatives (not under symmetry) and 0s elsewhere.

shapes=tap(λ d: tap(λ s: (λ s: permutationState(λ i: s[i]))(tuple(decompose(ORsum(map(λ c: 1<<sum(starmap(λ i,f: (f+dims)*size**i,enumerate(c))),s)),boardCells))),d),duplicatedFacets)

For instance, for 2D, each one would look like a square in this list

     |     |     | o o 
     | o o |  o  |o   o
  o  |     | o o |     
     | o o |  o  |o   o
     |     |     | o o 

And for 3D, each one would be a cube (represented by a horizontal list of squares) in this vertical list

                                                 
                                                 
                                                 
                        o                        
                                                 
                                                 
                                                 

                                                 
                                                 
                o o           o o                
                                                 
                o o           o o                
                                                 
                                                 

                                                 
                                                 
                 o     o o     o                 
                o o           o o                
                 o     o o     o                 
                                                 
                                                 

                                                 
                                                 
                        o                        
                 o     o o     o                 
                        o                        
                                                 
                                                 

                                                 
         o o   o   o         o   o   o o         
        o   o                       o   o        
                                                 
        o   o                       o   o        
         o o   o   o         o   o   o o         
                                                 

                                                 
                o o           o o                
         o o   o   o         o   o   o o         
                                                 
         o o   o   o         o   o   o o         
                o o           o o                
                                                 

                                                 
                 o     o o     o                 
          o           o   o           o          
         o o   o   o         o   o   o o         
          o           o   o           o          
                 o     o o     o                 
                                                 

         o o   o   o         o   o   o o         
  o o         o     o       o     o         o o  
 o   o o     o                     o     o o   o 
                                                 
 o   o o     o                     o     o o   o 
  o o         o     o       o     o         o o  
         o o   o   o         o   o   o o         

Now, we will make a function that intakes a maximum period (p), λ (l) and shape (shape), and returns the numbers of cycles of cells of subperiods.

cycles=(λ p,l,shape: (λ f: reduce(λ f,i: f[:i+1]+lap(λ e: (e[0]-f[i][0]*(not e[1]%f[i][1]),e[1]),f[i+1:]),range(len(f)),f))(lap(λ i: (sum(map(λ a,b: a==b!=0,λs[exps[l][i]](shape),shape)),i),factorise(p))))

It first finds, for each factor of the classification's period, the number of cells in the shape that are equal under the sample action to the power of the factor, and generates this zipped with the periods, as f, which is passed to the function that subtracts, from each period, the number of cells invariant under subperiods.

Now we can use this to find, for each action classification, the number of distinct points (that are not mapped to each other under its exponentiation) in each shape (which, remember, represents a disdyakis facet type) by dividing each period's occurrences by its value again.

shapeDistincts=tap(λ m: tap(λ s: (tap(λ s: sum(map(λ x: x[0]//x[1],cycles(m[1][-1][0],m[2][0],s))),s)),shapes),RLEmatrices)

Now we are very close, however apg's explanation gave us our result equation with the exponents in terms of choose functions (with the first value changing and the second remaining constant), we would like to convert them to polynomials.

choosopoly=(λ k: (reduce(λ a,k: tarmap(λ a,b: a-k*b,pairwise((0,)+a+(0,))),range(k),(1,)),fact(k)))

For a given value of k, this returns a tuple representing a fraction, the second element being the integer denominator, the first the coefficients of powers of n. Remember that in Pascal's triangle (left-justified), each column (n choose k with respect to n) is the accumulated sum of the one to its left, so it describes the side-length-n right-angled simplex that fills every cell it lies partially within, meaning it will be asymptotic to n**k/k!, so the last term of the numerator will be 1. We can find the polynomial for each k by beginning with 1 and recursively applying the summation formula, on the kth iteration the coefficient of each exponent is the preceding exponent's coefficient in the previous iteration minus k times its own (which, bearing in mind that the entire thing will be divided by k each iteration by the factorial, is equivalent to using the preceding's previous value divided by k minus its own). The table of values with respect to iterations is

0,((1,), 1)
1,((0, 1), 1)
2,((0, -1, 1), 2)
3,((0, 2, -3, 1), 6)
4,((0, -6, 11, -6, 1), 24)
5,((0, 24, -50, 35, -10, 1), 120)
6,((0, -120, 274, -225, 85, -15, 1), 720)
7,((0, 720, -1764, 1624, -735, 175, -21, 1), 5040)

, the coefficients in the numerator comprise the terms of A048994, the Stirling numbers (about which Mokhtar Mohamed commented "the triangle gives coefficients T(n,k) of x^k in the expansion of C(x,n) = a(k)*x^kn!" in 2012, so it is correct). Note that slightly interestingly, the sum of the absolute values of the coefficients is equal to the factorial in the denominator.

Now, because the initial equation apg's algorithm returns contains coefficients of choose functions on successive values of k in a run, we can sum these coefficients times polynomials together by casting the polynomials to fractions first, to sum to get a minimal form (with one term per exponent of n in each exponent of 2).

polyopoly=(λ d: tap(sum,zip_longest(*map(λ d,i: (λ c: (λ f: tap(λ c: c*f,c[0]))(Fraction(sum(d),c[1])))(choosopoly(i)),d,range(dims+1)),fillvalue=0)))

We calculate the choosopoly of the index of each coefficient and pass it into our nested function as c, then create the fraction of the coefficient (from the algorithm) divided by the denominator (from the choosopoly's factorial), and multiply each term of the polynomial by this. We use zip_longest(*l) instead of transpose (equivalent to zip(*l)), because they're of varying degree, then sum together each term (which contains a list of fractional coefficients for each exponent).

Now we can construct the final equation.

return(shapeDistincts,RLEmatrices,actionCounts,'λ n: ('+'+'.join(map(λ d,m: str(m[0])+'*2**('+'+'.join(str(c.numerator)+('*n'+('**'+str(i))*(i!=1))*bool(i)+(λ d: ('/'+str(d))*(d!=1))(c.denominator) for i,c in enumerate(polyopoly(d)) if c)+')',shapeDistincts,RLEmatrices))+')//'+str(actions(dims)))

We return shapeDistincts,RLEmatrices and actionCounts as well as this equation, so that the program outside can use them to calculate the number for specific widths itself in a more straightforward way (for debugging purposes).

Remember, by the enumeration theorem we are taking the average of (2^{number of cycles} for each action), so we multiply each classification's equation by its occurrence number and divide by the number of actions (dims!*2**dims). Inside the classification, it works in much the same way as the previous equation generator, albeit with our more nicely-formatted input.

Now, for any value of dims, we can take

(actionCounts,func)=disdyakis(dims)

and then we can test them with

print(func)
func=eval(func)
for s in range(1,17,2):
    setBoard(s,dims,l=False)
    shapeNumbers=tap(λ i: choose(dims,i)*(choose((size-1)//2,i),),range(dims+1))
    print('over',tap(bind(zip,tuple),actionCounts,shapeNumbers))
    multiplicities=tap(dot,actionCounts,shapeNumbers)
    print('size',size)
    print(polya(size))
    print(sum(map(λ d,m: m[0]*2**sum(tap(dot,d,shapeNumbers)),shapeDistincts,RLEmatrices))//actions(dims))
    print(func((size-1)//2))

(because the function returns values for ranges, not board widths), and they should each return the same value (nonequivalents(size)).

External links

• the program (it is presented here in its entirety, but explained very quickly, the program itself has many commented-out lines for debugging, from which you may benefit if you don't understand all of the explanations here :-)