q = 2 F = GF(q) R = PolynomialRing(F, 'x') def Centr(p,a): s,sp,P = 0,0,1 for k in range(len(p)): sp += 2*s*p[-k-1]*(len(p)-k)+p[k]**2*k s += p[-k-1] P *= prod(q^(p[k]*a)-q^(i*a) for i in range(p[k])) return P*q^(sp*a) N = 30 Facts = [R(x^k-1).factor() for k in range(1,N+1)] for n in range(1,N+1): S = 0 for P in Partitions(n): D = {} for k in P: for p in Facts[k-1]: if p[0] in D: if p[1]>len(D[p[0]]): D[p[0]].extend((p[1]-len(D[p[0]])-1)*[0]+[1]) else: D[p[0]][p[1]-1] += 1 else: D[p[0]] = (p[1]-1)*[0]+[1] S += prod(Centr(D[p],p.degree()) for p in D)/P.aut()^2 print(n,S)