# How do I develop numerical routines for the evaluation of my own special functions?

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.

• I would have to say that I think the responses you got on MSE from Hagen von Eitzen and from njuffa are probably right — I think it's rare for there to be a "special" technique. If you look at the multitude of "Methods of Computation" sections of the DLMF, it often recommends reducing a problem to evaluating an ODE/integral/series/etc. I also agree with njuffa that polynomial/rational approximation and argument transformations is the best first thing to try. For your example function the first step is to transform it to get rid of the singularity at $\alpha=0$. – Kirill Nov 24 '15 at 23:13
• The question seems somewhat broad. Would you say that argument transformation together with polynomial/rational approximation is the kind of a systematic technique that is a suitable answer to this question? – Kirill Nov 24 '15 at 23:15
• @Kirill I was hoping there would be more powerful general-purpose techniques for constructing fast function approximations, but if that's the best one can do at this level of generality, then I'd gladly accept that answer. – David Zhang Nov 24 '15 at 23:18
• @Kirill In particular, what systematic techniques are there for constructing high-accuracy polynomial or rational approximations to a function for a given error tolerance over an interval? I know these approximations are often easier to do piecewise; are there systematic means for deciding how to subdivide the interval of interest? – David Zhang Nov 24 '15 at 23:24
• If this question was answered over on MSE, would you mind adding a note to this question here? – Kirill Nov 26 '15 at 1:08