znak=1  (* 1 dla fermionów, -1 dla bozonów *)

Energia=p c (* Energia cząstki bezmasowej lub skrajnie relatywistycznej *)

(*Energia = Sqrt[p^2 c^2 + (m c^2)^2]*)

Integrate[1/h^3 Energia 4Pi p^2/(Exp[(Energia-μ)/kT]+znak),{p,0,Infinity}]

ϵ=Refine[%,c/kT>0]

P=znak Refine[kT/h^3 Integrate[Log[1+znak*Exp[(-Energia+μ)/kT]]4Pi p^2,{p,0,Infinity}],c/kT>0&&μ/kT<0]

P/ϵ