// Test cases with theta=0
// Fermi(k,eta,0)=exp(eta)*LerchPhi(exp(eta),k+1,1)
Fermi dirac = new Fermi();
for (i=0;i<N;i++) {
zz[i]=dirac.val(k[i],eta[i],0.0);
expect[i]=exp(gammln(k[i]+1.0))*exp(eta[i])*lerchphi[i];
// System.out.printf(setw(15) << expect[i] << setw(15) << zz[i];
// System.out.printf(setw(15) << (zz[i]/expect[i]-1.0));
}