%I #16 Jun 24 2023 16:23:15
%S 1,3,5100,305727048000,7748770873210669158912000,
%T 476007332700693200670745550306381336371200000,
%U 272661655519533773844144991586798737775635133552905539740860416000000000
%N a(n) = Product_{i=1..n, j=1..n} (i^3 + j^3 + 1).
%C Product_{i=1..n, j=1..n} (1 + 1/(i + j)) = A324444(n) / A079478(n) ~ 2^(2*n + 1) / (sqrt(Pi)*n^(3/2)).
%C Product_{i=1..n, j=1..n} (1 + 1/(i^2 + j^2)) = A324443(n) / A324403(n) ~ c * n^(Pi/2), where c = A306398 * 2^(3/4) * exp(-Pi/12) * Pi^(1/4) * Gamma(3/4) = 0.36753062884677326134620846786416595535234038999313315993144237973600...
%F a(n) ~ A307209 * A324426(n).
%F a(n) ~ c * A * 2^(2*n*(n+1) + 1/4) * exp(Pi*(n*(n+1) + 1/6)/sqrt(3) - 9*n^2/2 - 1/12) * n^(3*n^2 - 3/4) / (3^(5/6) * Pi^(1/6) * Gamma(2/3)^2), where c = A307209 = Product_{i>=1, j>=1} (1 + 1/(i^3 + j^3)) = 3.504782999339728375891120570... and A is the Glaisher-Kinkelin constant A074962.
%p a:= n-> mul(mul(i^3+j^3+1, i=1..n), j=1..n):
%p seq(a(n), n=0..7); # _Alois P. Heinz_, Jun 24 2023
%t Table[Product[i^3 + j^3 + 1, {i, 1, n}, {j, 1, n}], {n, 1, 8}]
%Y Cf. A307209, A324403, A324426, A324443, A324444.
%K nonn
%O 0,2
%A _Vaclav Kotesovec_, Mar 28 2019
%E a(0)=1 prepended by _Alois P. Heinz_, Jun 24 2023