function [psi,psiprime]=digamma(a) % outputs psi(a) and psiprime(a) where psi is the digamma function k=100000; n=1:k; gamma=.5772156649; psi=-gamma+(a-1)*(1/(k+.5+a)+sum(1./(n.*(n+a-1)))); psiprime=(1/(k-.5+a))+sum(1./((n-1+a).*(n-1+a)));