// Test Fermi
System.out.println("Testing Fermi");
// 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));
}
sbeps=1.e-13;
System.out.printf("Fermi: Maximum discrepancy = %f\n", maxel(vecsub(zz,expect)));
localflag = maxel(vecsub(zz,expect)) > sbeps;
globalflag = globalflag || localflag;
if (localflag) {
fail("*** Fermi: Incorrect function values");
}
// Limiting cases for large eta and k=1/2
maxerr=0.0;
for (i=0;i<N;i++) {
arg=1000.0;
theta=0.1*(i+2);
zz[i]=dirac.val(0.5,arg,theta);
u=1.0+arg*theta;
y=sqrt(SQR(u)-1.0);
x=log((y+sqrt(SQR(y)+4.0))/2.0);
expect[i]=(y*u-x)/pow(2.0*theta,1.5);
err=abs(zz[i]/expect[i]-1.0);