login
Number of n X n X n 3-D matrices over symbol set {1,...,n} equivalent under any permutation of rows, columns, stacks or the symbol set.
2

%I #19 Jul 27 2023 12:17:39

%S 1,1,30,5887172299,1025638889935309161425800069552344,

%T 11337715575060644543768059040620852798601554512475786008052416860406653219312996

%N Number of n X n X n 3-D matrices over symbol set {1,...,n} equivalent under any permutation of rows, columns, stacks or the symbol set.

%H Philip Turecek, <a href="/A091511/b091511.txt">Table of n, a(n) for n = 0..8</a>

%F a(n) = Sum_{1*s_1+2*s_2+...=n, 1*t_1+2*t_2+...=n, 1*u_1+2*u_2+...=n, 1*v_1+2*v_2+...=n} (fixA[s_1, s_2, ...;t_1, t_2, ...;u_1, u_2, ...;v_1, v_2, ...]/ (1^s_1*s_1!*2^s_2*s_2!*... *1^t_1*t_1!*2^t_2*t_2!*... *1^u_1*u_1!*2^u_2*u_2!*... *1^v_1*v_1!*2^v_2*v_2!*...)) where fixA[...] = Product_{i, j, k>=1} ( (Sum_{d|lcm(i, j, k)} (d*v_d))^(s_i*t_j*u_k *lcm(i, j, k)/(i*j*k))).

%F a(n) asymptotic to n^(n^3)/(n!^4).

%o (Sage)

%o Pol.<x> = InfinitePolynomialRing(QQ)

%o @cached_function

%o def Z(n):

%o if n == 0: return Pol.one()

%o return sum(x[k]*Z(n-k) for k in (1..n))/n

%o @cached_function

%o def monprod(M):

%o p = Pol.one()

%o V = [m.variables() for m in M]

%o T = cartesian_product(V)

%o for t in T:

%o r = [Pol.varname_key(str(u))[1] for u in t]

%o j = [Pol(M[u[0]]).degree(u[1]) for u in enumerate(t)]

%o lcm_r = lcm(r)

%o p *= x[lcm_r]^(prod(r)/lcm_r*prod(j))

%o return p

%o @cached_function

%o def pol_isotop(n,k):

%o P = Z(n)

%o p = Pol.zero()

%o coeffs = P.coefficients()

%o mons = P.monomials()

%o C = cartesian_product(k*[mons])

%o Csorted = [tuple(sorted(u)) for u in C]

%o Cset = set(Csorted)

%o for c in Cset:

%o p += Csorted.count(c)*prod([coeffs[mons.index(u)] for u in c])*monprod(c)

%o return p

%o @cached_function

%o def rule_sub(r,m):

%o D = 0

%o for d in divisors(r):

%o try: D += d*m.degrees()[-d-1]

%o except: break

%o return D

%o def a(n,k=3):

%o P = Z(n)

%o coeffs = P.coefficients()

%o Q = pol_isotop(n,k)

%o inds = [Pol.varname_key(str(u))[1] for u in Q.variables()]

%o p = 0

%o for mon in enumerate(P.monomials()):

%o m = Pol(mon[1])

%o p += coeffs[mon[0]]*Q.subs({x[i]:rule_sub(i,m) for i in inds})

%o return p

%o # _Philip Turecek_, Jun 17 2023

%Y Cf. A091058, A091510.

%K nonn

%O 0,3

%A _Christian G. Bower_, Jan 16 2004

%E a(2) corrected by _Philip Turecek_, Jun 13 2023