I couldn't find any implementations I liked, so here you can find some programs I wrote in BASIC (with included source code) to compute various special functions. Unless otherwise specified, these programs are subject to the GNU General Public License as published by the Free Software Foundation, either version three or, at your option, any later version. These files contain their own source code and are the preferred form for modification.
Files with the extension .8Xp are designed for use on calculators in the TI-83/TI-84 Plus family,
and files ending in .v2f are designed for calculators like the Voyage 200, the TI-89 Titanium,
and others in the same family.
GAMMA.8Xp and main.gamma.v2f:
Gamma function with near-exact precision for real and complex values
The gamma function is perhaps the most well-known special function which obeys the definition \( \Gamma(z) = \int_0^\infty t^{z - 1} e^{-t} \, \mathrm{d}t \) for complex \( z \) with a real part greater than \( 0 \). This function is defined elsewhere by analytic continuation. This program uses a specially-crafted version of Stirling's approximation that works even for medium-sized values of the parameter, and uses the defining property of the gamma function \( x \Gamma(x) = \Gamma(x + 1) \) to scale up smaller values of the parameter so this approximation can still be used. Specifically, we use a correction term to Stirling's approximation which is the Padé approximant to the Stirling series, in particular the formula \[ \Gamma(z) = \sqrt{\frac{2 \pi}{z}} \left(\frac{z}{e}\right)^z e^{\frac{210z^2 + 53}{2520z^3 + 720z}} (1 + O(z^{-7})) \quad (z \to \infty) \]
EIEXPINT.8Xp and main.ei.v2f:
Exponential integral with near-exact precision for real parameters greater than zero
These programs compute the exponential integral \[ \mathrm{Ei}(x) = ⨍_{-\infty}^x \frac{e^t}{t} \, \mathrm{d}t \] for real positive \( x \), where the integral is defined to be Cauchy principal value integral around the singularity at \( t = 0 \). The computation is done with a slightly modified version of the exponential integral that is appropriate for direct numerical integration for small values, and the asymptotic expansion is used for larger values.
Support for the exponential integral \( E_1(z) = \int_z^\infty \frac{e^{-t}}{t} \, \mathrm{d}t \) as well as the sine and cosine integrals for real and complex parameters will be coming soon.
The same conventions as before apply to these programs.
main.eulersum.v2f: Euler summationThis program performs Euler summation, also known as the Euler transform, to sum infinite series having terms with changing signs. It's useful both for slowly-converging series as well as divergent series. The transformation is \[ \sum_{n = 0}^\infty a_n = \sum_{n = 0}^\infty \frac{\sum_{k = 0}^n {n \choose k} a_k}{2^{n + 1}} \] Here is a fun exercise to "confirm" the validity of the transformation. (If you're familiar with symbolic expressions for the Bernoulli numbers, this should be familiar.) Pretend that the subscripts on \( a_n \) can be switched out for superscripts, like \( a_n \to a^n \), and then try to simplify the expression. You will find that the Euler transform expression is "equivalent" to the geometric series \( \frac{1}{1 - a} \) which, when re-expanded, is equal to the sum \( \sum_{n = 0}^\infty a_n \).