function [d,pd] = mgdf(pn,y,py)
% MGDF [d,pd] = mgdf(pn,y,py)  Function version of mgD
% Version of 5/19/97
% Uses m-functions mgsum and csort
% Do not forget zeros coefficients for missing
% powers in the generating function for N 
n  = length(pn);               % Initialization
a  = 0;
b  = 1;
d  = a;
pd = pn(1);
for i = 2:n
  [a,b] = mgsum(y,a,py,b);
  d  = [d a];
  pd = [pd b*pn(i)];
  [d,pd] = csort(d,pd);
end
a  = find(pd>1e-13);          % Location of positive probabilities
pd = pd(a);                   % Positive probabilities
d  = d(a);                    % D values with positive probability