You can take a generating function approach here. That is, you let $x$ be an indeterminate and find the $40!$ times the coefficient of $x^{40}$ in $$\left(\frac{x^8}{8!} + \frac{x^9}{9!} + \cdots + \frac{x^{16}}{16!}\right)^4.$$
The reason this works is that when you expand the product in all possible ways, factoring out a factorial of the degree, each product has a multinomial coefficient. For example, taking $$\frac{x^{11}}{11!} \frac{x^9}{9!} \frac{x^8}{8!} \frac{x^{12}}{12!} = \frac{40!}{11!\, 9!\, 8!\, 12!}\frac{x^{40}}{40!} $$ and multiplying by $40!$ gives the coefficient $$\frac{40!}{11!\, 9!\, 8!\, 12!}$$ which is the number of ways of choosing groups where the first supervisor has $8$ students, the second supervisor has $9$ students, etc. You only have to go up to degree $16$ in each factor because no group can have more than $16$ people.
Of course, expanding this out completely is equivalent to adding up all the possible multinomials corresponding to how many students each supervisor takes as JMoravitz mentions. But software packages like Sage are optimized to do this kind of convolution pretty quickly.
Edit: In Sage,
p = sum([x^n / factorial(n) for n in range(8, 17)])
q = (p^4).expand()
print q.coefficient(x,40) * factorial(40)
prints $435451605680654896510320$.