The probability that, after drawing the lottery \(K\) times, each item \(i\) is obtained \(c_i\) times is \(\displaystyle \frac{K!}{\prod_{i=1}^{N}c_i!}\prod_{i=1}^{N}p_i^{c_i}\). For a fixed \(c_i\), the contribution for \(i\) can be calculated, so the following DP (Dynamic Programming) can be used to solve the original problem.

Let \(dp[x][y][z]\) be the sum, for all sequences of non-negative integers \(c=(c_1,\ldots,c_x)\) such that \(|\{i\mid c_i\neq 0\}|=y\) and \(\sum_{i=1}^{x} c_i=z\), of \(\displaystyle \frac{1}{\prod_{i=1}^{x}c_i!}\prod_{i=1}^{x}p_i^{c_i}\). What we want to find is \(dp[N][M][K]\times K!\).

By inspecting all possible values for \(c_{x+1}\), the transition from \(dp[x][y][z]\) to \(dp[x+1][*][*]\) can be done. For more details, please refer to the sample code.

Since the number of states is \(O(NMK)\) and each state has \(O(K)\) destination of transitions, the overall time complexity is \(O(NMK^2)\).

Sample code (C++)