Tutorials/Coding Life simulators/eightfold reducer

So you've made a Life simulator and decided to make a search program with symmetry reduction? Excellent!

Prerequisites

This will be in Python, so you should understand the functions from the previous article and the Pólya enumeration theorem (otherwise you will have to accept the ideas therein as given).

Idea

This page will be concerning the algorithm for generating canonical representations of states in n*n square boards "under action of dihedral group of the square D_4," as per OEIS sequence A054247 (which also gives us equations explaining how many exist for each width, that our program ought to be matching), ie. once you have the program's output, for every one of the 2**n**2 trivially generable states not under symmetry, only one of the eight orientations of it is in the output.

You would imagine that this article would be very short, and would explain only that if you would like to exhaustively search states in an isotropic rule, you should iterate over all 2**n**2 states and check for each whether it's the minimum of its eight orientations. However, we have a more efficient convention that has more interesting properties, that may seem at first to be unnecessarily complicated but by the end you will realise to be a more natural way.

We divide the square state into "layers," which will hereafter refer to collections of cells of the same displacement (under symmetry) from the origin, and count over each of the cells in these in binary (each layer carrying the 1 onto the beginning of the next upon overflowing, and returning to the outermost layer upon the one-carrying finishing), with the layers of lowest precedence being on the outside (from the corners to the cells at the centres of the edges, then again for the edges of the (n-2)*(n-2) square inside, repeating until terminating in the centre).

This all seems very arbitrary, however for each layer we store the cumulative intersection of the "symmetries" of those preceding it (further inwards).

To describe a subset of the group actions under which the state is invariant (ie. all cycles are comprised of cells of the same state), we define a "symmetry byte" as an integer from 0 to 255 inclusive, each of the eight bits of the binary representation of which describes whether the state is equivalent to itself under it. Where our original state is the R-pentomino,

| oo
|oo 
| o 

and the bit at an index represents the action of reflecting it in the x axis if the index's last bit is 1, then the y if its second-last is 1, then the y=x line of symmetry if the third-last is.

|index   7  |  6  |  5  |  4  |  3  |  2  |  1  |  0
|binary 111 | 110 | 101 | 100 | 011 | 010 | 001 | 000
|      -----+-----+-----+-----+-----+-----+-----+-----
|        o  | o   |  o  |   o |  o  |  o  | oo  |  oo
|repres ooo | ooo | ooo | ooo |  oo | oo  |  oo | oo 
|       o   |  o  |   o |  o  | oo  |  oo |  o  |  o 

As you can see, every transformation corresponds with an index. (This is not optimal, there are 10 possibilities for the symmetry class it falls under (C1, C2, C4, D2|, D2-, D2\, D2/, D4+, D4X and D8) that can be encoded in 4 bits, and even if we didn't want to do that, the first bit is redundant (since it cannot be different from itself under the identity) and the ones representing 90º anticlockwise and clockwise rotations are equivalent, however this way we determine and store their properties directly, that make it easier to implement (ie. the intersection of two is their bitwise AND).)

For each type of layer (containing the central cell or those orthogonally, diagonally or obliquely displaced from the centre), for each possible state we can store the symmetry byte of transformations under which it itself abides, and another under which it is minimal. (We will also store a tuple containing the lexicographic rank of each bit among the sorted symmetry byte, where tied ones share the same value.) Upon its increment, if there exist any actions that all of its predecessors (further inwards) preserve symmetry in but that it does not, it will check whether its own state under the action it has violated, is smaller than its current state, and if so will increment itself again instead of returning to the lowest-precedence (outermost) bit, skipping over a swathe of 'non-canonical' states instead of iterating over them all to determine that they're not minimal. (Note that this method requires use of its layer-precedence comparison function to compare indexes (being that the values it generates are, under this function, a strictly increasing sequence), and that simply searching lexicographically for a disparity will yield nonsense.)

This thus far is slightly interesting but, alone, probably not enough to make you decide to implement yourself. The other advantage is that each layer (and moreover, the union of all layers outside a given one) comprises a group, in that actions upon the square board will permute elements in-place, thus we can use the Pólya enumeration theorem to determine, for a given state of inner layers, the number of nonequivalent states that can be obtained by changing the outer layers. This allows, for a given input state, the number of states to be counted that have the innermost layer fixed at each value preceding that in the inputted state (and the sum for all preceding values to be computed), then fixed as its inputted value to have this computation occur again for the second-innermost layer, and so on until the layer is reached. These sums of predecessors allow, in linear time with respect to the number of layers (proportional to board squares, and the logarithm of the number of states), a state's index to be generated from its value without the use of a lookup table. By a similar method (incrementing each layer until the accumulated sum of Pólya enumeration theorem outputs exceeds an inputted integer), states can be obtained from their indexes in logarithmic time with respect to search space. Thus, we can implement methods for emulating Python's __index__ and __getitem__ methods, equivalent in time complexity to list bisection but without storing anything (if you would like to write a program of your own to conduct a census of all oscillators that emerge from the 2305843028004192256 nonequivalent 8*8 bounded boards, you can do so without paying for 16 EiB (exbibytes) of storage), and the aforementioned __next__ method that runs in constant time (constant with respect to the number of layers and thus board squares (which are a logarithm of the size of the list) because the proportion of times it will carry a 1 n layers deep decays exponentially with respect to n so the average converges to a constant).

First steps

We will implement our reducer as a class, 'reducer', within a library, 'eightfold_reducer.py', that can be trated (by programs using it) as a structure, like a list or tuple. We will be able to import the library and set a variable in a different program, say r, to reducer.reducer(), then check its length with len(r), get indexes from elements with r.index() and elements from indexes with r[index], and slices with r[l:u] (with incrementing being more efficient than indexing).

For the purposes of compatibility,, when creating reducer instances we will have an optional parameter "bitwise", determining whether to store states as integers (with margins, as was the convention in our SWAR program) or lists of booleans (also flatly encoded, row-by-row). (This doesn't cause as much complication as you may imagine, because it converts them into its layer structure then works with that.)

First, for our convenience later, we import

from itertools import groupby
from functools import reduce

All functions defined hereafter are to be created within the class (indented and within a "class reducer:" statement).

Due to the typical people who make classes in Python not being the kind to care about elegance, all variables we declare or retrieve that are specific to the class must be preceded by "self." and functions (and lambdas) must accept "self" as a parameter (due to r.f(*a) being equivalent to f(r,*a)).

This explanation will be slightly out-of-order due to functions being defined directly within the class but variables (and functions that are set to one of two definitions (unravelled for the specified board width), that require accessing local variables once upon initialisation) will be put inside an __init__ function (that Python runs for each instance of the class that is created, and passes parameters into). We will begin it with

def __init__(self,boardWidth,reduction=True,bitwise=True,cellWidth=4,cellHeight=1,OT=True,byteINT=True):
    self.bitwise=bitwise;self.cellWidth=cellWidth;self.cellHeight=cellHeight;self.cellBits=cellWidth*cellHeight;self.OT=OT;self.byteINT=byteINT
    self.symmetryReduction=reduction
    self.setBoardWidth(boardWidth,True)

(symmetryReduction, if disabled for whatever reason, makes it only iterate over states in binary as usual.)

Now, onto other functions outside __init__ (directly within the class itself).

conditionalReverse=(lambda self,reversal,toBe: reversed(toBe) if reversal else toBe)
conditionalTranspose=(lambda self,transpose,x,y: x*boardWidth+y if transpose else y*boardWidth+x)
ORsum=(lambda self,l: reduce(int.__or__,l,0))

conditionalReverse is necessary because there is no functional exponentiation in Python, it and conditionalTranspose only appear twice, and ORsum is the same as in our bitwise simulator, like sum but when we know the elements' binary representations not to intersect (or don't care but wouldn't like overflow due to it).

boardStr=(lambda self,board,inner=True,b='b',o='o',j='\n': j.join((lambda i: "".join(o if i>>self.cellWidth*j&1 else b for j in range(inner,i.bit_length()//self.cellWidth+2-inner)))(board>>(self.cellWidth*self.WIDTH*i if self.OT or self.byteINT else 3*self.WIDTH*(i*3+1)+1)&(1<<self.cellWidth*self.WIDTH)-1) for i in range(inner,self.WIDTH-inner)))
RLE=(lambda self,board,inner=True,includeRule=False: ("x ="+str(boardWidth)+", y = "+str(boardHeight)+", rule = "+RLErulestring+'\n' if includeRule else '')+''.join((lambda l,q: l*2 if self.niemiec and q==2 else l if q==1 else str(q)+l)(l,len(list(q))) for l,q in groupby(self.boardStr(board,inner,'b','o','$') if self.bitwise else '$'.join(''.join('o' if l else 'b' for l in board[b*boardWidth:(b+1)*boardWidth]) for b in range(boardWidth-1,-1,-1))))+'!')

These are the same as the ones explained in the previous article, albeit RLE also supports the list mode.

def print(self,board,chequerboard=False,delimit=True,multiple=False,spaces=False):
    part=(lambda board,b: (" "*spaces).join(("_" if letter==" " else letter+'\u0332') if (i+b)%2==0 and chequerboard else letter for i,letter in enumerate(["o" if s else " " for s in map(((lambda s: board>>self.BIAS+s//boardWidth*self.cellBits*self.WIDTH+s%boardWidth*self.cellWidth&1) if self.bitwise else board.__getitem__),range((b-1)*boardWidth,b*boardWidth))])))
    print(( "".join("["*(b==boardWidth)+(" [" if b==boardWidth else "] " if b==1 else "  ").join(part(subboard,b) for subboard in board)+(("]\n" if b==1 else "\n ") if delimit else "\n") for b in range(boardWidth,0,-1))
            if multiple else
            "["*delimit+"".join(part(board,b)+(("]\n" if b==1 else "\n ") if delimit else "\n") for b in range(boardWidth,0,-1))+"\033["+str(boardWidth)+"B"),end="")

This is similar to our previous printBoard function, except we are allowed to use the same name as a builtin function because it will be accessed as r.print() instead. Due to having been created originally for a chess program, it prints rows in reverse order (such that upwards is positive), and it allows printing multiple boards side-by-side.

def setBoardWidth(self,n,new=False):
    global boardWidth,boardSquares,halfWidth,centre,unreflectedLayerIndices,layerTypes,stateNumber
    boardWidth=n
    boardSquares=boardWidth**2
    halfWidth=(boardWidth-1)//2
    centre=(boardSquares-1)/2
    if boardWidth%2:
        centre=int(centre)
    (unreflectedLayerIndices,layerTypes)=zip(*[(self.layerIndices(x,y),self.layerType(x,y)) for x in range(halfWidth+1) for y in range(x+1)])

(The setting of centre's type depending on the parity is slightly suspicious but will become clear shortly, it is only to be used internally, and if bitwise is enabled, the variables WIDTH, COLSHIFT, WRAPSHIFT, unitNeighbourhood, lastColumn and BIAS (from the previous article) should all be set also, as well as lastRow (equivalent to WRAP_MASK but on the bottom margin row instead of inside it).)

This is where it begins to become exciting. We convert a triangular array of displacements to a flat one of layers (though each one contains 1, 4 or 8 cells depending on whether it's the centre or orthogonally or diagonally displaced, or obliquely), the variables will be explained.

layerIndices=(lambda self,x,y: ([centre] if x==0 else
                                [centre+x*p*(boardWidth if skewedness else 1) for skewedness in range(2) for p in range(-1,3,2)] if y==0 else
                                [centre+x*(xp+yp*boardWidth) for xp in range(-1,3,2) for yp in range(-1,3,2)] if y==x else
                                [centre+self.conditionalTranspose(skewedness,x*xp,y*yp) for skewedness in range(2) for xp in range(-1,3,2) for yp in range(-1,3,2)]) if boardWidth%2 else
                               ([int(centre+(x+0.5)*xp+(y+0.5)*yp*boardWidth) for xp in range(-1,3,2) for yp in range(-1,3,2)] if y==x else
                                [int(centre+self.conditionalTranspose(skewedness,(x+1/2)*xp,(y+1/2)*yp)) for skewedness in range(2) for xp in range(-1,3,2) for yp in range(-1,3,2)]))
layerType=(lambda self,x,y: (1 if x==0 else 2 if y==0 else 3 if y==x else 0) if boardWidth%2 else (3 if y==x else 0))

layerType is 1 for the centre, 2 for orthogonally-displaced cells, 3 for diagonal ones and 0 for all other (oblique) ones. layerIndices contains, for each layer, the indexes of its cells on the board. The order doesn't matter, as will be explained shortly. Inside the __init__ function, put

self.typeLengths=(8,1,4,4)
self.typeBits=(3,0,2,2)
self.instanceReflectionAllowednesses=tuple(tuple((i>>(8*l))&((1<<8)-1) for l in range(2**t)) for t,i in zip(self.typeLengths,(0xff104050203060108090c040a02080f001114010010101108101401020202010020242102202021080028010202020200301011102030201020103012202020304104440241004400404044080808040051040500111401004104440a0208050061040402002664080960440202020600101011002010211020101102002220108180840202008408808084080808080091040102001691080990810808080900a0208502022201080808840a02080a0020101100202020102010211202220020c040444080c040408080c048808080c04040440080404400808044480880804080404400808044408080804808088080f0104500203061008090c40a02080ff,0xffff,0xff50a0f0031122030c44880c0f50a0ff,0xff1144502203661188990c44a02288ff)))
self.instanceReflectionPriorities=tuple(tuple(tuple(j>>(self.typeBits[m]*l)&(self.typeLengths[m]-1) for l in range(8)) for j in i) for m,i in enumerate(((0,1635557,3959123,799370,12818092,4829193,5087304,10154067,15141658,6894081,7152192,12477633,6390865,13881882,16205448,2396160,5361783,569427,1889510,1635557,9122351,5169255,5943399,5927005,6361619,11362855,12402782,10055773,12822117,13353501,13869597,10154067,9089598,5365885,827482,1667820,4297242,7033404,7807100,5927019,14849716,13226556,16130611,10055787,12818092,13353515,13869611,13881882,2433051,5361783,3522743,569427,9089598,2433051,5419134,5361783,10892350,9032247,4755483,7430255,4297242,9089598,9380413,2433051,7687617,3438995,2892993,3959123,4555976,5437313,5945800,6186817,13511624,11364801,12138945,10315585,15403353,15419208,15935304,12477633,267475,1111390,3950485,799370,6493543,340107,1979285,1700629,12945877,9205070,4985538,8152396,6390410,13091029,15930691,5058625,296154,1111411,1886122,1893731,8881523,2838843,294984,1958698,15075242,2101320,12130272,8152417,13076266,12825386,13600097,5087304,571767,841015,1357111,1635557,6460734,2773239,3547390,569427,10589502,8966391,7611191,5796581,10720500,8827070,4554842,5361783,11415432,2635393,7950145,6056715,10980634,9166024,7810824,6187713,6620808,13229825,14003976,10316481,15141658,15420104,15936200,16205448,2102931,3177118,3951829,3700949,8624798,4646943,295425,1701973,14818517,2101761,13938372,7895692,14883484,14891093,15665804,6894081,2131610,846524,3686186,799825,8624819,2204697,7572145,3831338,15076586,14797930,6850128,10283672,6390865,12826730,15665825,6922760,4299582,841911,1358007,1373862,6461630,4638270,5412414,3265591,10590398,10831415,11339902,2634259,12818092,4297242,13338220,9089598,4757184,7396802,7687617,2892993,9346960,2434752,7744968,5884865,11415432,11358081,4757184,7687617,6620808,13254472,11415432,4757184,2895333,2907604,3423700,3959123,6721428,646604,3550659,1927499,10850196,8970115,9743811,2892993,15109395,6362753,11411330,7687617,6623148,2907618,3423714,3955098,6721442,4374433,5414360,828616,10850210,10833816,11607960,7654864,15141658,14887705,6620808,11415432,585,571767,2895333,799370,4299582,38043,296154,1635557,6623148,2102931,2361042,3959123,6390865,12818092,15141658,0),(0,0),(0,13158,52377,21760,23055,10011,36174,23055,42480,29361,55524,42480,85,13158,52377,0),(0,10011,29361,13158,36174,23055,5140,10011,55524,16705,42480,29361,52377,36174,55524,0))))

The typeLengths tells you the number of cells in a layer of each type, and typeBits the log base 2 of this. For a layer with the type t, there are 2**typeLengths[t] possible states. We store in instanceReflectionAllowednesses, for each layer type, a list of states, and within each of these a symmetry byte describing whether each reflection makes them minimal or not (which can be determined in one of many ways but we use an empty state containing only the layer and compare lexicographically). In instanceReflectionPriorities, we similarly have a list of states for layers, except each state containing the rank of the result of each reflection in the sorted list of reflections. The idea is that we will be able to begin by going through the layers, beginning with an all-1s symmetry byte, and recursively bitwise AND it with each layer's corresponding instanceReflectionAllowednesses element, until we reach the end of the layers or a point at which doing so would make the symmetry byte zero, because all actions that would make the current layer minimal are ones we've already disqualified in prior layers, so we switch to instead disqualifying those that don't lead to the resultant state's rank being minimal of those left allowed. The lists are decoded from integers that will be explained later.

layerPriorities=(lambda layer,type: self.instanceReflectionPriorities[type][self.ORsum(l<<i for i,l in enumerate(layer))])
def constructCellular(self,layers):
    if self.bitwise:
        return(self.ORsum((l>>i&1)<<(j%boardWidth*self.cellWidth+j//boardWidth*self.COLSHIFT+self.BIAS) for k,l in zip(unreflectedLayerIndices,layers) for i,j in enumerate(k)))
    else:
        output=[False]*boardSquares
        for k,l in zip(unreflectedLayerIndices,layers):
            for i,j in enumerate(k):
                output[j]=l>>i&1
        return(tuple(output))
setLayers=(lambda self,state: (lambda layers: (layers,[self.ORsum(m<<i for i,m in enumerate(l)) for l in layers]))([list(map((lambda l: state>>self.BIAS+l%boardWidth*self.cellWidth+l//boardWidth*self.COLSHIFT&1) if self.bitwise else state.__getitem__,l)) for l in unreflectedLayerIndices]))

layerPriorities converts an input layer to a symmetry byte (with our little-endian indexing convention), and constructCellular converts a list of layers back into a state. setLayers converts a state to two lists, layers and layerNumbers (the latter being only generating an integer for indexing in the instance lists, the first step of layerPriorities, for our convenience). Now, onto the reflection functions.

Firstly, functions to change indexes to correspond with reflections through a single axis (where x is 0, y is 1 and the x=y line of symmetry is 2). These will go in the __init__ function, as you will see the reason for shortly.

self.positionReflect=(lambda position,axis: boardWidth+~position+position//boardWidth*boardWidth*2 if axis==0 else position%boardWidth+(boardWidth+~position//boardWidth)*boardWidth if axis==1 else position//boardWidth+(position%boardWidth)*boardWidth if axis==2 else position)

This is all mostly simple, the index is moduloed by the width to get the x position and floordiv'd to get the y. We can do a small reduction, however, in x reflections we bitwise NOT the position and add it to board width (causing it to be reflected to a negative location 'beneath' the board) and add twice its y position (remember, floordiv and multiplication are of the same precedence so are evaluated left-to-right), achieving it in only a single modular operation instead of two as you may expect. Now, for enacting compositions of the three (remember, ordered x, y, x=y).

self.compoundPositionReflect=(lambda position,axes,reverse=False: ( ( ( (position+1)*boardWidth+~((boardSquares+1)*(position//boardWidth)))
                                                                       if axes[reverse] else
                                                                        ( boardWidth*(boardWidth+~position)+(1+boardSquares)*(position//boardWidth)))
                                                                     if axes[0]^axes[1] else
                                                                      ( ( ((1-boardSquares)*~(position//boardWidth)-position*boardWidth))
                                                                       if axes[0] else
                                                                        position*boardWidth-(boardSquares-1)*(position//boardWidth)))
                                                                 if axes[2] else
                                                                  ( ( (boardSquares+~position)
                                                                     if axes[0] else
                                                                      position+boardWidth*(boardWidth+~(2*(position//boardWidth))))
                                                                   if axes[1] else
                                                                    ( boardWidth*(1+2*(position//boardWidth))+~position
                                                                     if axes[0] else
                                                                      position)))

This is similar except with a tree of ternary operators instead. (Some may be reducible but it isn't used in this program and is only included here for completeness, so isn't a concern.) Now, enough with indexes, onto reflecting entire boards!

As mentioned earlier, we will generate functions with their iteration unravelled, like the manifolds with Möbius edges in the simulator except unable to do logarithmic-time bitwise reversal (due to it requiring each row overflow into a footprint of the next power of 2 up from its number of bits, thus causing overlap (a problem that the simulator circumvents by having there be a margin of at least one row between each edge)), so we will have to do so in linear time (though this is easier because it doesn't require intermediate phases so lets us generate a closed form). First, for our convenience, we will make a shift function, that converts an integer to a string equivalent to a signed shift rightwards.

self.shift=(lambda i: '>>'+str(i) if i>0 else '' if i==0 else '<<'+str(-i))
self.boardReflect=eval("lambda board,axis: board if axis==-1 else "+( ''.join('|'.join(f)+(' if axis=='+str(i)+' else ')*(i!=2) for i,f in enumerate(expressions))
                                                                     if self.bitwise else
                                                                      ''.join('('+','.join('board['+str(j)+']' for j in l)+')'+(' if axis=='+str(i)+' else ')*(i!=2) for i,l in enumerate(indexes))))

The bitwise mode will intake lists of strings that, if each enacted upon the board, would need only be ORsummed to reflect it, and generates a loopless function equivalent to their ORsum (having to be 'recompiled' when the board width is changed), and the other (boolean list) one concatenates the callings of the indexes to form an expression for a new tuple. The lists themselves are generated as follows.

expressions=(("(board&"+str(self.lastColumn<<self.cellWidth*i)+')'+self.shift(self.cellWidth*(i-(self.WIDTH+~i))) for i in range(self.WIDTH)),
             ("(board&"+str(self.lastRow<<self.COLSHIFT*i)+')'+self.shift(self.COLSHIFT*(i-(self.HEIGHT+~i))) for i in range(self.HEIGHT)),
             ("(board&"+str(self.ORsum(1<<self.cellWidth*(i+(self.WIDTH+1)*j)&self.lastRow<<self.COLSHIFT*j for j in range(self.WIDTH) if i+(self.WIDTH+1)*j>0))+')'+self.shift(self.cellWidth*(-i*(self.WIDTH-1))) for i in range(1-self.WIDTH,self.WIDTH)))
indexes=((x for y in range(boardWidth) for x in range((y+1)*boardWidth-1,y*boardWidth-1,-1)),
         (x for y in range(boardWidth,0,-1) for x in range((y-1)*boardWidth,y*boardWidth)),
         (x for y in range(boardWidth) for x in range(y,y+boardSquares,boardWidth)))

indexes is simple enough, and the first two expressions are only doing the same thing as the toroidal scrolling in the simulator (ANDing with row/column masks and shifting them one at a time), the third is similar except upon diagonals instead. Now, for compositions.

self.compoundReflect=((lambda board,axes: board if axes==(0,0,0) else reduce(self.boardReflect,(i for i,a in enumerate(axes) if a),board)) if self.bitwise else eval("lambda board,axes: "+''.join('('+','.join("board["+str(x)+']' for y in self.conditionalReverse(axes[not axes[2]],range(boardWidth)) for x in self.conditionalReverse(axes[axes[2]],(range(y,y+boardSquares,boardWidth) if axes[2] else range(y*boardWidth,(y+1)*boardWidth))))+')'+(' if axes==('+','.join(map(str,axes))+') else')*(i!=7) for i,axes in ((i,tuple(i>>j&1 for j in range(3))) for i in range(8)))))

This is another unravelled function definition generated on-the-fly for the non-bitwise mode, the same as the index list except more complicated due to being more general, though for the bitwise one, for now it suffices to use only a reduce and run the masks up to thrice (because it's only queried once per call). For converting states' uninflated bitwise representations (ie. one bit per cell instead of four or eight) to their canonical ones (either an inflated integer or list of booleans), we define another function, intState.

self.intState=(lambda state: self.ORsum((state&1<<i)<<(self.BIAS+i%boardWidth*self.cellWidth+i//boardWidth*self.COLSHIFT-i) for i in range(boardSquares)) if self.bitwise else tuple(bool(state>>i&1) for i in range(boardSquares)))

Now, for the length, we use the equation from A054247 (from the Pólya enumeration theorem), ensuring we use floordivs in each exponent as well as the outermost one (because for consistency, regular division always returns floats even when there is no remainder),, or 2**boardSquares otherwise.

self.length=( ( (((2**boardSquares+2**((boardSquares+1)//2))//2+2**((boardSquares+3)//4))//2+2**((boardSquares+boardWidth)//2))//2
               if boardWidth%2 else
                ((2**(boardSquares-1)+3*2**(boardSquares//2-1))+2**(boardSquares//4)+2**((boardSquares+boardWidth)//2))//4)
             if self.symmetryReduction else
              2**boardSquares)

And elsewhere in the class we can put

__len__=(lambda self: self.length)

, allowing us to override Python's builtin length function for our class.

Algorithms

Firstly, we would like one for canonicalising states (ie. returning the single orientation of them that appears in our emulated structure) for our indexing purposes.

def symmetry(self,state,reflectionMode=0):
    if self.symmetryReduction:
        (layers,layerNumbers)=self.setLayers(state)
        newReflections=(1<<8)-1
        l=0
        while l<len(layers):
            newReflections&=self.instanceReflectionAllowednesses[layerTypes[l]][layerNumbers[l]]
            if newReflections:
                reflections=newReflections
                l+=1
            else:
                break
        while l<len(layers):
            m=min(p for i,p in enumerate(self.instanceReflectionPriorities[layerTypes[l]][layerNumbers[l]]) if reflections>>i&1)
            newReflections=reflections&self.ORsum((p==m)<<i for i,p in enumerate(self.instanceReflectionPriorities[layerTypes[l]][layerNumbers[l]]))
            reflections=newReflections
            if not reflections&reflections-1:
                break
            l+=1
        reflections=next(tuple(i>>a&1 for a in range(3)) for i in range(8) if reflections>>i&1)
    else:
        reflections=(False,)*3
    return reflections if reflectionMode==2 else (self.compoundReflect(state,reflections),reflections) if reflectionMode==1 else self.compoundReflect(state,reflections)

We keep ANDing the symmetry byte with the sets of actions under which each layer is minimal until the value it would become by doing so has none remaining, then we use the minimum rank instead, as explained. There, we instead terminate if it ANDs with itself minus 1 to yield 0 (which evaluates to False), which only powers of 2 and 0 satisfy (and 0 cannot be reached because we are checking which actions' resultant states match the minimum rank that hasn't been disqualified). Now for our function for incrementing states (providing the next one, given the current).

Note that symmetrical states will have multiple bits remaining on in their symmetry bytes (all of which reflect it to the same state), so if we wanted to receive the order of symmetry, before converting it to a reflection tuple on the last line within the symetryReduction, we could use reflections.bit_count().

def generateCellular(self,state,ind):
    if self.symmetryReduction:
        (layers,layerNumbers)=self.setLayers(state)
        branchReflections=reduce(lambda li,i: li+[(lambda li,arm,t,n: (li[-1]&(lambda pr,m: self.ORsum((p==m)<<i for i,p in enumerate(pr)))(*(lambda pr: (pr,min(p for i,p in enumerate(pr) if arm==0 or li[-1]>>i&1)))(self.instanceReflectionPriorities[t][n]))))(li,i[0],*i[1])],enumerate(zip(layerTypes,layerNumbers)),[(1<<8)-1])[1:]#[(1<<8)-1]*len(layerNumbers)
        arm=len(layerNumbers)-1
        #elbow=arm #layer at which it switches from instanceReflectionAllowednesses to self.instanceReflectionPriorities
        for i in range(ind[1]-ind[0]):
            arm=len(layerNumbers)-1
            yield(self.constructCellular(layerNumbers))
            while True:
                while arm>=0 and layerNumbers[arm]==(1 if layerTypes[arm]==1 else 2**self.typeLengths[layerTypes[arm]]-1):
                    layerNumbers[arm]=0
                    del branchReflections[-1]
                    arm-=1
                if arm<0:
                    break
                else:
                    layerNumbers[arm]+=1
                    m=min(p for i,p in enumerate(self.instanceReflectionPriorities[layerTypes[arm]][layerNumbers[arm]]) if len(branchReflections)<=1 or branchReflections[-2]>>i&1)
                    if (self.instanceReflectionPriorities[layerTypes[arm]][layerNumbers[arm]][0]==m):
                        break
            if arm<0:
                    break
            branchReflections[-1]=((1<<8)-1 if arm==0 else branchReflections[-2])&self.ORsum((p==m)<<i for i,p in enumerate(self.instanceReflectionPriorities[layerTypes[arm]][layerNumbers[arm]]))
            branchReflections+=[branchReflections[-1]]*(len(layerNumbers)+~arm)
    else:
        for i in range(ind[1]-ind[0]):
            yield(state)
            j=0
            state[0]^=True
            while not state[j] and j<boardCells-1:
                j+=1
                state[j]^=True

It begins by converting to layer form with setLayers, then sets its list, branchReflections, which is a symmetry byte for each layer, being the intersection of its predecessor's byte and that of the reflections that result in it being minimal (of those under the actions its predecessor allows). This allows us to check later (in each increment) whether the identity action results in the overall state being minimal (of rank m within the priorities), and skip over values of layers that aren't (and avoiding all iteration over their succeeding layers (which become exponentially larger as the number of layers increases), hence the skipping of swathes of non-canonical states). arm is the index of the layer the binary incrementer is carrying the 1 through, and branchReflections grows and shrinks with it (containing the layers preceding and including it), and is filled in with more copies of the last element when it is reset to the end (being that layers the 1 has been carried through will be all-zero and won't change the allowed actions).

Without symmetry reduction, it is much simpler, and only increments and carries 1s in the state's binary representation (which is faster than using a list comprehension to convert each element of a range() each time).

Note that in both cases, it uses yield statements instead of constructing and returning a list, so it will generate elements on-the-fly as they're requested instead of having to pause everything else for a potentially very expensive computation (that could consume a great deal of memory), though it means it can only be iterated over, not subscripted (unless it is converted to a list or tuple first).

polya will be a function for applying the Pólya enumeration theorem to all cells in layers succeeding the arm instead of preceding it, it can be a great deal easier than some programs because we know the numbers of cells of each layer type in each cycle.

    def polya(self,layer,reflections):
        cellTypes=[0]*4
        for l in layerTypes[layer+1:]:
            cellTypes[l]+=self.typeLengths[l]
        n=0
        s=0
        for i in range(8):
            if reflections>>i&1:
                n+=1
                s+=2**(sum(cellTypes) if i==0 else cellTypes[1]+((3*cellTypes[2]+2*cellTypes[3]+2*cellTypes[0])//4 if i==1 or i==2 else (cellTypes[2]+cellTypes[3]+cellTypes[0])//2 if i==3 else (2*cellTypes[2]+3*cellTypes[3]+2*cellTypes[0])//4 if i==4 or i==7 else (cellTypes[2]+cellTypes[3]+cellTypes[0])//4))
        return(s//n if n else 0)

For an inputted state, by running polya for each state of the innermost layer preceding its inputted value, then fixing it at the inputted one and doing the same for the next layer, and so on, and summing its results, we will obtain the index within the emulated structure, which is what the index function returns.

def index(self,state):
    if self.symmetryReduction:
        state=self.symmetry(state)
        (layers,layerNumbers)=self.setLayers(state)
        reflections=(1<<8)-1
        stateIndex=0
        l=0
        while l<len(layers):
            for n in range(layerNumbers[l]+1):
                m=min(p for i,p in enumerate(self.instanceReflectionPriorities[layerTypes[l]][n]) if reflections>>i&1)
                if self.instanceReflectionPriorities[layerTypes[l]][n][0]==m:
                    newReflections=reflections&self.ORsum((p==m)<<i for i,p in enumerate(self.instanceReflectionPriorities[layerTypes[l]][n]))
                    if n<layerNumbers[l]:
                        stateIndex+=self.polya(l,newReflections)
                    else:
                        break
            reflections=newReflections
            l+=1
        return(stateIndex)
    else:
        return(reduce(int.__or__,(1<<i&s for i,s in enumerate(state))))

Now we need one more such function containing an algorithm, for doing the inverse, getting a state from an index.

def getter(self,ind):
    if self.symmetryReduction:
        layerNumbers=[0 for l in range(len(layerTypes))]
        newReflections=reflections=(1<<8)-1
        stateIndex=0
        for l,t in enumerate(layerTypes):
            for n in range(2**self.typeLengths[t]):
                m=min(p for i,p in enumerate(self.instanceReflectionPriorities[t][n]) if reflections>>i&1)
                if self.instanceReflectionPriorities[t][n][0]==m:
                    newReflections=reflections&self.ORsum((p==m)<<i for i,p in enumerate(self.instanceReflectionPriorities[t][n]))
                    nextIndex=stateIndex+self.polya(l,newReflections)
                    layerNumbers[l]=n
                    if nextIndex>ind:
                        break
                    else:
                        stateIndex=nextIndex
            if stateIndex==ind:
                break
            reflections=newReflections
        return(self.constructCellular(layerNumbers))
    else:
        return(self.intState(ind))

The idea to keep in mind is that it's doing the same thing except incrementing each layer until doing so would make the state's index exceed the inputted value before moving on, instead of until the layer's value exceeds that in the state. Now, we can combine this with generateCellular, to make a __getitem__ function (that is called when subscripting with square brackets, and inputted an integer or slice).

def __getitem__(self,ind):
    if type(ind)==slice:
        ind=ind.indices(len(self))
        if ind[2]!=1:
            raise(ValueError("you WILL NOT use increments other than 1"))
        return(self.generateCellular(self.getter(ind[0]),ind))
    else:
        return(self.getter(ind))

We would usually be able to inline generateCellular to avoid splitting the program across more functions than is necessary, but Python seems to interpret any function that contains both yield and return statements as a generator object instead of evaluating it when the input isn't a slice.

Generating instanceReflectionAllowednesses and instanceReflectionPriorities

The reason for encoding to integers after their generation and decoding afterwards is to avoid having to run this, which would require setting boardWidth to an odd number (at least 5) so that a layer of every type may occur.

def layerBoard(layer,layerIndices):
    board=[False]*boardSquares
    for n,m in enumerate(layerIndices):
        board[m]=bool(layer>>n&1)
    return(board)
instances=[[layerBoard(n,i) for n in range(2**l)] for i,l in ((unreflectedLayerIndices[layerTypes.index(i)],l) for i,l in enumerate(self.typeLengths))]
self.instanceReflectionPriorities=[[[l.index(f) for f in j] for j,l in zip(i,[sorted(set(j)) for j in i])] for i in [[[self.compoundReflect(j,[f>>r&1 for r in range(3)]) for f in range(8)] for j in [tuple(layerBoard(n,unreflectedLayerIndices[layerTypes.index(i)])) for n in range(2**l)]] for i,l in enumerate(self.typeLengths)]]
instanceReflectionAllowednesses=[[self.ORsum((l==0)<<i for i,l in enumerate(j)) for j in i] for i in self.instanceReflectionPriorities]
def equation(name,ins,nestedness):
    print(name+"=("+(str(ins) if nestedness==0 else ",".join(map(str,ins) if nestedness==1 else (("("+",".join(map(str,i) if nestedness==2 else ("("+",".join(map(str,j))+")" for j in i))+")") for i in ins)))+")")
equation("instanceReflectionAllowednesses",[self.ORsum(l<<(8*i) for i,l in enumerate(j)) for j in instanceReflectionAllowednesses],0)
print(tuple(map(tuple,instanceReflectionAllowednesses))==tuple(tuple((i>>(8*l))&((1<<8)-1) for l in range(2**t)) for t,i in zip(self.typeLengths,[self.ORsum(l<<(8*i) for i,l in enumerate(j)) for j in instanceReflectionAllowednesses])))
equation("self.instanceReflectionPriorities",self.instanceReflectionPriorities,3)
equation("self.instanceReflectionPriorities",tuple(tuple(tuple(j>>(self.typeBits[m]*l)&(self.typeLengths[m]-1) for l in range(8)) for j in i) for m,i in enumerate(((self.ORsum(k<<(self.typeBits[m]*l) for l,k in enumerate(j)) for j in i) for m,i in enumerate(self.instanceReflectionPriorities)))),3)
equation("self.instanceReflectionPriorities",((self.ORsum(k<<(self.typeBits[m]*l) for l,k in enumerate(j)) for j in i) for m,i in enumerate(self.instanceReflectionPriorities)),2)
print(tuple(map(lambda n: tuple(map(tuple,n)),self.instanceReflectionPriorities))==tuple(tuple(tuple(j>>(self.typeBits[m]*l)&(self.typeLengths[m]-1) for l in range(8)) for j in i) for m,i in enumerate(((self.ORsum(k<<(self.typeBits[m]*l) for l,k in enumerate(j)) for j in i) for m,i in enumerate(self.instanceReflectionPriorities)))))

layerBoard is like constructCellular but only for a single layer (keeping all other cells off), equation prints a line for setting a variable to an encoding of a nested structure (in an inelegant manner (repeating code) but it doesn't need to work for arbitrary nestedness), we pack the elements' bitwise representations in each tuple, being that we know for the allowednesses that each element (symmetry byte) is in range(1<<8) and for the priorities each element (priority of an action within a layer state) is in range(8) (hence the different nestednesses).

Demonstration

You can import it and run

(lambda r: (lambda r,l: (print(l),r.print(tuple(r[l//2:l//2+4]),multiple=True)))(r,len(r)))(reducer(8))

and it will return

2305843028004192256
[    o    [    o    [o   o    [o   o   
  oooooo    oooooo    oooooo    oooooo 
  oo  ooo   oo  ooo   oo  ooo   oo  ooo
  oooooo    oooooo    oooooo    oooooo 
  oo  oo    oo  oo    oo  oo    oo  oo 
   ooooo     ooooo     ooooo     ooooo 
  oooooo    oooooo    oooooo    oooooo 
         ] o       ]         ] o       ]

(telling you its length and the four states in the middle). Note that storing these 2305843028004192256 64-bit states optimally would require 16 EiB (exbibytes), plus a corrective factor (from the Pólya enumeration theorem) of 1146881/140737488355328ths (140 GiB) due to the symmetrical states not being divided by 8.

External Links

• the program