Sorry if the thoughts are a bit disjoint
I think I have a solution which is faster than the one I propose above. Essentially, I view each die as a finite dimension, in a geometry where only movements parallel to the directions of the dimensions are valid. Then, each value of the dice is a constant distance from the origin, and the sum of those and lesser values can be obtained through the integral under the curve of all points of that constant distance (back to the origin).
The advantage of using such a geometry in N-space is largely that it makes visualization of what components are useful from one calculation to the next....
Wait - I'm being stupid, aren't I - these are just the N-spacial icosahedral numbers, aren't they? In 2-space, with evenly sized dice, they'd be:
X*(X+1)/2
In 3-space, they'd be:
X*(X+1)*(X+2)/6
So in general, they're:
Where A = the number of dice,
(X*(X+1)*(X+2)*...*(X+A-3)*(X+A-2)*(X+A-1))/(A*(A-1)*(A-2)*...3*2*1)
This obviously only works directly to the mid-way point, but the same series can be used to calculate beyond this point- it merely needs to be used multiple times - once for a superset of the N-space we want to use, then once for each portion of that set we aren't planning to use. Since the maxima is inside the set we're planning to use, we don't need to worry about the overlap beyond that point.
I think that this can even be used for non-evenly sized dice, but I'll have to work at that more. In any case, a direct equation around the order of the square of the number of dice is far better than I expected to be able to obtain.