Now that I am just talking about numbers and not drawing any pictures, I think this thread may be losing interest. Also we haven't heard from the relevant experts as to why / how I am out-competing a usually-strong Mma function. But I talked breifly with Peter Luschny who found an interesting connection to the Pochammer function. Working from there and searching the OEIS more I obtained the following closed form in terms of MultiFactorial:

MultiFactorial[n_, nDim_] := 
Times[n, If[n - nDim > 1, MultiFactorial[n - nDim, nDim], 1] ]

AlternateGeneralT[n_, nLead_] := 
Table[MultiFactorial[i + nLead j + 2, nLead]/
MultiFactorial[i + 2, nLead], {i, 0, n}, {j, 0, i}] // TableForm

And the following Prints:

Grid[Partition[(Abs@GeneralT[5, #] // TableForm) & /@ Range[12], 4], 
 Frame -> All]

Grid[Partition[AlternateGeneralT[4, #] & /@ Range[12], 4], 
 Frame -> All]

Yield the same triangles

Thank you Wolfram Community Moderation Team, for showing me a little bit of respect.

Conclusion

To reiterate my points:

• What is the Mathematica Algorithm?
• Why is Mathematica going so slow?
• Can we change the Mma Source Code to use one of my definitions or something similar?
• Why aren't any of the triangles above already in OEIS?
• Is this reduction known?

I find the discovery enjoyable, and look forward to finding out just why it works this way. More work ahead . . .

Next, I'm going back to physics research, where I apply these triangles of numbers to perturbation calculations of planetary motion. Furthermore, I have another idea for applying series inversion to make it easy to compute integrals along Rotational Energy Surfaces.

HAKMEM: General Algorithms for Leading Term $x^n$

Series

$$ y = x^n + c_{n+1}\; x^{n+1}+ c_{n+2} \;x^{n+2} + c_{n+3}\; x^{n+3} + \ldots $$

Mathematica Built In

MmaExpand[m_, nLead_] := Expand[Refine[Normal[
     InverseSeries[ Series[f[x], {x, 0, m + nLead}] /. Join[
        {f[0] -> 0, (D[f[x], {x, nLead}] /. x -> 0) -> (nLead!)},
        (D[f[x], {x, #}] /. x -> 0) -> 0 & /@ Range[1, nLead - 1],
        ((D[f[x], {x, #}] /. x -> 0) -> #! c[#]) & /@ 
         Range[nLead + 1, m + nLead]]]] /. x -> y^nLead, y > 0]]

Conceptually Simple Multinomial Method

RExp[n_] := Expand[y Plus[R[0], Total[y^# R[#] & /@ Range[n]]]]

RCalc[0, _] = 1;
RCalc[n_, nLead_] := With[{basis = Subtract[Tally[Join[Range[n + nLead], #]][[All, 2]],
Table[1, {n + nLead}]] & /@ Select[IntegerPartitions[n + nLead], Length[#] > nLead - 1 &][[2 ;; -1]]},
Total@ReplaceAll[Times[-1/nLead, Multinomial @@ #, c[Total[#]],
Times @@ Power[RSet[# - 1] & /@ Range[n + nLead], #]] & /@     basis, {c[nLead] -> 1}]]

MultinomialExpand[n_, nLead_] := ( Clear@RSet; Set[RSet[#], Expand@RCalc[#, nLead]] & /@ Range[0, n];
Expand[RExp[n] /. R[m_] :> RCalc[m, nLead]])

Frobenius-Stirling Conjectural Magic

Also see Peter Luschny, Stirling-Frobenius Cycle Numbers.

SF[0, 0, m_] = 1;
SF[n_, k_, m_] := SF[n, k, m] = If[Or[k > n, k < 0], 0,
If[And[n == 0, k == 0], 1,
SF[n - 1, k - 1, m] + (m*n - 1)*SF[n - 1, k, m]
]]

GeneralPolys[n_, nLead_] := MapThread[Dot, {
Table[SF[i, j, nLead], {i, 0, n - 1}, {j, 0, i}],
Table[ x^j, {i, 0, n - 1}, {j, 0, i}]}]

GeneralT[n_, nLead_] := With[{ps = GeneralPolys[n, nLead]},
Table[(-1)^j  ps[[j]] /. {x -> i + 2}, {i, 1, n}, {j, 1, i} ]
]

Basis[n_, nLead_] := Total /@ Map[
Times[nLead^(-Length[#]) 1/Times @@ (Tally[#][[All, 2]]!),
Times @@ (# /. x_Integer :> c[x + nLead])] &, Function[{a},
Select[IntegerPartitions[n], Length[#] == a &]] /@ Range[n], {2}]

MagicExpand[n_, nLead_] := Expand[y + Dot[MapThread[Dot,
{GeneralT[n, nLead], Basis[#, nLead] & /@ Range[n]}], 
y^Range[2, n + 1]]];

Compare

In[]:= ColumnForm[MmaExpand[2, #] /. y -> (# y) & /@ Range[5]]
ColumnForm[MultinomialExpand[2, #] /. y -> (# y) & /@ Range[5]]
ColumnForm[MagicExpand[2, #] /. y -> (# y) & /@ Range[5]]

Out[]= ColumnForm[{
y - y^2 c[2] + 2 y^3 c[2]^2 - y^3 c[3], 
2 y - 2 y^2 c[3] + 5 y^3 c[3]^2 - 4 y^3 c[4], 
3 y - 3 y^2 c[4] + 9 y^3 c[4]^2 - 9 y^3 c[5], 
4 y - 4 y^2 c[5] + 14 y^3 c[5]^2 - 16 y^3 c[6], 
5 y - 5 y^2 c[6] + 20 y^3 c[6]^2 - 25 y^3 c[7]}]

Out[]= ColumnForm[{
y - y^2 c[2] + 2 y^3 c[2]^2 - y^3 c[3], 
2 y - 2 y^2 c[3] + 5 y^3 c[3]^2 - 4 y^3 c[4], 
3 y - 3 y^2 c[4] + 9 y^3 c[4]^2 - 9 y^3 c[5], 
4 y - 4 y^2 c[5] + 14 y^3 c[5]^2 - 16 y^3 c[6], 
5 y - 5 y^2 c[6] + 20 y^3 c[6]^2 - 25 y^3 c[7]}]

Out[]= ColumnForm[{
y - y^2 c[2] + 2 y^3 c[2]^2 - y^3 c[3], 
2 y - 2 y^2 c[3] + 5 y^3 c[3]^2 - 4 y^3 c[4], 
3 y - 3 y^2 c[4] + 9 y^3 c[4]^2 - 9 y^3 c[5], 
4 y - 4 y^2 c[5] + 14 y^3 c[5]^2 - 16 y^3 c[6], 
5 y - 5 y^2 c[6] + 20 y^3 c[6]^2 - 25 y^3 c[7]}]

In[]:= AbsoluteTiming[SameQ[
MagicExpand[15, #] & /@ Range[5],
MultinomialExpand[15, #] & /@ Range[5]
]]

Out[]= {4.61176, True}

In[]:= AbsoluteTiming[SameQ[
MmaExpand[10, #] & /@ Range[5],
MultinomialExpand[10, #] & /@ Range[5]
]]

Out[]= {9.79921, True}

Timing test for $n=1,2,3,4,5$

Time Vs. order of magnitude.

Not only is the mystery method the fastest, apparently it is also the most stable under change of leading term. I may try to prove this Stirling-Frobenius connection in my dissertation, but before I get into that I'm wondering: Is this algorithm already proven somewhere? Did Ferdinand Frobenius know about this? Are there Mma experts around here with an answer? Maybe Daniel Lichtblau will care to comment? Considering that Mma is going relatively slow and the optimization code is easy, I'm thinking that there's a possibility the methods here are not well known, possibly new.