How to compute poles and values of L-functions? Let $L(s)$ be an L-function, given by its series expansion and admitting an Euler product, say, for $s$ of large enough real part,
$$L(s) = \sum_n \frac{a_n}{n^s} = \prod_p \prod_{i=1}^k (1-\alpha_i p^{-s})^{-1}.$$
Assume that $L(s)$ admits a simple pole at $s=x$ with residue $R$, so that it admits a local expansion of the form
$$L(s) = \frac{R}{s-x} + H(s),$$
where $H$ is an analytic function around $s=x$. My question is therefore: 

Assuming we now explicitly the value of the residue $R$, is there any standard
  way of computing the value of $H(x)$?

(everything here is quite explicit, but I am still stuck when I try to do something with the cumbersome expression $L(s) - R/(s-x)$...)
NB: If it can help, I am interested in some kind of symmetric square L-functions.
 A: The computation of the Taylor expansion of an $L$-function is implemented in Pari/GP. Here is an example for the Riemann zeta function at $s=1$:
Z=lfuncreate(1);
lfun(Z,1+x+O(x^4))

which returns
1.0000000000000000000000000000000000000*x^-1 + 0.57721566490153286060651209008240243104 + 0.072815845483676724860586375874901319140*x + O(x^2)

For the theory behind, you can look at Henri Cohen's course. The documentation for Pari/GP is available here.
A: I don't claim my method is very clever, but let's look at this for the Riemann $\zeta$-function.
A standard trick is integrate around a symmetric rectangle, infinite in the vertical direction.
For instance with $\xi(s)=\Gamma(s/2)\zeta(s)/\pi^{s/2}$ consider
$$\biggl(\int_{(2)}-\int_{(-1)}\biggr)(s-1/2)\xi(s){1\over s(s-1)}{ds\over 2\pi i}.$$
Here I've made there be a double pole at $s=0,1$, so that the desired term should fall out of a residue calculation.
On the one hand, Cauchy's Residue Theorem says this is twice (by symmetry) the residue at $s=1$ of
$$((s-1)+1/2))({1\over s-1}+c_0+\cdots){1\over s-1}\sum_{n=0}^\infty (1-s)^n$$
and if my math is correct, this residue is $1+c_0/2-1/2$.
On the other hand, by symmetry this difference of integrals is twice the integral on the 2-line, where there is a Mellin expansion:
$$\sum_{n=1}^\infty\int_{(2)}{\Gamma(s/2)(s-1/2)\over s(s-1)}{1\over\pi^{s/2}}{1\over n^s}{ds\over 2\pi i},$$
and this latter integral can be computed if desired to high precision for each $n$ (and by general theory the sum over $n$ converges). Perhaps one can be more clever, and actually get the result for $\zeta(s)$ theoretically, being able to compute the transform explicitly, and sum over $n$, but I don't know how offhand.
In general there will be coefficients depending on $n$, and likely one can't do more than just compute to some precision.
