Consider a two-variable generating function \(f(x,y)\) where \([x^Ny^M]f(x,y) =\) (the number of wonderful integer sequences of length \(M\) with sum \(N\)).

We apply the inclusion-exclusion theorem. We have \((-1)^{j-1} x^{ij} y^{j}\) for sequences containing \(i\) consecutive copies of \(j\)’s. We have \((-1)^{j-1}\) because it violates the condition at \((j-1)\) positions. Let \(g(x,y)\) be the sum of this terms for all \(i,j \ge 1\); then \(f(x,y) = \frac{1}{1-g(x,y)}\).

We have \(g(x,y) = \sum_{i=1}^{\infty} \sum_{j=1}^{\infty} (-1)^{j-1}x^{ij}y^{j} = \sum_{i=1}^{\infty} \frac{x^iy}{1+x^iy}\), but naively finding \(f(x,y)\) from this is still difficult.

Now, consider a generating function \(x\) such that \([x^N] h(x)=\) (the sum of lengths of wonderful integer sequences with sum \(N\)); then \(h(x)\) equals \(\frac{d}{dy} f(x,y)\) where \(y=1\).

We have \(\frac{d}{dy} f(x,y)=\frac{d}{dy} \frac{1}{1-g(x,y)} = \frac{\frac{d}{dy} (g(x,y) - 1)}{(1-g(x,y))^2}\).

Consider \(\frac{d}{dy} (g(x,y)-1)\). Since \(\frac{d}{dy} \frac{x^iy}{1+x^iy} = \frac{x^i}{(1+x^iy)^2}\), we have \(\sum_{i=1}^{\infty} \frac{x^{i}}{(1+x^iy)^2}\).

Assigning \(y=1\), we obtain \(h(x) = \displaystyle \frac{\sum_{i=1}^{\infty} \frac{x^{i}}{(1+x^i)^2}}{\left(1-\sum_{i=1}^{\infty} \frac{x^i}{1+x^i}\right)^2}\).

Expanding the terms in the numerator and denominator is fast enough because there are only \(\mathrm{O}(N\log N)\) terms with degrees at most \(N\). By finding the inverse of the polynomial in an \(\mathrm{O}(N\log N)\) time, the problem can be solved in a total of \(\mathrm{O}(N\log N)\) time, which is fast enough.

The way of finding the inverse of a polynomial is explained in the editorial of ABC260-Ex.