Generating Function for sum of N dice [or other multinomial distribution] where lowest N values are “dropped” or removed

Question

Background

I found this interesting question Formula for dropping dice (non-brute force) and excellent answer https://stats.stackexchange.com/a/242857/221422, but couldn't figure out how to generalize a generating function for when more than one die is dropped. Similarly, I'm having difficulty figuring out a related mechanic for when the highest roll is dropped.

Description of the Problem

Suppose you have $N$ fair dice each with $S$ sides. Roll all the dice and then remove the lowest [or highest, alternatively] $M$ (where $M > 0$ and $M < N$) dice and then sum the remainder. What is the probability distribution of the sum? Specifically, how does one go about finding the generating polynomial?

Implementation of whuber's answer

I found whuber's answer to be incredibly thorough. I thought it might be nice to see how to actually implement it in code, so I've pasted that below.

from functools import reduce

from numpy.polynomial import polynomial as p


def generating_function(k, d, n):
    return p.polypow(
        [0] * k + [1] * (d - k + 1),
        n
    )


def drop_one_die(n, d):
    tmp = [
        generating_function(k, d, n) for k in range(1, d + 2)
    ]

    differences = (
        (tmp[i] - tmp[i + 1])[i + 1:] for i in range(d)
    )

    return reduce(p.polyadd, differences)


print(
    drop_one_die(4, 6)
)

Other considerations / Multinomial distribution

To generalize even further, instead of a fair die where each outcome is equally likely, what if you start with a general multinomial distribution?

So instead of

$$(1/6)x + (1/6)x^2 + (1/6)x^3 + (1/6)x^4 + (1/6)x^5 + (1/6)x^6$$

you start with

$$p_0 + {p_1}{x} + {p_2}{x^2} + ... + {p_n}{x^n}$$

Thanks!