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)));
