code: def _psi_n(x, n, b):
      
      return 2**(b-1) / gamma(b) * (-1)**n * \
      np.exp(gammaln(n+b) -
             gammaln(n+1) +
             np.log(2*n+b) -
             0.5 * np.log(2*np.pi*x**3) -
             (2*n+b)**2 / (8.*x))
docs: Compute the n-th term in the infinite sum of
      the Jacobi density.