I would suggest using a generating polynomial to do this. The idea is to represent the contribution of one die towards the total by the power series $(1 + x + \dots + x^n)$$(x^1 + x^2 + \dots + x^n)$ where $n$ is the number of faces on the die, then multiplying together one series per die and calculating the coefficient of the power of the sum you're looking for, e.g., the coefficient of $x^{41}$ in your example. Then you can divide by the number of different possible die rolls to get the result as a probability.
This works because expanding all the terms in the product is analogous to listing all the possible combinations of die rolls on your dice, e.g., rolling a 1 on all 6 dice in the question is represented by choosing the $x^1$ term from the series for each die and multiplying these together to get $x^6$. So with the 6 dice in the question, there is only 1 way to roll a total of 6. Similarly, there will be many different rolls that give a sum of 41 and number of these will be equal to the coefficient of $x^{41}$.
You can reduce the degree of combinatorial explosion by discarding terms where the power of $x$ exceeds the target, or cannot reach it based on the remaining terms in the product.
Wolfram Alpha pseudocode:
Coefficient(sum(x^n, n=1 to 10)^2 * sum(x^n, n=1 to 8)^4, x^41) / (10^2 * 8^4)
