What is the state of the art way of implementing double precision special functions? I need the following integral: $$ F_m(t) = \int_0^1 u^{2m} e^{-tu^2} d u = {\gamma(m+{1\over 2}, t)\over 2 t^{m+{1\over 2}}} $$ for $m=0, 1, 2, ...$ and $t>0$, which can be written in terms of the lower incomplete gamma function. Here is my Fortran and C implementation:
https://gist.github.com/3764427
which uses series expansion, sums up the terms until the given accuracy, and then uses recursion relations to efficiently obtain values for lower $m$. I tested it well and I get 1e-15 accuracy for all values of parameters that I need, see the Fortran version's comments for details.
Is there a better way to implement it? Here is a gamma function implementation in gfortran:
https://github.com/mirrors/gcc/blob/master/libgfortran/intrinsics/c99_functions.c#L1781
it is using rational function approximation instead of summing up some infinite series that I am doing. I think that's a better approach, because one should obtain uniform accuracy. Is there some canonical way to approach these things, or does one have to figure out a special algorithm for each special function?
Update 1:
Based on the comments, here is the implementation using SLATEC:
https://gist.github.com/3767621
it reproduces values from my own function, roughly on the level of 1e-15 accuracy. However, I noticed a problem that for t=1e-6 and m=50, the $t^{m+{1\over2}}$ term gets equal to 1e-303 and for higher "m" it simply starts giving incorrect answers. My function doesn't have this problem, because I use a series expansion/recurrence relations directly for $F_m$. Here is an example of a correct value:
$F_{100}$(1e-6)=4.97511945200351715E-003,
but I can't get this using SLATEC because the denominator blows up. As you can see, the actual value of $F_m$ is nice and small.
Update 2:
To avoid the above problem, one can use the function dgamit (Tricomi's incomplete Gamma function), then F(m, t) = dgamit(m+0.5_dp, t) * gamma(m+0.5_dp) / 2, so there is no problem with $t$ anymore, but unfortunately the gamma(m+0.5_dp) blows up for $m\approx 172$. This however might be high enough $m$ for my purposes.
