It's claimed on Page 54 of Bakalov and Kirillov's Lectures on Tensor Categories and Modular Functors that the formula in Scott's post above is first proven in Moore and Seiberg's Classical and Quantum Conformal Field Theory). I've only just grabbed that paper, but equation (A.7) gives a generalization of the formula above to "n-point function characters at genus g" as
$$
dim\, V(g,i_1,\ldots,i_n)=\sum_p\frac{S_{i_1 p}\ldots S_{i_n p}}{S_{0p}\ldots S_{0p}}\left(\frac{1}{S_{0p}}\right)^{2g-2}
$$
which for the case $n=3$,$g=0$ looks like $\mathcal D^2$ times Scott's formula since $N_{ij}^k=dim \, (V_{ij}^k)$ is defined on page 8.
Separately, in section (2.8.3) of DGNO's On Braided Fusion Categories I they provide what they call a Verlinde formula for a premodular category $\mathcal C$
$$
S_{XY}S_{XZ}=d(X)\sum_{W\in\mathcal O(\mathcal C)}N_{YZ}^W S_{XW}, \,\,\,\,\,\,\,\,\,\,\,\, X,Y,Z\in \mathcal O(\mathcal C)
$$
equivalent to the one in Scott and Marcel's answers when the S matrix is invertible. Their references for this are 1 as theorem 3.1.12 (which I think should be proposition 3.1.12 in the online edition and prop 3.1.13 in the printed edition, given that the formula is obtained in (3.1.31) as rewriting (3.1.26)) and Muger's On the structure of Modular Categories as Lemma 2.4, though it look like Marcel's Reference from Rehren may be the first to prove this in the pre-modular setting.
I think this is worth mentioning in its own right because later in 3 section 3.4.2 DGNO prove a "non-spherical analogue of the verlinde formula" as stated in the previous equation, via
$$
\tilde S_{XY}\tilde S_{XZ}=d(X)\sum_{W\in\mathcal O(\mathcal C)}N_{YZ}^W \tilde S_{XW}, \,\,\,\,\,\,\,\,\,\,\,\, X,Y,Z\in \mathcal O(\mathcal C)
$$
where $\mathcal C$ is a braided
$$
\tilde S_{XY}=\frac{Tr_{-}\otimes Tr_{+}(c_{Y,X}c_{X,Y})}{d_{-}(X)d_{+}(Y)},
$$
$Tr_{-}/Tr_{+}$ and $d_{-}/d_{+}$ are the left/right trace and dimension defined as in their section 2.4.2. This would then be yet another version in a slightly more general context. In the context of a premodular category we have that $S_{XY}=d(X)d(Y)\tilde S_{XY}$.