In the interesting range, the integrand is always monotonic, so numerical integration is (or should be) straightforward. I have to admit that in zipfR I need finite gamma integrals, which I currently compute as diff of incomplete gamma function for end points … exactly, catastrophic cancellation.