Let \(E_k\) be the expected number of occurrences of \(k\) after all operations.

Fix \(k \geq 2\) and consider the event “\(k\) occurs for the first time by the \(i\)-th operation.” Then by the \((i-1)\)-th operation, operation \(1\) is chosen exactly \((k-2)\) times, and in the \(i\)-th operation, operation \(1\) is chosen. The probability that this happens is \(\binom{i-1}{k-2}p^{k-1}(1-p)^{i-k+1}\). Now assume that this event has happened.

Let \(f(N-i-1)\) be the expected number of \(k\) chosen during the remaining \((N-i-1)\) operations (which is independent of \(k\)). If \(A\) has \(n\) elements, \(c\) elements of which equal \(k\), then the expected number of \(k\) after one operation is \(c+(1-p)\frac{c}{n}=\frac{n+1-p}{n}c\). Thus, \(f(N-i-1)=\frac{(N-p)\cdots(i+2-p)}{(N-1)\cdots(i+1)}\). By the recurrence relation \(f(0)=1, f(x)=f(x-1)\frac{N-x+1-p}{N-x}\), with appropriate precalculation we can compute \(f(0),f(1),\dots,f(N)\) in \(O(N)\) time.

Using this, \(E_k\) can be computed as follows:

\[ E_k = \sum_{i=k-1}^{N-1} \binom{i-1}{k-2}p^{k-1}(1-p)^{i-k+1}f(N-i-1) \\ = \sum_{i=0}^{N-k} \binom{i+k-2}{k-2}p^{k-1}(1-p)^{i}f(N-k-i)\\ =\frac{p^{k-1}}{(k-2)!} \sum_{i=0}^{N-k} \frac{(1-p)^i}{i!}\cdot f(N-k-i)(i+k-2)!\\ =\frac{p^{k-1}}{(k-2)!} \sum_{i=0}^{N-k} \frac{(1-p)^i}{i!}\cdot f(N-k-i)(N-(N-k-i)-2)! \]

Therefore, if we define polynomials \(F(x)\) and \(G(x)\) as \(F(x)=\sum_{j=0}^{N-2} f(j)(N-j-2)! x^j\) and \(G(x)=\sum_{j=0}^{N-2} \frac{(1-p)^j}{j!}x^j\), the answer for \(k\) is represented as \(\frac{p^{k-1}}{(k-2)!}[x^{N-k}]F(x)G(x)\). The polynomials \(F(x)\) and \(G(x)\) can be constructed in \(O(N)\) time, and their product can be found in \(O(N\log N)\) time with Fast Fourier Transform. If you can use AtCoder Library, atcoder::convolution is a convenient function to find the product.

\(E_1\) can be found as \(E_1 = N - \sum_{k=2}^N E_k\).