Pólya enumeration theorem applier for closed forms for A054247 in arbitrarily many dimensions

dimensional_INT_enumerator.py

#initially made for confirming there to be 1426144 INT transitions in 3D for https://conwaylife.com/wiki/Isotropic_non-totalistic_cellular_automaton (which was computed faster by hand (assumedly with the Pólya enumeration theorem), but twice with different results), now finds closed forms with the enumeration theorem (for hypercubes and orthoplexes)
from operator import __add__,__mul__#,__call__ #not sure how else to do this for call or how to have multi-input add
from functools import reduce #the land before time
from itertools import starmap,accumulate,groupby,product,combinations,chain,pairwise,zip_longest,islice,count #we manipulate iterables
from math import isqrt,lcm,floor,ceil
ceilsqrt=(lambda x: (lambda s: s+(s**2<x))(isqrt(x)))

A002024=(lambda n: isqrt(n<<3)+1>>1)
A002260=(lambda n,b=False: (lambda s: (lambda o: (o,s) if b else o)(n-s*(s-1)//2))(A002024(n))) #1-indexed antidiagonal coordinates
A003056=(lambda n: isqrt((n<<3)+1)-1>>1)
A002262=(lambda n,b=False: (lambda s: (lambda o: (o,s-o) if b==2 else (o,s) if b else o)(n-s*(s+1)//2))(A003056(n))) #0-indexed antidiagonal coordinates

from fractions import Fraction
stratrix=(lambda m,dims=2,strict=True,hidezero=False: (lambda m: '\n'.join((lambda s: (lambda c: starmap(lambda i,r: (' ' if i else '(')+(','+'\n'*(dims==3)).join(starmap(lambda i,s: ' '*(c[i]-len(s))+s,enumerate(r)))+(',' if len(m)+~i else ')'),enumerate(s)))(tap(lambda c: max(map(lambda t: bool(t) and len(t),c)),zip_longest(*s))))(tap(taph(lambda f: stratrix(f,2,strict) if dims==3 else str(f) if f else ' '),m))))(m if dims==2 else (m,)))
from time import time #nevermore

dbg=(lambda x,*s: (x,print(*s,x))[0]) #debug
transpose=(lambda l: zip(*l))
revange=(lambda a,b=None,c=1: range(b-c,a-c,-c) if b else range(a-c,-c,-c)) #(lambda *a: reversed(range(*a)))
dot=(lambda a,b: sum(map(__mul__,a,b)))
permute=(lambda a,b: map(lambda n: b[n],a)) #APL's from but that name is taken by Python
maph=(lambda f: lambda *i: map(f,i)) #in the convention of the hyperbolic functions, for currying (just like Haskell used to make :-)
tap=(lambda f,*i:tuple(map(f,*i)));taph=(lambda f: lambda *i:tuple(map(f,*i)));tarmap=(lambda f,*l:tuple(starmap(f,*l)))
lap=(lambda f,*i: list(map(f,*i)));laph=(lambda f: lambda *i: list(map(f,*i)));larmap=(lambda f,*l: list(starmap(f,*l)))
sap=(lambda f,*i:  set(map(f,*i)));saph=(lambda f: lambda *i:  set(map(f,*i)));sarmap=(lambda f,*l:  set(starmap(f,*l)))
redumulate=(lambda f,l,i=None: accumulate(l,f,initial=i)) #I did (lambda f,l,i=None,accum=False: (accumulate(l,f) if accum else reduce(f,l)) if i==None else (accumulate(l,f,initial=i) if accum else reduce(f,l,i))) but then I didn't use the reduce mode (I don't suppose it comes up much)
bind=(lambda *f: lambda *a: reduce(lambda a,f: (lambda a,i,f: (f(a) if i else f(*a)))(a,*f),enumerate(f),a))
construce=(lambda f,l,i=None: reduce(lambda a,b: f(*a,b),l,i)) #reduce for constructing tuple, typically with range (you cannot unpack parameters directly in lambda arguments)
def shortduce(f,l=None,i=None,o=None,b=None): #with different function depending on shortcut (second element is used only as whether to proceed)
    if i==None: i=next(l)
    i=(i,True)
    if l==None:
        while True:
            if i[1]: i=f(i[0])
            else: return(i[0] if b==None else b(i[0]))
    else:
        for j in l:
            if i[1]: i=f(i[0],j)
            else: return(i[0] if b==None else b(i[0]))
    return((lambda f,i: i if f==None else f(i))(o if i[1] else b,i[0]))
#(funcxp,expumulate)=map(lambda f: eval("lambda f,l: lambda i: "+f+"(lambda x,i: f(x),range(l),i)"),("reduce","redumulate")) #unfortunately has overhead
funcxp=(lambda f,l: lambda i: reduce(lambda x,i: f(x),range(l),i)) #short for funcxponentiate
consxp=(lambda f,l: lambda i: reduce(lambda x,i: f(*x),range(l),i)) #short for consxponentiate
expumulate=(lambda f,l: lambda i: accumulate(range(l),lambda x,i: f(x),initial=i)) #inlined, expumulate(f,l)(i) is equivalent to map(lambda n: funcxp(f,n)(i),range(l))
ORsum=(lambda l: reduce(int.__or__,l,0))
moddiv=(lambda a,b: divmod(a,b)[::-1])
fact=(lambda n: reduce(int.__mul__,range(1,n+1)) if n else 1)
choose=(lambda n,*k: (lambda n,*k: fact(n)//reduce(int.__mul__,map(fact,k)) if all(map(lambda k: 0<=k,k)) else 0)(n,*k,n-sum(k))) #more optimised closed form exists but not used too much
#product=(lambda *args,repeat=1: reduce(lambda result,pool: (x+(y,) for x in result for y in pool),lap(tuple,args)*repeat,((),))) #itertools' definition
decompose=(lambda n,l=None: (n>>i&1 for i in range(n.bit_length() if l==None else l)) if type(n)==int else chain(*n)) #not sure whether to do this or (lambda n: funcxp(lambda t: (lambda t,m,d: (t+(m,),d))(t[0],*moddiv(t[1],2)),n.bit_length())(((),n))[0]) #not related to compose
recompose=(lambda i: reduce(int.__or__,(k<<j for j,k in enumerate(i))))
mapdec=(lambda dims: expumulate(lambda l: (lambda i: (0,)*i+(1,)+l[i+1:])(l.index(0)),2**dims-1)((0,)*dims)) #map decompose, equivalent to (lambda dims: map(lambda n: decompose(n,dims),range(2**dims)))

def firstchoose(k,i): #n such that choose(n,k)>=i
    n=k;c=1
    while c<=i:
        n+=1;c=c*n//(n-k)
    return(n)

#bindex=(lambda i,n: reduce(lambda m,i: (lambda s,m,i,b: (s,m-1) if b else (s+choose(n+~i,m),m))(*m,*i),enumerate(decompose(i,n)),(0,i.bit_count()-1))[0]) #from https://kociemba.org/math/UDSliceCoord.htm
#bindex=(lambda i: reduce(lambda m,i: (lambda s,m,i,b: (s+choose(i,m),m) if b else (s,m+1))(*m,*i),enumerate(decompose(i)),(0,-1))[0]) #I played with it a bit, this one seems to return the sorted rank (ie. A079071)
bindex=(lambda i: reduce(lambda m,i: (lambda 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]) #longer version but more efficient (without choose)

class lexbin:
    def __init__(self,m,n=0):
        self.m=m;self.n=self.m
        self.length=choose(n,m)
    __len__=(lambda self: self.length)
    lexinc=(lambda self,n: (lambda u: ((n^n+u)//u>>2)+n+u)(n&-n)) #from OEIS (minorly beats (lambda n: (lambda t: t|((t&-t)//(n&-n)>>1)-1)((n|n-1)+1)) from Stanford I think)
    index=(lambda self,i: bindex(i)) #doesn't consider self, index() is the same for all lexbins
    getter=(lambda self,i: construce(lambda o,m,i,n: (o<<1,m,i) if i<choose(n,m) else (o<<1|1,m-1,i-choose(n,m)),revange(firstchoose(self.m,i)),(0,self.m,i))[0])
    __getitem__=(lambda self,i: expumulate(self.lexinc,i.stop-i.start-1)(self.getter(i.start)) if type(i)==slice else self.getter(i))
    __iter__=(lambda self: expumulate(self.lexinc,choose(self.n,self.m)-1)((1<<self.m)-1)) #not __getitem__ call because more efficient initialisation
class hexer:
    def __init__(self,r):
        self.r=r
        self.cells=3*self.r*(self.r+1)+1
    __len__=(lambda self: self.cells)
    #storing width as range (to remain consistent with A003215), ie. ranges 1 and 2:
    '''(1,2) (2,2)     5 6
    (0,1) (1,1) (2,1) 2 3 4
       (0,0) (1,0)     0 1
          (2,4) (3,4) (4,4)         g h i
       (1,3) (2,3) (3,3) (4,3)     c d e f
    (0,2) (1,2) (2,2) (3,2) (4,2) 7 8 9 a b
       (0,1) (1,1) (2,1) (3,1)     3 4 5 6
          (0,0) (1,0) (2,0)         0 1 2'''
    #row widths are tuple(range+1+y for y in range(range+1))+tuple(range*2-y for y in range(range))
    #if y<=range: #index<(∑_{y=0}^{self.r} (self.r+1+y)=(self.r+1)*(3*self.r+2)//2):
        #index=y*(y+2*self.r+1)//2+x
        #y=(sqrt(8*index+(2*self.r+1)**2)-3)//2-self.r
        #x=index-(y+1)*(y+2*self.r+2)//2
    #else:
        #sum(range+1-y for y in range(range))    =∑_{y=self.r+1}^{2*self.r} (self.r*2-y)  =self.r*(self.r-1)//2
        #sum(range+1-y for y in range(range+1,y-1))=∑_{y=self.r+1}^{y-1} (self.r+1-y)=(self.r-y-1)*(self.r-y+2)//2=index-(self.r+1)*(3*self.r+2)//2 (where x=0)
        #index=(y-self.r)*(y-5*self.r-1)//2+(self.r+1)*(3*self.r+2)//2+x
        #     =((y-self.r)*(y-5*self.r-1)+(self.r+1)*(3*self.r+2))//2+x
    #width=(lambda y: self.r+1+y if y<=self.r else self.r+1-y)
    #__index__=(lambda x,y: y*(y+2*self.r+1)//2+x if y<=self.r else (self.r-y-1)*(self.r-y+2)//2+2*(y-1-self.r)+(self.r+1)*(3*self.r+2)//2+x)
    __iter__=(lambda self: ((x,y) for y in range(2*self.r+1) for x in (range(self.r+1+y) if y<=self.r else range(y-self.r,2*self.r+1))))
    __index__=(lambda self,c: (lambda x,y: (y*(y+2*self.r+1)//2 if y<=self.r else ((self.r+1)*(3*self.r+2)+(y-self.r-1)*(5*self.r-y))//2)+x)(*c))
    #print(tuple(((x,y),self.index(x,y)) for y in range(2*self.r+1) for x in (range(self.r+1+y) if y<=self.r else range(y-self.r,2*self.r+1))),len(tuple(1 for y in range(2*self.r+1) for x in (range(self.r+1+y) if y<=self.r else range(y-self.r,2*self.r+1)))),3*self.r*(self.r+1)+1)
    #sum(1 for y in range(self.r+1) for x in (range(self.r+1+y)))=∑_{y=0}^{self.r} (∑_{x=0}^{self.r+y} (1))=(self.r+1)*(3*self.r+2)//2
    #(1 for y in range(self.r+2,2*self.r+1) for x in (range(y-self.r,2*self.r+1)))=∑_{y=self.r+2}^{y} (∑_{x=y-self.r}^{2*self.r} (1))=(y-self.r-1)*(5*self.r-y)//2, i=(y-self.r-1)*(5*self.r-y)//2+(self.r+1)*(3*self.r+2)//2
    #print(tuple(self.index(x,y) for y in range(2*self.r+1) for x in (range(self.r+1+y) if y<=self.r else range(y-self.r,2*self.r+1))))
    #print(tuple(len(tuple(1 for y in range(2*self.r+1) for x in (range(self.r+1+y) if y<=self.r else range(y-self.r,2*self.r+1)))) for self.r in range(8)))
    __getitem__=(lambda self,i: (lambda y: (i-y*(y+2*self.r+1)//2,y))((isqrt(8*i+(1+2*self.r)**2)-1)//2-self.r) if i<(self.r+1)*(3*self.r+2)//2 else (lambda y: (i+1-((self.r+1)*(3*self.r+2)+(y-self.r-1)*(5*self.r-y))//2,y))(3*self.r+(3-ceilsqrt(9-8*i+28*self.r*(1+self.r)))//2))
    #print(tap(self.__getitem__,range(len(self))))
    #print(tuple(((x,y),self[i]) for i,(x,y) in enumerate(coords(self.r))))
class permutation: #from quaternion_solids
    __repr__=(lambda p: 'permutation('+(','*(len(p)>9)).join(map(str,p.internal))+')')
    __iter__=(lambda p: iter(p.internal))#(lambda p: chain(iter(p.internal),count(len(p))))
    __len__=(lambda p: len(p.internal))
    __getitem__=(lambda p,i: p.internal[i] if type(i)==slice or type(i)==int and i<len(p) else i) #islice(iter(p),i)
    __add__=(lambda a,b: permutation(a.internal+tap(lambda i: i+len(a),b.internal)))
    conjugate=(lambda p: permutation(map(p.index,range(len(p)))))
    #__pow__=(lambda p,n: p.conjugate() if n==-1 else reduce(lambda r,i: p*r,range(n-1),p) if n else range(len(p)))
    #__pow__=(lambda p,n: (lambda n: p.conjugate()**-n if n<0 else reduce(lambda r,i: p*r,range(n-1),p) if n else range(len(p)))(lambda o: ((lambda n: n-o*(n>o>>1))(n%o))(order(p))))
    def __pow__(p,n): #any asymptotically faster than this would require factorising n to enact recursively, probably
        c=p.cycles()
        '''m=tap(lambda c: (lambda e: c[-e:]+c[:-e])(n%len(c)),c)
        o=lap(lambda _: None,p) #do not multiply lists by integers (worst mistake of my life)
        for c,m in zip(c,m):
            for c,m in zip(c,m):
                o[c]=m'''
        o=lap(lambda _: None,p)
        for c in c:
            #for c,m in zip(c,(lambda e: c[-e:]+c[:-e])(n%len(c))):
            for c,m in zip(c,(lambda e: c[e:]+c[:e])((lambda o: (lambda n: n-o*(n>o>>1))(n%o))(len(c)))): #very elegant I think
                o[c]=m
        return(permutation(o))
    comp=(lambda a,b: (a,permutation(b.internal+tuple(range(len(a)-len(b))))) if len(a)>=len(b) else permutation.comp(b,a))
    __eq__=(lambda a,b: __eq__(*map(tuple,permutation.comp(a,b))))
    __gt__=(lambda a,b: __gt__(*permutation.comp(a,b)))
    __lt__=(lambda a,b: __lt__(*permutation.comp(a,b)))
    __ge__=(lambda a,b: __ge__(*permutation.comp(a,b)))
    __le__=(lambda a,b: __le__(*permutation.comp(a,b)))
    __mul__=(lambda a,b: permutation(permute(a,b)))
    __rmul__=__mul__
    index=(lambda p,i: p.internal.index(i))
    def __init__(p,t): #my feeling when it cannot be a lambda :-(
        p.internal=((lambda p: reduce(lambda t,i: ((len(p)+~t[1].pop(i),)+t[0],t[1]),p,((),list(range(len(p)))))[0])(shortduce(lambda t: (lambda m,d: (((m,)+t[0],d,t[2]+1),d))(*moddiv(t[1],t[2])),i=((),t,1))[0]) if type(t)==int else tuple(t))
    __int__=(lambda p: reduce(lambda r,i: (r[0]+(len(p)+~(i[1]+sum(map(i[1].__lt__,r[1]))))*fact(len(p)+~i[0]),r[1]+(i[1],)),enumerate(p[:0:-1]),(0,()))[0])
    
    def cycles(p,o=0): #len(cycles) if o==2 else (oscillatory period) if o else (cycles themselves) (idea to use sets is from https://stackoverflow.com/a/75823973 :-)
        if p.internal==(0,): return(0 if o==2 else 1 if o else [])
        pi={i: p for i,p in enumerate(p[:len(p)])}
        cycles=(o!=2 if o else [])
        while pi:
            nexter=pi[next(iter(pi))] #arbitrary starting element
            cycle=(not(o) and [])
            while nexter in pi:
                if o: cycle+=1;curr=nexter;nexter=pi[nexter];del(pi[curr]) #a little bit of tessellation (very suspicious)
                else: cycle.append(nexter);nexter=pi[nexter];del(pi[cycle[-1]])
            if o==2: cycles+=1
            elif o: cycles=lcm(cycles,cycle)
            else: cycles.append(cycle) #inelegant (however I am not quite deranged enough for eval('cycles'+('+=1' if o==2 else '=lcm(cycles,cycle)' if o else '.append(cycle)')) :-)
        return(cycles)
    order=(lambda p: p.cycles(1))
    maxder=(lambda p: A000793(len(p)))
    modder=(lambda p: A003418(len(p))) #could be used instead of order, perhaps (depending)
    parity=(lambda p,b=None: len(p)^p.cycles(2)&1 if b==None else permutation.parity(tap(p.index,b)))

try:
    from sympy import divisors#factorint
    factorise=lambda m: tuple(divisors(m)[1:]) #tuple(factorint(m).keys())+(m,)
except:
    factorise=lambda m: tuple(filter(lambda n: not m%n,range(1,m//2+1)))+(m,) #this one is quite terrible in comparison (many trivial optimisations) but the program doesn't ever have to factorise such large integers

extend=lambda i,n: chain(i,range(len(i),n))

#similar explicit program for permutable hypercubes (finds minimum divisor that produces an oscillatory subperiod, the number with period p is a multiple of p so dividing sum of n//p's by n produces output)
permutatenumerator=lambda n,k: (2**k if k<2 else sum(starmap(lambda i,j: (lambda f: (lambda p: (lambda d: 2**(sum(map(lambda c: 1 if p==1 else max(filter(lambda f: all(map(lambda i: d[(p//f+i)%p][c]==d[i%p][c],range(p))),(1,)+factorise(p))),range(k**n)))//p))(tuple(expumulate(f,p-1)(tuple(range(k**n))))))((lambda g: g(g))(lambda g: lambda i,t: i if i and t==tuple(range(k**n)) else g(g)(i+1,f(t)))(0,tuple(range(k**n)))))(lambda t: tap(lambda c: t[sum(starmap(lambda o,p: permutation(j[o])[c//k**o%k]*k**p,enumerate(extend(permutation(i),n))))],range(k**n))),product(range(fact(n)),product(range(fact(k)),repeat=n))))//(fact(n)*fact(k)**n)) if n else 2
#print(tuple(map(lambda n: permutatenumerator(*A002262(n,2)),range(20))))
'''
(2,2,  2,     2,  2,2,2,
 1,2,  3,     4,  5,6,
 1,2,  6,    26,192,
 1,2, 22,111618,
 1,2,402,
 1,2,
 1,
'''

strate=(lambda state: ''.join(map(' o'.__getitem__,(decompose(state) if bitwise else state)))) #like str but for states
def printState(state):
    print(''.join( ('\n'+'\n'.join(''.join('o' if s else ' ' for s in map(((lambda r: s>>r&1) if bitwise else state.__getitem__),decompose(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)))
                  if dims>=3 else
                   (('-' if i%size else '=')*size+'\n'+'\n'.join(''.join('o' if s else ' ' for s in map(((lambda r: s>>r&1) if bitwise else state.__getitem__),range((i*size+j)*size,(i*size+j+1)*size))) for j in range(size))+'\n' for i in range(size**(dims-2)))))
actions=(lambda dims: fact(dims)*2**dims)
convert=(lambda m,i,a: (lambda 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): #size, dimensions, mode (0 for hypercube, 1 for orthoplex, 2 for hexagon (not infinite family)), bitwise, lambdas
    global size,dims,mode,boardCells,bitwise,unravel,matrices,exps,matrixPeriods,lambdas
    size=s
    dims=d
    mode=o
    boardCells=size**dims
    if l:
        bitwise=b
        unravel=True
        if mode==2:
            hexer(size-1)
        else:
            matrindex=(lambda 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))) #equivalent to matrices.index
            matrices=sorted((tuple((j,n>>i&1) for i,j in enumerate((lambda t: tuple(reduce(lambda 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(lambda m,i: moddiv(m[1],i),range(dims,1,-1),(0,m)))))))) for m in range(fact(dims)) for n in range(2**dims)),key=matrindex) if dims else [] #will optimise later
        compose=(lambda a,b: tap(lambda i: (lambda j,k: (j,k^i[1]))(*b[i[0]]),a))
        #print(lap(len,(matrices,lap(lambda i: matrindex(compose(matrices[0],i)),matrices))))
        #print('\n'.join(map(str,zip(matrices,lap(lambda i: compose(matrices[0],i),matrices),lap(matrindex,matrices)))))
        exps=tap(lambda m: lap(matrindex,expumulate(lambda i: compose(m,i),lcm(4,fact(dims)))(matrices[0])),matrices)
        matrixPeriods=tap(lambda e: e[1:].index(0)+1,exps)
        #print(matrixPeriods)
        lambdas=tuple(eval("lambda s: "+( ( '|'.join((lambda l: map((lambda e: ''.join("(s&("+'|'.join(starmap(lambda j,ind,d: '1<<'+str(ind),e[1]))+'))')+(lambda d: ">>"+str(d) if d>0 else '' if d==0 else "<<"+str(-d))(e[1][0][2])),l))(map(lambda x: (x[0],tuple(x[1])),groupby(sorted(((lambda j,ind: (j,ind,ind-j))(j,eval('+'.join(map(lambda i: convert(m,i,str(j)),range(dims))))) for j in range(boardCells)),key=(lambda j: j[2])),key=(lambda j: j[2])))))
                                 if unravel else #bitwise unravel tries to group together common shifts (outperformed by non-bitwise unravel on my machine for 3*3*3 cube, not sure about very large polytopes but you ought not to be using bitwise mode for those anyway)
                                  "ORsum("*(not unravel)+('|' if unravel else ',').join((lambda j,ind: '(s&1<<'+str(ind)+')'+(">>"+str(ind-j) if ind>j else '' if ind==j else "<<"+str(j-ind)))(j,eval('+'.join(map(lambda i: convert(m,i,str(j)),range(dims))))) for j in range(boardCells))+")"*(not unravel))
                               if bitwise else
                                ( '['+','.join('s['+str(eval('+'.join(map(lambda i: convert(m,i,str(j)),range(dims)))))+']' for j in range(boardCells))+']'
                                 if unravel else
                                  '[s'+'['+'+'.join(map(lambda i: convert(m,i,'i'),range(dims)))+"] for i in range("+str(boardCells)+")]"))) for m in matrices)
        '''print('\n'.join("lambda s: "+('|'.join((lambda l: map((lambda e: ''.join("(s&("+'|'.join(starmap(lambda j,ind,d: '1<<'+str(ind),e[1]))+'))')+(lambda d: ">>"+str(d) if d>0 else '' if d==0 else "<<"+str(-d))(e[1][0][2])),l))(map(lambda x: (x[0],tuple(x[1])),groupby(sorted(((lambda j,ind: (j,ind,ind-j))(j,eval('+'.join(map(lambda i: convert(m,i,str(j)),range(dims))))) for j in range(boardCells)),key=(lambda j: j[2])),key=(lambda j: j[2]))))))+'\n'+
                           "lambda s: "+('|'.join((lambda j,ind: '(s&1<<'+str(ind)+')'+(">>"+str(ind-j) if ind>j else '' if ind==j else "<<"+str(j-ind)))(j,eval('+'.join(map(lambda i: convert(m,i,str(j)),range(dims))))) for j in range(boardCells)))+'\n'+'\n'.join(' '*eval('+'.join(map(lambda i: convert(m,i,str(j)),range(dims))))+'o' for j in range(boardCells)) for m in matrices)) #for testing purposes '''
#setBoard(4,4,o=0,b=False,l=True)
#print('\n'.join(map(str,matrices)))
permutationState=(lambda key=None: lap(lambda i: (i+1)*(key(i) if key else (sum(abs(x[0]-(size-1)/2) for x in tuple(expumulate(lambda a: moddiv(a[1],size),dims)((0,i)))[1:])<=(size+1-size%2)/2) if mode else True),range(boardCells)))
#symmetries=[]
def shapedec(test=None): #like mapdec but for binary states of input shape
    nextJ=(lambda 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=(lambda: range(2**boardCells) if bitwise else shapedec() if mode else mapdec(boardCells))

def nonequivalents(size,dims,mode=False):
    setBoard(size,dims,mode)
    distincts=0
    interval=2**max(16,boardCells-10)
    t=time()
    #print(len(tuple(generateStates())))
    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((lambda f: f(state)),lambdas))))))
            t=time()
        #print('states\n'+'\n'.join(map((lambda f: strate(f[1](state))+' '+str(f[0])),enumerate(lambdas))))
        #symmetries.append(len(set(map(tuple,lap((lambda f: f(state)),lambdas)))))
        distincts+=(state==min(map((lambda f: f(state)),lambdas)))
        '''if state==min(map((lambda f: f(state)),lambdas)):
            printState(state)'''
    #print(set(symmetries))
    return(distincts)

def cycleLengths(size,dims,mode=False):
    setBoard(size,dims,mode)
    test=permutationState()
    #printState(test)
    factors=tap(factorise,matrixPeriods) #I wonder how this could be optimised
    #subperiods=tuple(tuple(sum(map(((lambda a,b: a==b!=0) if mode else int.__eq__),test,funcxp(l,a)(test))) for a in f) for l,f in zip(lambdas,factors))
    subperiods=tuple(tuple(sum(map(((lambda a,b: a==b!=0) if mode else int.__eq__),test,lambdas[exps[i][a]](test))) for a in f) for i,f in enumerate(factors))
    #print(lap(lambda a,b: tuple(zip(a,b)),factors,subperiods))
    #print('returning',size,dims) #useful when it is slow to act as a loading bar with gradual slowdown of results (albeit only polynomial)
    return((lambda t: lap(lambda t: reduce(lambda r,t: r+((t[0],t[1]-sum(i[1] for i in r if not t[0]%i[0])),),t,()),t))(tap(lambda a,b: tuple(zip(a,b)),factors,subperiods)))
#print('attention',lap(lambda s: nonequivalents(s,dims,mode),range(1,4+2*mode)))
octahedronCells=(lambda n: (n**3+5*n)//6 if n%2 else (n**3+8*n)//6+n**2);lowerBound=(lambda n: 2**octahedronCells(n)//actions(dims))
specific=(lambda cycles: int(bool(cycles) and sum(2**sum(i[1]//i[0] for i in c) for c in cycles)//len(cycles)))

fastNonequivalents=lambda s,d,m=False: 2**s if s<2 else specific(cycleLengths(s,d,m))
#print(tuple(map(lambda n: A002262(n,2),range(8))))
#print(tuple(map(lambda n: (lambda n,k: (n,k,fastNonequivalents(k,n,3)))(*(lambda n,k: (n+1,k))(*A002262(n,2))),range(8))));exit()
'''for n in range(16):
    print((lambda n,k: (n,k,fastNonequivalents(k,n,0)))(*(lambda n,k: (n+1,k))(*A002262(n,2))))'''
#print(stratrix(tuple(map(lambda n: tuple(map(lambda k: fastNonequivalents(k,n,0),range(9-n))),range(1,9)))))
#print((tuple(map(lambda n: fastNonequivalents(n,1,0),range(0,16)))))

#print([(lowerBound(n),specific(cycleLengths(n,3,mode))) for n in range(8)])
#print(lap(lambda n: specific(cycleLengths(n,2,mode)),range(12)))
def generateEquation(rle=True):
    typeCounts=tap(lambda i: lap(lambda g: (len(tuple(g[1])),g[0]),groupby(sorted(cycleLengths(i,dims,mode)))) if rle else lap(lambda c: (1,c),cycleLengths(i,dims,mode)),range(parities,parities*(dims+2)))
    #print('count done','\n'+'\n'.join(map(str,typeCounts)));exit()
    #print('\n'.join(map(lambda l: str(lap(lambda x: (x[0],lap(lambda n: n[0],x[1])),l)),typeCounts)))
    #print('\n'.join(map(lambda l: ','.join(lap(lambda x: '('+str(x[0])+',('+','.join(lap(lambda x: ''.join(map(str,x)),x[1]))+'))',l)),typeCounts)))
    #print('\n'.join(tuple(str(tuple(tuple(t[j] for t in typeCounts[i::2]) for j in range(len(typeCounts[i])))) for i in range(2,4))))
    #print('\n'.join(tuple(str(tuple((typeCounts[i][j][0],tuple(t[j][1] for t in typeCounts[i::2])) for j in range(len(typeCounts[i])))) for i in range(2,4))))
    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))
    print('togethers done')#\n','('+'\n '.join(map(lambda t: '('+'\n  '.join(map(str,t))+')',togethers))+')');exit()
    terms=tuple(tuple((f,)+tuple((i[0][0],(lambda a: tap((Fraction if mode else int),reduce(lambda a,b: lap(__add__,a,b),(tap(j.__mul__,(lambda 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(lambda i,d: (lambda a,i,d: (lambda c: ((c,)+a,tuple(k-c*j**d for j,k in enumerate(i[:-1]))))(funcxp(lambda i: tap(int.__sub__,i[1:],i[:-1]),d)(i)[0]//fact(d)))(*i,d),revange(dims+1),((),lap(lambda n: n[1],i)))[0])) for i in transpose(s)) for f,s in t) for b,t in enumerate(togethers))
    print('terms done')
    #print('('+'\n '.join(map(lambda t: '('+'\n  '.join(map(str,t))+')',terms))+')');exit()
    #print('\nif n%2 else\n'.join('+'.join(str(i[0])+'*('+'2**('+'+'.join(str(k)+'*x**'+str(j) for j,k in enumerate(i))+')' for i in t)+')' for t in terms[::-1]))
    #print('( '+'\n if n%2 else\n  '.join(map(lambda t: str(t).replace(' ',''),terms[::-1]))+')/'+str(fact(dims)*2**dims));exit()
    return('lambda n: ( '+''.join('+'.join(map(lambda t: (str(t[0])+'*')*(t[0]!=1)+'2**('+(lambda x: x if x else '0')('+'.join((lambda 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)))
setBoard(1,4,o=False)
parities=2*lcm(*range(1,1+dims*mode))
polya=generateEquation(False)
print(polya)
exit()
polya=eval(polya)
print('\n'.join(map(bind(polya,str),range(7))))
#print('\n'.join(map(lambda n: '('+str(specific(cycleLengths(n,dims,mode)))+','+str(eval(polya)(n))+')',range(7))))

enact=(lambda m,f: tap(lambda m: f[m[0]]*(-1)**m[1],m)) #enact matrix on coordinates
def disdyakis(dims):
    #vertices=tuple(mapdec(dims)) #bitwise
    weight=(lambda d: tap(lambda w: tuple(lexbin(w,d)),range(d+1)))
    weave=(lambda a,b: consxp(lambda a,b,j: (lambda j: (a|(b&1)<<j,b>>1,j+1))(next(filter(lambda i: ~a>>i&1,range(j,j+a.bit_length()+1)))),b.bit_length())((a,b,0))[0]) #weaves b's bits into the 0 bits of a
    #f=(lambda n,k: {tuple(sorted(reduce(lambda o,r: ((o[0]^((c>>o[1]&1)<<r),o[1]+1) if i>>r&1 else o),range(n),(o,0))[0] for c in range(1<<k))) for o in range(1<<n) for i in lexbin(k,n)}) #suspicious function #use o in lexbin(n-k,n) for A141905
    #print(lap(len,(f(d,i) for d in range(dims+1) for i in range(d+1))))
    #A141905=(lambda n,k: sum(choose(n,k,j) for j in range(k+1)))
    #print(lap(lambda d: lap(lambda k: A141905(d,k),range(d+1)),range(dims+1)))
    #print(tuple(A141905(d,k) for d in range(dims+1) for k in range(d+1)))
    #weights=weight(dims) #elements are lists of instances of binary numbers of weight of indexes, used for vertex partitions
    #print(lap(lambda w: lap(bin,w),weights))
    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))) #all from origins and increasing
    #print('('+',\n '.join(map(lambda a,b: str(tuple(zip(a,b))),weights,agonals))+')')
    #print('('+',\n '.join(map(lambda a: str(tuple(a)),agonals))+')')
    facets=tap(lambda 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))
    #print('n'+'\nn'.join(map(lambda x: '\n '.join(map(str,x)),zip(facets,map(lambda k: f(dims,k),range(dims+1))))))
    #print(lap(len,facets))
    #print(lap(lambda i: len(f(dims,i)),range(dims+1)))
    #A038207=(lambda n,k: choose(n,k)*2**(n-k))
    #print(lap(lambda d: lap(lambda k: A038207(d,k),range(d+1)),range(dims+1)))
    #print(tuple(A038207(d,k) for d in range(dims+1) for k in range(d+1))) #:-)
    #metaEdges=lap(,facets) #of the disdyakis polytope
    #vertices=tuple(tuple(tuple(sum(map(lambda v: v>>d&1,n))*2//len(n)-1 for d in range(dims)) for n in f) for f in facets)
    vertices=tap(lambda n: tarmap(lambda c,i: construce(lambda x,j,d: ((x+((-1)**(i>>j&1),),j+1) if d else (x+(0,),j)),decompose(c,dims),((),0))[0],product(lexbin(n,dims),range(1<<n))),revange(dims+1))
    unsignedVertices=tuple(((1,)*n+(0,)*(dims-n),) for n in range(dims,0,-1))#tuple(tap(lambda x: tuple(decompose(x,dims)),lexbin(n,dims)) for n in revange(dims+1))
    #print(tuple(transpose(lap(lambda n: lap(len,n),(facets,vertices)))))
    unsignedMode=True #do not use False for more than 3 dimensions (too slow)
    #disdyakisFacets=((((0,)*dims,),),)+tuple(redumulate(lambda x,i: lap(lambda x: tuple(lap(sum,transpose(map(lambda s: s[1],s))) for s in product(*map(lambda x: map(lambda y: (x[0],y),x[1]),x)) if all(sum(map(lambda a,b: abs(a-b),a,b))<=abs(k-l)*2**(i-2) for j,(k,a) in enumerate(s) for l,b in s[:j])),combinations(enumerate(x),i)),range(2,dims+1),vertices[:-1])) #do not use (requires function for determining whether facets are connected from centres) #map(int.__ne__,a,b))<=abs(i-j) doesn't work either
    #disdyakisFacets=tuple(redumulate(lambda x,i: tarmap(lambda e,x: tuple(tap(((lambda s: tap(sum,transpose(s[1]))) if i-2 else (lambda s: s[1])),s) for s in product(*map(lambda x: map(lambda y: (x[0],y),x[1]),x)) if all(any(starmap(tuple.__eq__,product(a,b))) if i-2 else all(map(int.__le__=,pairwise(map(lambda i: i[1],sorted(enumerate(b),key=(lambda i: a[i[0]])))))) for j,(k,a) in enumerate(s) for l,b in s[:j])),enumerate(combinations(enumerate(x),i))),range(2,dims+1),(unsignedVertices if unsignedMode else vertices)[:-1])) #condition was originally sum(map(bind(int.__sub__,abs),a,b))<=abs(k-l)
    #disdyakisFacets=((((0,)*dims,),),disdyakisFacets[0])+tap(lambda f: tap(lambda f: tap(lambda f: tap(sum,transpose(f)),f),f),disdyakisFacets[1:])
    disdyakisFacets=tap(lambda i: tarmap(lambda e,x: tuple(tap(sum,transpose(map(lambda s: s[1],s))) for s in product(*map(lambda x: map(lambda y: (x[0],y),x[1]),x))),enumerate(combinations(enumerate(unsignedVertices),i))) if i else (((0,)*dims,),),range(dims+1)) #now I should look quite the fool
    if unsignedMode:
        #print('f',disdyakisFacets)
        #print(tap(lambda d: tap(lambda d: tap(len,d),d),disdyakisFacets))
        #print('wibble','\n'.join(map(str,matrices)))
        actionCounts=tap(lambda f: tuple(reversed(tuple(decompose(map(lambda f: lap(lambda f: choose(dims,*tap(lambda x: len(tuple(x[1])),groupby(sorted(f))))*2**sum(map(bool,f)),f),f))))),disdyakisFacets) #that but without the matrices
        duplicatedFacets=tap(lambda f: tap(lambda f: tuple(set(tap(lambda f: enact(*f),product(matrices,f)))),f),disdyakisFacets) #re-sign them (not necessary)
    #print('['+']\n['.join(starmap(lambda i,f: ',\n '.join(map(lambda x: str(len(x))+' '+str(lap(lambda x: lap(lambda x: x/2**i,x),x)),f)),enumerate(disdyakisFacets)))+']')
    #print('['+']\n['.join(starmap(lambda i,f: ',\n '.join(map(bind(len,str),f)),enumerate(disdyakisFacets)))+']')
    #print(tap(len,disdyakisFacets))
    #print(tap(disdyakis,range(1,7))) #very efficient
    #print(tap(lambda x: actionCounts[x.bit_count()][bindex(x)],range(1<<dims)))
    #encoding=tuple(decompose(actionCounts))
    def period(matrix,facet):
        f=enact(matrix,facet)
        i=1
        while f!=facet:
            f=enact(matrix,f)
            i+=1
        return(i)
    print(disdyakisFacets)
    RLEmatrices=tap(lambda g: (lambda g: (len(g),*g[0]))(tuple(g[1])),groupby(sorted(zip(cycleLengths(3,dims,mode),enumerate(matrices))),lambda x: x[0])) #each one is a tuple for an equivalence class, containing its number of instances, the cycle counts under it for a size-3 hypercube (there are two with a single instance (identity and eversion) and a power of 3 cells are invariant), the index of the first instance that appears and te contents
    print('\n'.join(map(str,RLEmatrices)),'matrices out')
    #cycleTypes=tap(lambda m: tap(lambda f: tap(lambda f: m[0]*len(f)//sum(map(lambda f: 1/period(m[1],f),f)),f),duplicatedFacets),RLEmatrices) #will not work due to intermediate floats making floordiv sometimes round down (I hate the IEEE)
    cycleTypes=tap(lambda m: (m[0],m[1][-1][0],tap(lambda f: tap(lambda f: sum(map(lambda g: g[0]//tuple(g[1])[0],groupby(sorted(map(lambda f: period(m[2][1],f),f))))),f),duplicatedFacets)),RLEmatrices)
    print('\n'.join(map(str,cycleTypes)),'it begins')
    print(tap(lambda x: len(tuple(x[1])),groupby(sorted(cycleTypes))))
    setBoard(2*dims+1,dims)
    print('\n'.join(map(lambda f: str(lap(lambda f: f[0],f)),duplicatedFacets)))
    shapes=tap(lambda d: tap(lambda s: (lambda s: permutationState(lambda i: s[i]))(tuple(decompose(ORsum(map(lambda c: 1<<sum(starmap(lambda i,f: (f+dims)*size**i,enumerate(c))),s)),boardCells))),d),duplicatedFacets)
    #tap(lambda s: tap(printState,s),shapes);exit()
    cycles=(lambda p,l,shape: (lambda f: reduce(lambda f,i: f[:i+1]+lap(lambda e: (e[0]-f[i][0]*(not e[1]%f[i][1]),e[1]),f[i+1:]),range(len(f)),f))(lap(lambda i: (sum(map(lambda a,b: a==b!=0,lambdas[exps[l][i]](shape),shape)),i),factorise(p)))) #similar to some others but for specific shapes
    shapeDistincts=tap(lambda m: tap(lambda s: (tap(lambda s: sum(map(lambda x: x[0]//x[1],cycles(m[1][-1][0],m[2][0],s))),s)),shapes),RLEmatrices)
    #print('\n'.join(map(str,shapeDistincts)))
    #choosopoly=(lambda k: '('+(lambda x: x if x else '1')('*'.join(map(lambda x: '(n-'+str(x)+')',range(k))))+')/'+str(fact(k))) #the product of a length-0 sequence is 1 because uhh (∫_0^∞ (e**-t*t**x dt))(0)=1
    #choosopoly=(lambda k: '('+'+'.join(starmap(lambda i,p: ((str(p)+'*')*(p!=1)+'n'+('**'+str(i))*(i!=1))*bool(p),enumerate(reduce(lambda a,k: tarmap(lambda a,b: a-k*b,pairwise((0,)+a+(0,))),range(k),(1,)))))+')//'+str(fact(k))) #coming to a game shop near you (much better than Monopoly)
    choosopoly=(lambda k: (reduce(lambda a,k: tarmap(lambda a,b: a-k*b,pairwise((0,)+a+(0,))),range(k),(1,)),fact(k))) #coming to a game shop near you (much better than Monopoly)
    #print('\n'.join(tap(lambda k: str(k)+','+str(choosopoly(k)),range(8))));exit()
    #func='lambda n: ('+'+'.join(map(lambda d,m: str(m[0])+'*2**('+'+'.join(tap(lambda d,i: 'choose(n,'+str(i)+')'+'*'+str(sum(d)),d,range(dims+1)))+')',shapeDistincts,RLEmatrices))+')/'+str(actions(dims))
    #print(func)
    polyopoly=(lambda d: tap(sum,zip_longest(*map(lambda d,i: (lambda c: (lambda f: tap(lambda c: c*f,c[0]))(Fraction(sum(d),c[1])))(choosopoly(i)),d,range(dims+1)),fillvalue=0)))
    return(shapeDistincts,RLEmatrices,actionCounts,'lambda n: ('+'+'.join(map(lambda d,m: str(m[0])+'*2**('+'+'.join(str(c.numerator)+('*n'+('**'+str(i))*(i!=1))*bool(i)+(lambda d: ('/'+str(d))*(d!=1))(c.denominator) for i,c in enumerate(polyopoly(d)) if c)+')',shapeDistincts,RLEmatrices))+')//'+str(actions(dims)))
(shapeDistincts,RLEmatrices,actionCounts,func)=disdyakis(dims)
print(func)
print(actionCounts)
func=eval(func)
for s in range(1,17,2):
    setBoard(s,dims,l=False)
    shapeNumbers=tap(lambda 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(boardCells,multiplicities,sum(multiplicities),'darkness')
    print('size',size)
    print(polya(size))
    print(sum(map(lambda d,m: m[0]*2**sum(tap(dot,d,shapeNumbers)),shapeDistincts,RLEmatrices))//actions(dims))
    print(func((size-1)//2))

Foreword (2023-06-15): See A361870, my OEIS sequence that eventually culminated from this.

I originally made this because I read the ConwayLife wiki article for isotropic non-totalistic (INT) cellular automata (a superset of Life, remaining with a binary infinite plane and using Moore neighbourhoods (each cell considering those a king move away to determine what it will become in the next iteration), in between the rule being a function of the number of neighbours that are on (that allows 2**18 different rules) and one of the neighbourhood's state itself (that allows 2**512), instead of the neighbourhood's state under symmetry (that allows 2**102)), there is a table therein that lists various variations that have been implemented or considered (ie. rules considering von Neumann neighbourhoods (cells displaced only orthogonally) or with three states instead of two), sorted by the number of transitions (ie. configurations that rules consider nonequivalent), it didn't seem too interesting except for one for the 3D von Neumann INT rulespace, that had a footnote explaining that the value of 1426144 in the table calculated by Milo Jacquet was an upper bound, not the true value, and that user wildmyron had done calculations that returned 1426132. (Note that for whatever reason, the table excludes the considering cell's own value so each is halved there from the outputs this will yield.)

This didn't seem like too large of a number to me, and I had thought much about a similar problem when implementing my eightfold reducer for my tablebase program's cellular automaton mode, in which I had used the Pólya enumeration theorem to make efficient indexing methods (explained somewhat in this forum post, I may make a wiki article about it later). I knew that trying to use the theorem on each group action for a 3D cube, or even on equivalence classes, would lead me to possibly making errors in my answer, so I decided to instead iterate over all 2**3**3 states of the 3**3 cube of cells, for each one enacting each of the 48 group operations and checking whether the initial value was the minimum of the results under each of these, and if so incrementing the counter, and it did indeed return 1426144. (It has a bitwise mode for the purposes of simulating it for its binary cubes as quickly as possible (that was sped up about 50% by making the lambda functions it generates group common bitshifts), but one of my typical flatly-encoded list implementations seemed to work many times faster, because I added unravelling of the list comprehensions, which Python seems to like a great deal more (you may have expected that iteration would be bad for performance, but arithmetic seems more computationally expensive than I expected, compared to indexing (and this problem wasn't very parallelisable by bitwise operators, I suppose, unlike some)).)

I realised that, having made most of it as general as possible already (to avoid implementing everything in the table manually myself), it wouldn't be very much more difficult to make the number of dimensions configurable. (The group operations for hypercubes only consist of elementwise permutation and combinations of inversion, so there are n!*2**n.) I found that it computed A000616 when you keep width=2 and vary dims, because irreducible binary functions are defined as those nonequivalent under elementwise permutation and inversion. (However, if you were to want to generalise it to k-ary functions, you couldn't only look at the program's outputs with width=k and k states, you would have to account for elementwise permutation in each dimension also.)

Of course, with the program as it was, it would quickly become infeasible to compute anything meaningful in higher numbers of dimensions (being that it has O(dims!*2**(dims+width**dims)) time complexity to run with this method, so being that the 2**27 of states took me slightly under an hour, it wouldn't be able to find beyond A054247(6) (requiring 2**6**2) or A000616(5) (requiring 2**2**5)).

I was contented for a while, having resolved this mathematical issue by programming and a quite elegant implementation of lambda functions and even some eval()s, but then I began thinking some more, and realised that the Pólya enumeration theorem (as was used to compute the closed-form equation given for A054247) could be implemented for arbitrary numbers of dimensions, I will explain here (in brief, but perhaps make a wiki article about this as well).

The theorem states that for a group (a structure for a set of elements with a set of actions under which they're permuted (for hypercubes, those under which the cube itself is invariant), specifically with the property that any composition of actions is described by another action), the number of distinct states (where two are considered equivalent if an action upon one converts it to the other (and, by the property of groups, vice versa)) is the average of ((the number of states per element)**(the number of nonequivalent elements) for each action) (where equivalence is here instead only under the action being considered and its exponents). (I don't understand why the theorem should work but from this information am able to blindly apply it.)

For instance, for an n*n square of binary cells, there are eight actions.
1 identity leaves all n**2
2 orthogonal reflections leave n*(n+1)/2 if n%2 else n**2/2
2 diagonal reflections leave n*(n+1)/2
2 90º rotations leave 1+(n+1)*(n-1)/4:=(n**2+3)/4 if n%2 else n**2/4
1 inversion leaves (n**2+1)/2 if n%2 else n**2/2
For the total, we take 2** each of these and average them, for

( 2**n**2+2*2**(n*(n+1)/2)+2*2**(n*(n+1)/2)+2*2**((n**2+3)/4)+2**((n**2+1)/2)
 if n%2 else
  2**n**2+2*2**(n**2/2)+2*2**(n*(n+1)/2)+2*2**(n**2/4)+2**(n**2/2))/8

This is in fact the aforementioned closed form for A054247, and all you need to know, however for 3D you must think about it not in terms of squares and triangles with jagged edges (nice) but tetrahedra with jagged faces (not very nice at all) and will become very error-prone, and at higher numbers of dimensions still, I imagine it would become quite infeasible for human brains to reason about geometrically.

However, I experimented a small amount with this program (the residue of which remains in some commented-out fragments) and made it generate tuples of actions grouped by their numbers of cells of each period, and was able to make it sort these consistently to align them across widths of the same parity (thankfully, there are only odd and even cases to consider (because only orthogonal lines of symmetry can intersect cells xor facets, diagonals and triagonals and so on do both either way)). It does this for widths in range(2+2*(dims+1)) and discards the first two (which, having only zero and one cell respectively, don't align with the others' tuples of cycles due to only having a single element), and generates fitting polynomials for cells of each type (being that their growth rates their degree can't exceed that of the cube due to being bounded by it), converts the basis from being in terms of (n-2)/2 or (n-3)/2, and adds them up within each action type and puts them in the 2**, and averages these. It's reasonably complicated so I'm not sure of its time complexity, but it computes the case for cubes now in under three seconds instead of about 55 minutes on my machine, in particular

lambda n: ( 8*2**((1)+(-1+n)//2+(-1*n+n**3)//6)+6*2**((1)+(-1+n)//2+(-1*n+n**3)//4)+2**((1)+(-1+n**3)//2)+6*2**((n)+(-1*n+n**3)//4)+9*2**((n)+(-1*n+n**3)//2)+8*2**((n)+(-1*n+n**3)//3)+9*2**((n**2)+(-1*n**2+n**3)//2)+2**((n**3))
           if n%2 else
            12*2**((n**3)//4)+8*2**((n)//2+(-1*n+n**3)//6)+13*2**((n**3)//2)+8*2**((n)+(-1*n+n**3)//3)+6*2**((n**2)+(-1*n**2+n**3)//2)+2**((n**3)))//48

And the 4D case in about 10 minutes,

lambda n: ( 48*2**((1)+(-1+n**4)//8)+12*2**((1)+(-1+n**4)//4)+32*2**((1)+(-1+n**2)//2+(-1*n**2+n**4)//6)+12*2**((1)+(-1+n**2)//2+(-1*n**2+n**4)//4)+2**((1)+(-1+n**4)//2)+32*2**((n)+(-1*n+n**2)//2+(-1*n**2+n**4)//6)+32*2**((n)+(-1*n+n**2)//2+(-1*n+n**3)//3+(n+-1*n**2+-1*n**3+n**4)//6)+96*2**((n)+(-1*n+n**2)//2+(-1*n**2+n**4)//4)+16*2**((n)+(-1*n+n**4)//2)+12*2**((n**2)+(-1*n**2+n**4)//4)+42*2**((n**2)+(-1*n**2+n**4)//2)+32*2**((n**2)+(-1*n**2+n**4)//3)+16*2**((n**3)+(-1*n**3+n**4)//2)+2**((n**4))
           if n%2 else
            48*2**((n**4)//8)+84*2**((n**4)//4)+96*2**((n**2)//2+(-1*n**2+n**4)//6)+51*2**((n**4)//2)+48*2**((n)+(-1*n+n**2)//2+(-1*n**2+n**4)//4)+12*2**((n**2)+(-1*n**2+n**4)//2)+32*2**((n**2)+(-1*n**2+n**4)//3)+12*2**((n**3)+(-1*n**3+n**4)//2)+2**((n**4)))//384

The sequences for numbers of dimensions beyond two aren't in the OEIS, I will add them upon my chess-related sequences being approved (and perhaps one with an infinite array (read by antidiagonals) with widths in one axis and dimensions in the other).

Edit (2023-01-06): Add support for orthoplexes (with the convention of including all hypercube cells that the orthoplex lies at least partly within, and the vertices alternating between falling on vertices and cell centres of the hypercube tiling, as is explained in the Pólya enumeration theorem article on the ConwayLife wiki) and some caching of permutation matrix exponentiation, make some parts a great deal more elegant.

Edit (2023-01-10): Add mode for computing it more efficiently, for hypercubes in particular, with an algorithm (enumerating disdyakis polytope facets) that Adam P. Goucher told me about (it seems to return noninteger results for higher than 2 dimensions but now has a complete (if broken) implementation, with some interesting parts), a ConwayLife wiki article (like for the bitwise SWAR Life) is forthcoming (addendum: here it is)

Edit (2023-06-15): Add few other things, mainly a new function (permutatenumerator) for the related task of enumerating hypercube 2-colourings under action of the permutation of indexes in axes (instead of their reversal (reflection)) as well as permutation of the axes. I will make a sequence for it also, but it has the nifty property that as well as agreeing with the main program's sequence (A361870) for k=2 (a(n,2)=A361870(n,2)=A000616(n)), it satisfies a(1,k)=k+1 (each equivalence class is a bit_count), and seemingly a(2,k) is described by A007139. (It is much less efficient than the surrounding program, and contains no such nifty eval tricks because I wanted it to be a one-liner :-)