I am finding many situations where I have to numerically integrate some function $f(x)$ of the form:
$$f(x)=f_{s}(x)-ax^{-n},$$
where $f_s$ is a special function with a finite-order pole that is canceled by the $-ax^{-n}$ term. I.e., $f_s(x)$ has Laurent expansion
$$f_s(x)=ax^{-n}+a_0+a_1x+a_2x^2+\cdots.$$
As an example, consider $f_s(x)=E_{-n}(x),$ or fs[x_]:=ExpIntegralE[-n,x]
, and $f(x)=E_{-n}(x)-(n!)x^{-n},$ which has a removable singularity at $x=0.$
In most cases I have found that integrating directly over $f(x)$ works, but sometimes I have to crank up WorkingPrecision
fairly high in order for the cancellation near the pole to work out. And I think this method of computation must be inefficient, because Mathematica must compute many digits of precision that only end up canceling.
So is there a better/more computationally efficient strategy for solving these integrals in general? In particular, I wonder if there is a more efficient way of computing $f(x)$ near the removable poles in general, even before integration comes into play. Since Mathematica has built-in algorithms to calculate special functions to arbitrary precision (presumably as efficiently as possible), I wonder if it has smart ways of computing only the "non-pole" part of special functions near their poles?