There are already a bunch of great answers: [1] gives formula for $f(p)$ in terms of infinite series and, moreover, explains its oscillating behavior.
I intend to elaborate on Markov chain approach used by [2], but in matrix form instead of recurrence relations, showing that it is possible to reproduce formula for $f(p)$ in terms of infinite series and in terms of finite series of rational functions.
Let's start from a coin toss game with $N$ coins (players) and only one absorbing state with zero coins.
Let $p$ be the probability of getting heads (or player staying after one round).
Transition probabilities are described by binomial distribution
$$T_{i j} = P(i\,\text{coins} \rightarrow j\,\text{coins}) = \binom{i}{j} p^j (1 - p)^{i - j},\quad i = 0, 1, \ldots N,\quad j = 0, 1, \ldots i$$
Transition matrix of this game is multiplicative in a sense that:
$$T(p_1) T(p_2) = T(p_1 p_2)$$
Proof is rather direct.
Notice that $T$ is lower-triangular block-diagonal matrix
$$T(p) = \begin{pmatrix} A & 0 \\ C & D \end{pmatrix}$$
where $A$ and $D$ are arbitrary square blocks.
Diagonal blocks do not depend on other blocks if matrices are multiplied.
$$\begin{pmatrix} A & 0 \\ C & D \end{pmatrix} \begin{pmatrix} A^\prime & 0 \\ C^\prime & D^\prime \end{pmatrix} = \begin{pmatrix} A A^\prime & 0 \\ C A^\prime + D C^\prime & D D^\prime \end{pmatrix}$$
Now consider coin toss game from the question, in which state with one coin is also absorbing.
Its transition matrix is almost the same, except that it has 2x2 identity matrix as top left block, reflecting that states with zero and one coins are absorbing states:
$$T^\prime(p) = \begin{pmatrix} I_2 & 0 \\ C & D \end{pmatrix}$$
Let $\pi_0 = \begin{pmatrix} 0 & 0 & \ldots & 1 \end{pmatrix}$ be vector of initial probability distribution (state with $N$ coins).
To obtain it's evolution after $l$ steps, we have to multiply it by power of transition matrix ($D^l$ vanishes because, as truncated transition matrix, its eigenvalues are all less than $1$ by absolute value):
$$\pi_l = \pi_0 T^{\prime \, l}$$
$$
T^{\prime \, l} = \begin{pmatrix} I_2 & 0 \\ \sum_{k=0}^{l-1} D^k C & D^l \end{pmatrix} \xrightarrow{l\rightarrow\infty} T^{\prime \, \infty} =
\begin{pmatrix} I_2 & 0 \\ \sum_{k=0}^{\infty} D^k C & 0 \end{pmatrix}
$$
So the $f(p)$ in question is given by the $(T^{\prime\,\infty})_{N, 1}$. But because blocks $D$, $C$ are the same as in original transition matrix, we can exploit multiplicativity of $T$: entries of $D^k C$ differ from corresponding ones in block inside $T^k T = T^{k + 1}$ by a term, namely
$$
T^k T = \begin{pmatrix} A^{k} A & 0 \\ \underbrace{\left(\sum_{m=0}^{k-1} D^m C A^{k - 1 - m}\right) A}_{\text{extra term}} + D^k C & D^{k} D \end{pmatrix}
$$
\begin{align}
(D^k C)_{N, 1} &=
\sum_{r = 2}^{N} (D^k)_{N, r} C_{r, 1} =
\sum_{r = 1}^{N} (T^k)_{N, r} T_{r, 1} - \underbrace{(T^k)_{N, 1} T_{1, 1}}_\text{from extra term} =\\
&= (T^{k + 1})_{N, 1} - (T^k)_{N, 1} T_{1, 1} = \binom{N}{1} p^{k + 1} (1 - p^{k + 1})^{N - 1} - \binom{N}{1} p^{k} (1 - p^{k})^{N - 1} p =\\
&= N p^{k + 1} \left((1 - p^{k + 1})^{N - 1} - (1 - p^{k})^{N - 1} \right)
\end{align}
Now we obtain $f(p)$ in terms of infinite series (apparently it converges thanks to $p^k$ factor) after some rearrangements:
\begin{align}
f(p) &= \sum_{k=0}^{\infty} N p^{k + 1} \left((1 - p^{k + 1})^{N - 1} - (1 - p^{k})^{N - 1} \right) =\\
% &= \sum_{k=0}^{\infty} N p^{k + 1} (1 - p^{k + 1})^{N - 1} - \sum_{k=1}^{\infty} N p^{k + 1} (1 - p^{k})^{N - 1} =\\
&= \sum_{k=1}^{\infty} N p^{k } (1 - p^{k })^{N - 1} - \sum_{k=1}^{\infty} N p^{k + 1} (1 - p^{k})^{N - 1} =\\
&= \sum_{k=1}^{\infty} N p^{k} (1 - p) (1 - p^{k})^{N - 1}
\end{align}
which passes cross-check with [1].
Let's expand power and take an infinite sum:
\begin{align}
f(p) &= \sum_{k=1}^{\infty} N p^{k} (1 - p) \sum_{r=0}^{N-1} \binom{N-1}{r} (-1)^r p^{r k} =\\
&= \sum_{r=0}^{N-1} N \binom{N-1}{r} (-1)^r \frac{1 - p}{1 - p^{r + 1}}
\end{align}
This finite sum is related to recurrence relations resulting in large, impractical expressions mentioned in [2]. Indeed, reduction of the fractions in the finite sum to a common denominator would yield a huge rational function.
We can obtain another series representation after rearranging terms in the infinite sum:
$$
f(p) = \sum_{k=1}^{\infty} N (1 - p) \sum_{r=0}^{N-1} a_{r + 1} p^{k (r + 1)}
$$
where $ a_{r + 1} = \binom{N-1}{r} (-1)^r $
After expanding both sums we get
\begin{alignat*}{5}
a_1 p^1 &{}+{}& a_1 p^2 &{}+{}& a_1 p^3 &{}+{}& a_1 p^4 &{}+{}& \ldots &{}+ \\
&{}+{}& a_2 p^2 & & &{}+{}& a_2 p^4 &{}+{}& \ldots &{}+ \\
& & &{}+{}& a_3 p^3 & & &{}+{}& \ldots & \\
\end{alignat*}
and finally
$$
f(p) = N (1 - p) \sum_{k=1}^{\infty} p^k \sum_{\substack{d \vert k \\ d \leq N}} a_d
$$
This one looks interesting because it features number theory-like summing over divisors.
Perhaps playing around with this expression can provide further insights, but I've reached limits of what I know at this point.
