This question was previously posted to Math.SE here and had received no answers at the time of this posting.
When performing computational work, I often come across a univariate function, defined in terms of an integral or differential equation, which I would like to rapidly evaluate (say, millions of times per second) over a specified interval to a given precision (say, one part in $10^{10}$). For example, the function $$ f(\alpha) = \int_{k=0}^\infty \frac{e^{-\alpha^2 k^2}}{k+1}\ \mathrm{d}k $$ over the interval $\alpha \in (0,10)$ came up in a recent project. Now it happens that this integral can be evaluated in terms of standard special functions (in particular, $\operatorname{Ei}(z)$ and $\operatorname{erfi}(z)$), but suppose we had a much more complicated function for which no such evaluation was known. Is there a systematic technique I can apply to develop my own numerical routines for the evaluation of such functions?
I am sure plenty of techniques must be out there, as fast algorithms seem to exist for basically all of the common special functions. However, I emphasize that the sort of technique I am looking for should not rely on the function having some particular structure (e.g. recurrence relations like $\Gamma(n+1) = n\Gamma(n)$ or reflection formulas like $\Gamma(z) \Gamma(1-z) = \pi \csc(\pi z)$). Ideally, such a technique would work for just about any (sufficiently well-behaved) function I come across.
You can take for granted that I do have some slow method of evaluating the desired function (e.g. direct numerical integration) to any precision, and that I'm willing to do a lot of pre-processing work with the slow method in order to develop a fast method.
