Numerical evaluation of the Petersson product of elliptic modular forms 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.
 A: 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*}
A: Let
$f(z) = \sum a(n) e(n z)$
and
$g(z) = \sum b(n) e(n z)$
be holomorphic modular forms
of weight $k \in 2 \mathbb{N}$ on $\Gamma =
\operatorname{SL}_2(\mathbb{Z})$
whose product decays rapidly.
Then
\begin{equation}
\int_{\Gamma \backslash \mathbb{H}} y^k \overline{f(z)}g(z) ~ \frac{dx ~ d y}{y^2} =
2
\sum_{n \in \mathbb{N} } \frac{  \overline{a(n)} b(n) }{  n^{k-1} } \sum_{d \in \mathbb{N}} \Phi(4 \pi d \sqrt{n}),
\end{equation}
where
\begin{equation}
\Phi(y)
= 2( \frac{y}{8 \pi})^{k-1}
(y K_{k-2}(y) - K_{k-1}(y)).
\end{equation}
Note that $\Phi(y)
\asymp_k
y^{k-1/2}
e^{-y}$ for $y \gg 1$.
For a general finite index subgroup
 $\Gamma$
 of $\operatorname{SL}_2(\mathbb{Z})$,
a correct formula
may be obtained
by
summing the RHS over the cusps
$\mathfrak{a}$ of $\Gamma$,
weighted by the width $w$ of $\mathfrak{a}$,
and taking for $a(n), b(n)$ ($n \in w^{-1} \mathbb{N}$) the Fourier coefficients at
$\mathfrak{a}$.
A reference
for
a general form of such identities
is Theorem 5.6 (p.24)
in my paper Evaluating modular forms on Shimura curves;
see also Example 5.7,
Remark 3.5,
and the discussion of Sections 5.3--5.6,
which includes a detailed comparison with the other approaches
mentioned in this thread that I will summarize briefly here.
Pros: no need to Hecke-decompose, unlike the "symmetric square"
approach; converges rapidly to the correct value,
unlike the vanilla "equidistribution of horocycles" approach;
generalizes to non-holomorphic forms
lacking a straightforward Hecke decomposition (e.g.,
certain theta series),
although perhaps this feature
is not important for your purposes.
Cons: requires the Fourier expansion at every cusp,
unlike either approach just mentioned
(although the "symmetric square" approach
is not devoid of such subtlety,
since it requires one to compute
the conductor and bad Euler factors
of the symmetric square of a newform).
Another approach (specific to the holomorphic case)
would be to exploit the connection with period polynomials,
for which
a search just now turned up
this article.
One variant of that method also requires knowing Fourier expansions at every cusp, and reduces the problem to evaluating a class of incomplete gamma functions some of which reduce to K-Bessel functions as above; another requires only that one be able to compute periods of a cusp form $f$ over split geodesics in $\Gamma \backslash \mathbb{H}$, which can apparently be done using the Fourier expansion at only one cusp.  Moreover, one can speed up the computation when $f = g$ is an eigenform by exploiting certain rationality results.
A: It's easy to reduce to the case of computing the Petersson product of a normalised new eigenform with itself. Here you can use the fact that the product is equal to the value at s=k of the symmetric square L-function of f, and this you can compute using e.g. Tim Dokchitser's algorithms. Here is a thread from the Sage developers mailing list with example code by Martin Raum: https://groups.google.com/forum/m/#!topic/sage-nt/EkBWOogY8yw
For elliptic curves there is also Mark Watkins' Sympow program, which will compute all the symmetric power L-functions.
