A similar but easier question is
“What is the probability that when a dice is rolled $3$ times, the sum of score is $13$ ?”
For which you can easily solved by listing down all the possible combination of three numbers $(\leq 6)$ to get $13$.
But if the dice is rolled a lot of times, it becomes very tedious to list down all the combinations. Is there easier method to solve this type of question?
Any help is much appreciated.
On
One way to go:
So we're looking for the number of $(a_1, a_2, \dots, a_6)$ for which $a_1 + a_2 + \dots + a_6 = 20$, with $1 \le a_i \le 6$ [we just divide that by $6^6$].
If the constraint $a_i \le 6$ wasn't there we could use stars and bars https://en.wikipedia.org/wiki/Stars_and_bars_(combinatorics). Define $b_i = a_i-1$. It's equivalent to count the number of $(b_1, b_2, \dots, b_6)$ with $b_1 + \dots + b_6 = 20-6 = 14$ and $0 \le b_i \le 5$. Without the $b_i \le 5$ constraint, we'd have $\binom{14 +6-1}{6-1} = 11628$ total solutions from stars and bars.
Let's subtract those solutions out of these which violate some $b_i \le 5$ constraint. Let "constraint 1" denote $b_1 \le 5$, "constraint 2" $b_2 \le 5$, etc.
Note that at most 2 constraints may be violated (in solutions among those 11628). This should simplify things.
We can use an inclusion-exclusion approach. The number of solutions which violate constraint 1 are those with $b_1 \ge 6$. Letting $c_1 = b_1-6$, this is the number of $(c_1, b_2, \dots, b_6)$ with $c_1 + b_2 + \dots b_6 = 14-6 = 8$; by stars and bars this is $\binom{8+6-1}{6-1} = 1287$ solutions.
Similarly, there are 1287 solutions that violate constraint $i$ for $2 \le i \le 6$ as well.
Now let's count the number of solutions which violate some pair of constraints. It's same for any pair; let's just count for constraint 1 and constraint 2. That is we want solutions with $b_1 \ge 6, b_2 \ge 6$. Setting $c_1 = b_1 - 6, c_2 = b_2-6$, this the number of solutions to $c_1 + c_2 + b_3 + \dots + b_6 = 14-2(6) = 2$. This is $\binom{2+6-1}{6-1} = 21$ solutions.
And we're ready to do the inclusion-exclusion. Subtract the number which violate one constraint, and add back in the number which violate any pair b/c those were double-subtracted.
$$11628 - 6(1287) + \binom{6}{2}(21) = 4221$$
You can solve such problems using generating functions, which does the counting of all the possibilities using a formal weight factor to keep track of the score. You can then sum over all possible dice throws without any restrictions as the weight factor does the job of keeping track of the score.
If we take the weight for a score of $k$ for a single throw to be $x^k$, then the sum of all possible outcomes is
$$w(x) = \sum_{k=1}^6 x^k = x \frac{1-x^6}{1-x}$$
The coefficient of $x^r$ of the series expansion of $w(x)$ is then the number of ways you can get a score of $r$, which is $1$ for $1\leq r\leq 6$ as we're only throwing a single die one time. If we throw 3 times, then the generating function for that will be $w(x)^3$ (it's a good exercise to try to understand why this is true).
In general, if we throw $n$ times we need to expand the function:
$$w(x)^n = x^n \left(\frac{1-x^6}{1-x}\right)^n \tag{1}$$
If $n$ is not too large, one can directly expand the function. E.g. in case of $n = 6$ the series expansion of the numerator in (1) is:
$$\frac{1}{(1-x)^6} = \sum_{r=0}^{\infty}\binom{r+5}{5}x^{r}$$
and multiplying that by $x^6 \left(1-x^6\right)^6 = x^6 - \binom{6}{1}x^{12} + \binom{6}{2}x^{18} - \cdots$ to extract the coefficient of $x^{20}$ isn't a lot of work, the result is:
$$\binom{19}{5} - \binom{6}{1} \binom{13}{5} + \binom{6}{2} \binom{7}{5} = 4221$$
For much larger $n$ we can approximately evaluate the series expansion coefficient as follows. It follows from Cauchy's integral formula that the coefficient $c_k$ of $z^k$ of an analytic function $f(z)$ is given by:
$$c_k = \frac{1}{2\pi i}\oint_{C}\frac{f(z)}{z^{k+1}}dz\tag{2}$$
where $C$ is a counterclockwise contour that encircles the origin. We can then estimate $c_k$ by evaluating the contour integral using the saddle point method. We find saddle points by setting the logarithmic derivative of the integrand equal to zero. If we deform the contour integral such that it moves through these points, then it will pick up its dominant contributions from these saddle points, and there will usually be one that will make the most important contribution. The contributions from each saddle point can be approximated by a Gaussian integral.
While this will yield good results for large $n$, we can see that it yields reasonably good approximations for not too small $n$. If we consider again the case considered above with $n = 6$ and a score of $20$, we find that the saddle point is at $z = 1$, the integrand (2) for $z = 1 + u$ can be approximated as:
$$\frac{1}{2\pi i} 6^6 \oint_{C}\exp\left(\frac{35}{4}u^2\right) du $$
This will then be integrated in the imaginary direction, and this then yields the Gaussian approximation of:
$$\frac{6^6}{2\pi}\sqrt{\frac{4\pi}{35}} = \frac{6^6}{\sqrt{35\pi}}\approx 4449$$