It is known how to compute the Fourier expansion of elliptic modular forms using modular symbols, and it is known how to get numerical evaluations of $L$-functions of various type ; it's possible to get explicit values on those matters with sage already.
It is also known that when an Eisenstein series is involved, it's possible to relate the Petersson scalar product to $L$-functions, and hence to evaluate them.
I have seen that sage bug about various pairings for modular forms, but it looks more like it's about the pairing between modular forms and modular symbols than the Petersson scalar product.
My question is: does there exist general formulas to compute the Petersson scalar product of two elliptic modular forms numerically?
EDIT(2012-12-23): I insist on the numerically: having an expansion with estimates on the order of the error with constants which depends on this or that (I'm thinking about those which can be found in chapter 5 of Iwaniek's "Topics in classical automorphic forms" for example) is very nice from a theoretical point of view, but doesn't help when one wants to actually compute with specific forms and to a given precision. In fact, I want to compute various things with the Petersson scalar product, so this question is to check whether I can directly work on them or if I should write something about the matter before.
Best Answer
There is a "quick and dirty" way to find the inner product of two cusp forms that are not necessarily Hecke eigenforms. I learned this from Akshay Venkatesh.
The formula is that \begin{equation*}\langle f, g \rangle = \lim_{y \rightarrow 0^+} y^k \int_0^{1} f(x+iy) \overline{g(x+iy)} dx, \end{equation*} where $f$ and $g$ are weight $k$ and the inner product is normalized via \begin{equation*} \langle f, g \rangle = \int_{\Gamma \backslash \mathbb{H}} y^k f(z) \overline{g(z) }\frac{1}{V} \frac{dx dy}{ y^2}, \end{equation*} where $V$ is the volume of $\Gamma \backslash \mathbb{H}$. The philosophy behind the proof is that the horocycle $x+iy: 0 \leq x \leq 1$ equidistributes in the fundamental domain as $y \rightarrow 0$. You can prove the formula by spectrally decomposition $f \overline{g}$. The projection onto the constant eigenfunction gives $\langle f, g\rangle$. The projections onto the cusp forms integrate out to zero. The projection onto the Eisenstein series leaves the constant terms which are bounded by $\sqrt{y}$, and hence have limit zero as $y$ tends to $0$.
If $f(z) =\sum_n a(n) e(nz)$ and $g(z) = \sum_n b(n) e(nz)$, then of course \begin{equation*} \int_0^{1} f(x+iy) \overline{g(x+iy)} dx = \sum_{n \geq 1} a(n) \overline{b(n)} \exp(-4 \pi n y). \end{equation*}