%PARFRAC
%  The matrix c of order (m,2*m) of coefficients of the partial
%  fraction decomposition (3.8)
%
function c=parfrac(dig,m,z);
digits(dig);
y=imag(z);
c=vpa(zeros(m,2*m));
if m==1
  c(1,1)=1/(2*i*y); c(1,2)=-c(1,1);
else
  for mu=1:m
    c(mu,mu)=1/(2*i*y)^m;
    if mu>1
      for nu=mu-1:-1:1
        c(mu,nu)=-(1+(m-1)/(mu-nu))*c(mu,nu+1)/(2*i*y);
      end
    end
    c(mu,mu+m)=(-1/(2*i*y))^mu;
    for kappa=m-1:-1:1
      c(mu,mu+kappa)=(1+(mu-1)/(m-kappa))*c(mu,mu+kappa+1)/(2*i*y);
    end
  end
end
