In 1978 or 1979, the Hungarian microelectronics company Híradástechnika began manufacturing an OEM version of the Commodore PR-100 programmable calculator. Until then, programmable calculators, privately imported, cost more than the average price of a new car in Hungary. With the PTK-1072, suddenly it was possible for a high-school geek like myself to acquire one. Wow!

But wait... here I was with my PTK-1072 and it couldn't do what my best friend's lousy non-programmable Commodore SR-4190R was capable of: it could not calculate the Gamma function.

Needless to say, in 1979 I had no access to resources such as the excellent book Numerical Recipes in C. However, I found another book on differential and integral calculus (A differenciál- és integrálszámítás elemei, by Pál Szász, Budapest, 1955) that contained an extensive discussion on the Gamma function.

As part of this discussion, the book provided a variation of the Gamma function's Maclaurin series:

\[2\ln\Gamma(1+x)=\ln\frac{1-x}{1+x}+\ln\frac{\pi x}{\sin\pi x}+2(1-C)x+2\sum\limits_{n=1}^\infty\left[\left(1-\sum\limits_{\nu=1}^\infty\frac{1}{\nu^{2n+1}}\right)\frac{x^{2n+1}}{2n+1}\right].\]

This approximation converges fairly quickly for $0.5\lt x\lt 1.5$.

$C$ is Euler's constant, defined as the infinite sum

\[\lim\limits_{N\rightarrow\infty}\left(\sum\limits_{n=1}^N\frac{1}{n}-\ln N\right)=0.577~215~664~901~532~860~607...\]

Values of this constant can be easily obtained from the literature. That is not the case for the remaining coefficients in this formula. My poor PTK-1072 ran for several weeks to obtain an approximation of \(\sum\limits_{\nu=2}^\infty 1/\nu^3\) to 8 digits of precision. In the end, I obtained the following values for the various coefficients:

\begin{align}x:&~0.8455686582,\\

x^3:&~{-0.1347045888},\\

x^5:&~{-1.477110203}\times 10^{-2},\\

x^7:&~{-2.38550782}\times 10^{-3},\\

x^9:&~{-4.463095166}\times 10^{-4},\\

x^{11}:&~{-8.985247346}\times 10^{-5}.\end{align}

As it turns out, my calculations were not entirely accurate. Being equipped with a dual-processor Pentium-II workstation, I no longer find it problematic to calculate these coefficients to a high degree of precision. Whereas the PTK-1072 ran for weeks to add a sum of a few thousand values, I can add several million values in a matter of seconds. Here are my modern results:

\begin{align}x:&~0.8455686702,\\

x^3:&~{-0.1347046021},\\

x^5:&~{-1.477110206}\times 10^{-2},\\

x^7:&~{-2.385507823}\times 10^{-3},\\

x^9:&~{-4.463095169}\times 10^{-4},\\

x^{11}:&~{-8.985247348}\times 10^{-5}.\end{align}

As I mentioned, the formula above can be used to obtain the logarithm of the Gamma function for arguments between 0.5 and 1.5. The function's value for other arguments can be calculated using the following relationship:

\[\Gamma(x)\Gamma(1-x)=\frac{\pi x}{\sin\pi x}.\]

With this and the recurrence relationship, it is now possible to calculate the Gamma function for any real argument. (I have never investigated whether the Maclaurin series can be used for any complex argument.)

Sadly, I no longer have my old notebook with my PTK-1072 programs. However, I do have a TI-59 magnetic card, still readable, that contains a variant of my Gamma function algorithm. From this, I was able to reconstruct my original program, shown below.

Note that you need to set memory registers 4-9 to the values specified before running the program (the coefficients shown in the listing below are the ones I calculated in 1979.) To enter values with ten-digit precision, use key sequences like this one:

0.1347045 +/− − 8.88 EE 8 +/− = M 4

In order for the program to work properly, you must first set the calculator to radian mode by pressing F rad. To use the function, enter the argument then hit R/S.

M4 = −.1347045888
M5 = −1.477110203e−2
M6 = −2.38550782e−3
M7 = −4.463095166e−4
M8 = −8.985247346e−5
M9 = 0.8455686582

85 00	−
81 01	1
74 02	×
51 03	M
91 04	0
95 05	=
51 06	M
81 07	1
74 08	×
52 09	MR
62 10	8
84 11	+
52 12	MR
61 13	7
74 14	×
52 15	MR
81 16	1
84 17	+
52 18	MR
73 19	6
74 20	×
52 21	MR
81 22	1
84 23	+
52 24	MR
72 25	5
74 26	×
52 27	MR
81 28	1
84 29	+
52 30	MR
71 31	4
74 32	×
52 33	MR
81 34	1
84 35	+
52 36	MR
63 37	9
74 38	×
52 39	MR
91 40	0
84 41	+
64 42	(
45 43	π
74 44	×
52 45	MR
91 46	0
75 47	÷
55 48	x-y
22 49	sin
74 50	×
64 51	(
81 52	1
85 53	−
52 54	MR
91 55	0
65 56	)
75 57	÷
64 58	(
52 59	MR
91 60	0
84 61	+
81 62	1
65 63	)
65 64	)
32 65	ln
95 66	=
21 67	F
32 68	ex
35 69	√