%I #67 May 31 2022 11:38:31
%S 0,0,8,251,2555,15055,63655,214918,616326,1561110,3586110,7612385,
%T 15139553,28506101,51229165,88438540,147420940,238291788,374813076,
%U 575377095,864177095,1272587195,1840775123,2619572626,3672629650,5078879650,6935344650,9360309933
%N a(n) = Sum_{1 <= i < j <= n} j^3*i^3.
%C a(n) is the sum of all products of two distinct cubes of positive integers up to n, i.e., the sum of all products of two distinct elements from the set of cubes {1^3, ..., n^3}.
%H Roudy El Haddad, <a href="https://arxiv.org/abs/2102.00821">Multiple Sums and Partition Identities</a>, arXiv:2102.00821 [math.CO], 2021.
%H Roudy El Haddad, <a href="https://doi.org/10.7546/nntdm.2022.28.2.200-233">A generalization of multiple zeta value. Part 2: Multiple sums</a>. Notes on Number Theory and Discrete Mathematics, 28(2), 2022, 200-233, DOI: 10.7546/nntdm.2022.28.2.200-233.
%H <a href="/index/Rec#order_09">Index entries for linear recurrences with constant coefficients</a>, signature (9,-36,84,-126,126,-84,36,-9,1).
%F a(n) = Sum_{j=2..n} Sum_{i=1..j-1} j^3*i^3.
%F a(n) = n*(n+1)*(n-1)*(21*n^5+36*n^4-21*n^3-48*n^2+8)/672 (from the generalized form of Faulhaber's formula).
%F From _Alois P. Heinz_, Jan 27 2022: (Start)
%F a(n) = Sum_{i=1..n} A000578(i)*A000537(i-1) = Sum_{i=1..n} i^3*(i*(i-1)/2)^2.
%F G.f.: -(x^5+64*x^4+424*x^3+584*x^2+179*x+8)*x^2/(x-1)^9. (End)
%e For n=3, a(3) = (2*1)^3+(3*1)^3+(3*2)^3 = 251.
%t CoefficientList[Series[-(x^5 + 64 x^4 + 424 x^3 + 584 x^2 + 179 x + 8) x^2/(x - 1)^9, {x, 0, 27}], x] (* _Michael De Vlieger_, Feb 04 2022 *)
%o (PARI) a(n) = sum(i=2, n, sum(j=1, i-1, i^3*j^3));
%o (PARI) {a(n) = n*(n+1)*(n-1)*(21n^5+36n^4-21n^3-48n^2+8)/672};
%o (Python)
%o def A347107(n): return n*(n**2*(n*(n*(n*(n*(21*n + 36) - 42) - 84) + 21) + 56) - 8)//672 # _Chai Wah Wu_, Feb 17 2022
%Y Cf. A000537, A000578.
%Y Cf. A346642 (for nondistinct cubes).
%Y Cf. A000217 (for power 0), A000914 (for power 1), A000596 (for squares).
%K nonn,easy
%O 0,3
%A _Roudy El Haddad_, Jan 27 2022