(* Mathematica program for A000620-A000625, Jean-Francois Alcover, Jul 18 2022, translated and adapted from N. J. A. Sloane's Maple code *) Ps = {0, 0, 0}; Pn = {1, 1, 1}; Ss = {0, 0, 0}; Sn = {0, 0, 1}; Ts = {0, 0, 0}; Tn = {0, 0, 0}; As = {0, 0, 0}; An = {1, 1, 2}; T = {1, 1, 2}; Do[AppendTo[Ps, As[[n - 1]]]; AppendTo[Pn, An[[n - 1]]]; t1 = Sum[2*T[[n - 1 - j]]*T[[j]], {j, 1, Floor[(n - 2)/2]}]; If[Mod[n, 2] == 1, i = (n - 1)/2; t1 = t1 + T[[i]]^2 - An[[i]]]; AppendTo[Ss, t1]; t2 = 0; If[Mod[n, 2] == 1, i = (n - 1)/2; t2 = An[[i]]]; AppendTo[Sn, t2]; t3 = 0; Do[Do[k = n - 1 - i - j; If[j < k, t3 = t3 + 2*T[[i]]*T[[j]]*T[[k]]], {j, i + 1, (n - 2)/2}], {i, 1, (n - 1)/3}]; t4 = 0; t5 = 0; Do[j = n - 1 - 2*i; If[j > 0 && i != j, t4 = t4 + (T[[i]]^2 - An[[i]])*T[[j]] + An[[i]]*As[[j]]; t5 = t5 + An[[i]]*An[[j]]], {i, 1, (n - 2)/2}]; t6 = 0; t7 = 0; If[Mod[n, 3] == 1, i = (n - 1)/3; t6 = (2*T[[i]] + T[[i]]^3)/3 - An[[i]]^2; t7 = An[[i]]^2]; AppendTo[Ts, t3 + t4 + t6]; AppendTo[Tn, t5 + t7]; AppendTo[As, Ps[[n]] + Ss[[n]] + Ts[[n]]]; AppendTo[An, Pn[[n]] + Sn[[n]] + Tn[[n]]]; AppendTo[T, As[[n]] + An[[n]]], {n, 4, 60}]; A000620 = Ps; A000621 = Pn; A000622 = Ss; A000623 = Ts; A000624 = Tn; A000625 = T;