I write here the function I had in mind, using the solution I found:

```
def br(x):
return sqrt(x)-sqrt(1/x)
import sympy as sp
k=3
q1,q2 = var('q1,q2')
X = sp.symbols('x_0:{}'.format(k))
rho = [1,q1,q2]
chi = 1/prod([ br(X[j])*br(q1/X[j])*X[j]*br(q2/X[j]) for j in range(k)])
for i in range(k):
chi = sp.residue(chi,X[i],rho[i])
chi
```

The problem is that even for a stupid example such as the one above, it takes forever to compute it.

Here's the same (with the actual function I need), using the above suggested answer:

```
k=3
br(x)=sqrt(x)-sqrt(1/x)
q1,q2,q3,m=var('q1,q2,q3,m')
# X=var('x_0,x_1,x_2')
# X=var('x_0:{}'.format(k))
X = [var("x%d" % i) for i in range(k)]
rho=[1,q1,q2]
# rho=[1,q1]
chi1 = prod([ br(X[j]/m)/(br(X[j])*X[j]) for j in range(k)])
chi2 = prod([prod([br(q1*q2*X[i]/X[j])*br(q1*q3*X[i]/X[j])*br(q2*q3*X[i]/X[j]) for i in range(k) if i != j]) for j in range(k)])
chi3 = prod([prod([br(q1*X[i]/X[j])*br(q2*X[i]/X[j])*br(q3*X[i]/X[j])*br(q1*q2*q3*X[i]/X[j]) for i in range(k) if i != j]) for j in range (k)])
chi = chi1*chi2/chi3
for xi,rhoi in zip(X,rho):
chi = chi.residue(xi==rhoi)
print chi.factor()
```

This gives zero (wrong) for k=3 (unless I made a stupid mistake) and it is still slow, but faster than sympy for k=2.