% MOMFERMIDIRAC  Generates the first 2N moments to dig 
%    decimal places of the Fermi-Dirac weight function 
%    w(x)=[1/(e^x+1)]^r on R_{+} with r a positive integer.
%
function mf=momfermidirac(dig,N,r)
digits(dig); 
m(1,1)=vpa('log(2)');
if r==1
  mf(1)=m(1,1);
else
  for rho=1:r-1
    m(1,rho+1)=vpa(m(1,rho)-1/(rho*(2^rho)));
  end
  mf(1)=m(1,r);
end
for k=2:2*N
  m(k,1)=(1-vpa(2^(1-k)))*gamma(vpa(k))*zeta(vpa(k));
  if r==1
    mf(k)=m(k,1);
  end
end
if r>1
  for k=2:2*N
    for rho=1:r-1
      m(k,rho+1)=vpa(m(k,rho)-(k-1)*m(k-1,rho)/rho);
    end
    mf(k)=m(k,r);
  end
end