What is a good way to compute $\Gamma(1/3)$ on a standard pocket calculator?

Question

This question is inspired by this one. The earlier question asks how to calculate a certain integral efficiently with a standard pocket calculator. A fine answer by Travis Willse gives a good result but begs the question of calculating $\Gamma(1/3)$, which appears in the derived series expansion.

So I explored a possible solution. As $n\to\infty$, we have an asymptotic relation $\Gamma(n+a)\sim n^a\Gamma(n)$; so I rendered $a=1/3$ and figured $n=10$ would make this asymptotic relation quite accurate. But, when I approximated $\Gamma(31/3)$ and stepped the recursion relation down to retrieve $\Gamma(1/3)$, I ended with $\approx2.71$ and actually $\Gamma(1/3)\approx2.68$. This was not suitably accurate to meaningfully contribute to Travis Willse's answer.

Now I'm calling in the experts. What are better approximations that (1) use only rational or easily calculator-rendered parameters, (2) yield more digits accurately than were obtained above, (3) do not require an excessive amount of stepping?

Answer 1

As mentioned by heropup, Borwein and Zucker found closed forms for $\Gamma(n/24)$, for integer $n$, using the complete elliptic integral of the first kind; see Wikipedia for details, with references.

Elliptic integrals can be rapidly evaluated to high precision using algorithms based on the AGM, the arithmetic-geometric mean, which converges quadratically.

The AGM is defined as follows. Let $x_0=x, y_0=y$. Then, $$\begin{align} x_{i+1} & = (x_i + y_i) / 2 \\ y_{i+1} & = \sqrt{x_i y_i} \end{align}$$ These sequences rapidly converge on the same value: $$\operatorname{AGM}(x, y) = \lim_{i\to \infty} \, x_i = \lim_{i\to \infty} \, y_i$$

Elliptic integral notation can be confusing. These integrals are sometimes written using the modulus $m$, and sometimes using the parameter $k$, where $m = k^2$.

The complete elliptic integral of the first kind can be calculated as $$K(m) = \frac{\pi/2}{\operatorname{AGM}(1, \sqrt{1-m})}$$

---

Wikipedia gives a couple of equivalent expressions for $\Gamma(1/3)$ in terms of $K()$. Here's a breakdown into a few logical steps.

$$\begin{align} m & = \frac{2 - \sqrt{3}}{4} \\ q & = \frac{\pi/2}{\operatorname{AGM}(1, \sqrt{1-m})} \\ g & = \frac{q \pi \sqrt[3]{128}}{\sqrt[4]{3}} \\ \Gamma\left(\frac{1}{3}\right) & = \sqrt[3]{g} \end{align}$$

As mentioned, the AGM converges quadratically. i.e., the number of correct digits doubles on each loop of the iteration. For the $\Gamma(1/3)$ calculation, two loops are adequate for the precision of a typical pocket scientific calculation, three loops are adequate for standard double-precision floating-point numbers.

Here's a simple implementation of the above algorithm in Python.

Code

""" gamma(1/3), using AGM
    Written by PM 2Ring 2024.06.10
    Pure Python
"""

from itertools import count
from math import sqrt, pi

def agm(g, a, eps=1e-9):
    for i in count(1):
        g, a = sqrt(g * a), (g + a) / 2
        print(i, g, a)
        if a - g < eps:
            break
    return (a + g) / 2

m = (2 - sqrt(3)) / 4
print(f"{m=}")

a = agm(sqrt(1 - m), 1)
print(f"{a=}")

q = pi/2 / a
print(f"{q=}")

g3 = pi * q * (128**(1/3) / 3**(1/4))
g = g3**(1/3)
print(g)

Output

m=0.0669872981077807
1 0.9828152554214186 0.9829629131445341
2 0.9828890815101808 0.9828890842829763
3 0.9828890828965786 0.9828890828965786
a=0.9828890828965786
q=1.59814200211254
2.6789385347077475

Here's a live version of that script, running on the SageMathCell server. And here's a Sage version, using arbitrary precision arithmetic. Of course, Sage also provides the gamma function, the AGM, and a full complement of elliptic integrals and functions.

Here's $\Gamma(1/3)$ using the Sage code, to 300 bits precision.

2. 67893 85347 07747 63365 56929 40974 67764 41286 89377 95730 11009 50428 32759 04176 10167 74381 95409 8289

---

Incidentally, it's possible to calculate arbitrary values of the gamma function using mostly elementary arithmetic (apart from a couple of exponentiations), via the continued fractions of the upper and lower incomplete gamma functions. However, that algorithm doesn't converge quickly, so it's a bit too tedious for manual calculation. I have details & code in this answer.

Answer 2

You can use a better approximation for the ratio of two gamma functions to obtain $$ \Gamma (a) = \frac{{\Gamma (n + a)}}{{a(a + 1) \cdots (a + n - 1)}} \approx \frac{{\Gamma (n)}}{{a(a + 1) \cdots (a + n - 1)}}\left( {n + \frac{{a - 1}}{2}} \right)^a , \quad 2n+1>a. $$ Taking $n=7$ and $a=\frac{1}{3}$ yields $\Gamma(\frac{1}{3})\approx 2.678197$ compared with $\Gamma(\frac{1}{3})=2.678938\ldots$.

Answer 3

You will want to read the article Fast evaluation of the gamma function for small rational fractions using complete elliptic integrals of the first kind by Borwein and Zucker, in which a quadratically convergent AGM is used to compute the elliptic integral.

Presently, I don't have time to work out the specific details for the algorithm for $\Gamma(\tfrac{1}{3})$ from the paper, so I encourage the reader to do so as an exercise. Based on some preliminary numerical computations, only a few AGM terms are needed to achieve high precision.

Answer 4

I have looked into improving the accuracy of my original method by combining the original asymptotic relation

$\Gamma(n+a)\sim n^a\Gamma(n)$

with the Reflection Formula

$\Gamma(a)\Gamma(1-a)=\pi/[\sin(\pi a)].$

The idea is to calculate both $\Gamma(n+a)$ and $\Gamma(n+1-a)$ using the asymptotic relation above for a moderately large $n$, here taken to be $n=10$ as was originally done. With a=1/3 the asymptotic relation above gives

$\Gamma(10+1/3)\approx 781801$

$\Gamma(10+2/3)\approx 1684339$

Stepping down with the recursion relation then gives

$\Gamma(1/3)\approx 2.70903$

$\Gamma(2/3)\approx 1.36915$

which are both correct to only two significant figures. But the errors are of the same sign and similar in magnitude, so they will be accumulated when the product is calculated and compared with the true one from the Reflection Formula:

$2.70703×1.36915=3.70908\tag{6 digits}$

but the Reflection Formula was supposed to give

$\Gamma(1/3)\Gamma(2/3)=\pi\sqrt{4/3}=3.62760\tag{6 digits}$

So if the errors are similar and of the same sign, we should divide the original estimates by the square root of the ratio of these two products, thus normalizing the product. Then

$\Gamma(1/3)\approx 2.67910$

$\Gamma(2/3)\approx 1.35403$

which are now correct to four significant digits (relative error <0.01%) instead of two. We gained an accuracy factor of roughly $100$ versus the original calculation without the reflection formula normalization.

Gamma functions of multiples of $1/24$ can be calculated more rapidly and with higher accuracy with elliptic functions, as described in other answers; the real objective is to explore calculation for some other arguments where the elliptic function option is unavailable/unknown. Multiples of $1/5$ are especially attractive to work with because the triginometric functions in the Reflection Formula may be resolved into combinations of rational operations and square roots. With $n=10$ I get numbers that are all accurate to within 0.01%. Significant digits calculated correctly are shown:

$\Gamma(1/5)\approx 4.591$

$\Gamma(2/5)\approx 2.218$

$\Gamma(3/5)\approx 1.489$

$\Gamma(4/5)\approx 1.164$

The method described here can be converted to a formula that avoids computing the factorials, although the step-down products must still be computed. To wit:

$\Gamma(a)=\lim\limits_{n\to\infty}\sqrt{\dfrac{\pi n^{2a-1}}{\sin(\pi a)}\prod \limits_{k=1}^n\dfrac{k-a}{k-1+a}}$

Truncating this limit at various $n$ values for $\Gamma(1/5)$ reveals that the error is $O(1/n^2)$, with a relatively small coefficient on the $1/n^2$ term because the formula is exact for all half-integral $a$.

Answer 5

Using the method of Borwein and Zucker mentioned by heropup and implemented by PM 2Ring in other answers (see also the answers to this related question) we have $$\Gamma\left(\frac13\right) = \frac{2^{\frac49} \pi^{\frac23}}{3^{\frac1{12}} \operatorname{AGM}\left(1, \sqrt{1 - m}\right)},$$ where $m = \frac{2 - \sqrt 3}{4}$ and $\operatorname{AGM}$ denotes the arithmetic-geometric mean.

Since $m = 0.0669872981\ldots \ll 1$, the arithmetic and geometric means of $1, \sqrt{1 - m}$ are already close, so either is a good approximation for $\operatorname{AGM}\left(1, \sqrt{1 - m}\right)$. Indeed, they are: \begin{alignat}{2} \frac12\left(1 + \sqrt{1 - m}\right) &= \frac{4 + \sqrt 2 + \sqrt 6}{8} &&= \color{#007f00}{0.982}8152554\ldots \\ (1 - m)^{\frac14} &= \frac{\sqrt{1 + \sqrt{3}}}{2^{\frac34}} &&= \color{#007f00}{0.982}9629131\ldots \\ \end{alignat} For the context of the motivating linked problem, the difference of $< 2 \cdot 10^{-4}$ here is more than sufficient.

Approximating the AGM with the geometric mean and simplifying gives the approximation $$\Gamma\left(\frac13\right) \approx \frac{2^{\frac{25}{36}} \pi^{\frac23}}{3^{\frac1{12}} (1 + \sqrt{3})^{\frac16}}.$$ Indeed, we have numerical values \begin{align} \Gamma\left(\frac13\right) &= \color{#007f00}{2.67}89385347\ldots \\ \frac{2^{\frac{25}{36}} \pi^{\frac23}}{3^{\frac1{12}} (1 + \sqrt{3})^{\frac16}} &= \color{#007f00}{2.67}90056121\ldots , \end{align} giving an absolute error of $< 10^{-4}$.

Answer 6

Applying the substitution $x = u^3$ to the usual integral representation of the gamma function gives $$\Gamma\left(\frac13\right) = \int_0^\infty x^{-\frac23} e^{-x} \,dx = 3 \int_0^\infty e^{-y^3} \,dy .$$ The latter integrand decays rapidly—in fact, applying integration by parts with $u = x^{-2}$ gives $\int_M^\infty e^{-y^3} \,dy < \frac13 M^{-2} e^{-M^3}$, so $$\Gamma\left(\frac13\right) - \int_0^2 e^{-y^3} \,dy = \int_2^\infty e^{-y^3} \,du < \frac1{12 e^8} < 10^{-4} .$$ But already for $n = 3$ Simpson's Rule gives $$\Gamma\left(\frac13\right) \approx \int_0^2 e^{-y^3} \,dy \approx \color{#007f00}{2.67}98249579,$$ which has absolute error $< 10^{-3}$, more than good enough for the context of the linked motivating question. One can modestly increase $M$ and $n$ to obtain higher accuracy some while keeping the computations in the realm of feasibility on a pocket calculator.