function [pdf,cdf,pdfo,cdfo] = poisson(mu,K);
% function [pdf,cdf,pdfo,cdfo] = poisson(mu,K);
% By comparing pdf0 to pdf you can see the increased accuracy
% that follows from computing in the "exp(log)" form.

% original method
pdfo = exp(-mu) *cumprod( [1 mu*ones(1,K)] ./ [1 1:K] );
cdfo = cumsum(pdfo);

% modified method using logs
pdf= exp(-mu + cumsum( [0 log(mu)+zeros(1,K)] - log([1 1:K]) ));
cdf = cumsum(pdf);