# $q$-multinomial series

I have the following code, which produces a $$q$$-multinomial coefficient, but selected randomly according to a Poisson distribution.

Consider a 3D lattice path selected uniformly at random from a cuboidal domain with random dimensions $$m1,m2,m3$$, then average over all dimensions, where each value has mean $$l1,l2,l3$$, respectively:

qmultinomial[q_, m1_, m2_, m3_] :=
Normal@Series[QFactorial[m1 + m2 + m3, q]/(
QFactorial[m1, q] QFactorial[m2, q] QFactorial[m3, q]), {q, 0,
m1 m2 + m2 m3}]
l1 = 3;
l2 = 3;
l3 = 3;
maxi = 10;
maxj = 10;
maxk = 10;
poissonqmult =
N@Series[
Total@Flatten@
Table[1/Multinomial[i, j, k] (Exp[-l1] l1^i)/i! (
Exp[-l2] l1^j)/j! (Exp[-l3] l1^k)/
k! qmultinomial[q, i, j, k], {i, 0, maxi}, {j, 0, maxj}, {k,
0, maxk}], {q, 0, 64}]; // AbsoluteTiming
ListPlot@Transpose[{Table[
i, {i, 1, Length@CoefficientList[poissonqmult, q]}],
CoefficientList[poissonqmult, q]}]


which gives

{32.9553, Null}


The problem is the code is slow, and I need to increase $$maxi,maxj,maxk$$ to finite values which are quite large (technically infinite). Is there some way I can speed it up? Perhaps using ParallelTable?

• One small note: your last line can be replaced by ListPlot@CoefficientList[poissonqmult, q]. ListPlot will automatically assume (1,2,3...) abscissae if none are provided. Unfortunately it won't make your code faster, but it certainly is more readable. – MarcoB Jun 24 '20 at 14:51

Drop the series expansion altogether and define the qmultinomial directly:

ClearAll[qmultinomial]
qmultinomial[q_, m1_, m2_, m3_] :=
QFactorial[m1 + m2 + m3, q]/(QFactorial[m1, q] QFactorial[m2, q] QFactorial[m3, q])


Then retain the rest of the code as written and you will obtain ostensibly the same result, but in 1.5 s (your original code took ca. vs. 20 s on my machine):

l1 = 3;
l2 = 3;
l3 = 3;
maxi = 10;
maxj = 10;
maxk = 10;
poissonqmult =
N@Series[
Total@Flatten@
Table[1/Multinomial[i, j, k] (Exp[-l1] l1^i)/
i! (Exp[-l2] l1^j)/j! (Exp[-l3] l1^k)/k! qmultinomial[q, i,
j, k], {i, 0, maxi}, {j, 0, maxj}, {k, 0, maxk}], {q, 0,
64}]; // AbsoluteTiming

ListPlot@CoefficientList[poissonqmult, q]

(* Out: {1.58515, Null} *)


• Ok nice, just remove the initial series. Thank you. Big help. – apkg Jun 24 '20 at 15:03