We will get the bivariate generating function $$$F(z, t) = \sum_{n\ge 1}\sum_{k=1}^n A_{n,k}\frac{z^nt^k}{n!}$$$, where $$$A_{n,k}$$$ is the sum of occurrences of $$$k$$$ in all good sequence of length $$$n$$$. Consider the inclusion-exclution principle. For all contribution with $$$\max k = m$$$, we have

This means that for all $$$i \in S$$$, we force all recurrences of $$$i$$$ are before $$$i + 1$$$, and the occurrences of $$$j$$$ are counted.

After some tough calculation, we will found out that this equates to $$$\frac{t(\mathrm e^{z(1-t)}-1)}{(1-z) (1-t \mathrm e^{z(1-t)})}$$$. I will post the details later.

Now our goal is to calculate

We consider $$$\left([z^n]\frac{t(\mathrm e^{z(1-t)}-1)}{(1-z) (1-t \mathrm e^{z(1-t)})}\right) + 1 = (1-t)[z^n] \frac 1{(1-z)(1-t\mathrm e^{z(1-t)})}$$$, which looks simpler. We let $$$z = \frac u{1-t}$$$, then

Hence we have

Noticed that $$$[u^m]\mathrm{e}^{ku} = \frac{k^m}{m!}$$$, this could be calculated straightly through multipoint evaluation with time complexity of $$$\Theta(n\log^2 n)$$$. And $$$[u^m] \frac1{(1-u)^k} = \binom{n+k-1}{n} = \frac{(n+k-1)!}{n!(k-1)!}$$$ so this part could be calculated through a convolution. It will pass if your implementation doesn't have a big constant.

It could also be improved to $$$\Theta(n\log n)$$$ through the Lagrange Inversion formula similar to the original solution, I leave this as an exercise.

UPD1: simplified some deduction.