T_{E} X

Numerical Recipe for Calculating the Gamma Function

Preface

I wanted to do some remedial mathematics. I enjoy doing math and I use it a little bit in my work but not enough to be sharp. To refresh my memory I started with Abramowitz and Stegun's Handbook of Mathematical Functions.

I started with Chapter 3 Elementary analytical methods and reviewed the binomial theorem and the rules of differentiation and integration because I needed to re-learn these subjects. I needed this background to understand the algorithms which I use for modeling electrical circuit and calibrating the test equipment at Saunders & Associates manufacturing firm. We design and manufacture automated test equipment. One specialty of our products is that the algorithms we use are the best possible.

I started using my numerical algebra program DERIVE ``A Mathematical Assistant Program'' to aid me in practicing elementary algorithms. Some of the solutions involved the gamma function. DERIVE generates the gamma function using the algorithm in NUMERICAL RECIPES by Press, Flannery, Teukolsy and Vetterling. The numerical values produced by that algorithm are good to thirteen decimal digits. I needed an algorithm to go to a larger precision. Two other programs which do numerical algebra, Maple and MACSYMA, do contain the gamma function with adequate precision. I did not find those programs as easy to use as DERIVE. Also I wanted to be able to do the algorithms myself in both BASIC and LISP so that I could use them as tools for making our products at Saunders.

Now the fun began. I called Albert Rich and David Stoutemyer the authors of DERIVE to ask if they knew of any methods for calculating the gamma function. They referred me to NUMERICAL RECIPES, Abramowitz and Stegun and A First Course in Numerical Analysis by Ralston and Rabinowitz. None of them had the algorithm I needed.

Then I asked the University of Waterloo who produced MAPLE what algorithms they used to produce the gamma function. After talking to several people I concluded none of them knew the algorithm. Back to square ONE.

My discussion with Richard Petti at Macsyma Inc. the producers of MACSYMA revealed that he did not know the methods used but he would ask the keeper, Bill Gosper. The outcome of that is Bill Gosper gave me hints on how to get started. Though the exact algorithms are proprietary to MACSYMA, Bill Gosper's hints were adequate to get me started. This monograph is the result.

The subject matter in this monograph is the numerical methods to generate the gamma function, related functions and the numbers and constants necessary. The related functions are digamma function, polygamma function, and the Riemann zeta function. The numbers required are the Bernoulli and Euler Numbers. The constants required are Pi and Euler's constant.

The techniques used include the Euler-Maclaurin summation formula and Taylor's series. Also methods to produce the logarithmic and exponential functions are included to make this monograph independent of the system on which it may be implemented.

A definition of the gamma function is:

Γ (z) = lim_{n \to \infty} \frac{n! n^{z}}{z (z + 1) \dots (z + n)}

Put it in a different form.

Γ (z) = \frac{1}{z} lim_{n \to \infty} \frac{n^{z}}{Π_{k = 1}^{n} \frac{z}{k} + 1}

Take the logramithm to convert the product into a series.

\ln (Γ (z)) = - \ln (z) + lim_{n \to \infty} z \ln (n) - \sum_{k = 1}^{n} \ln (\frac{z}{k} + 1)

This limit converges very slowly so it can be helped by the Euler-Maclaurin summation formula. The Euler-Maclaurin summation formula is:

\sum_{k = 0}^{m} F (a + k h) = \frac{1}{h} \int_{a}^{b} F (t) dt + \frac{1}{2} (F (b) + F (a)) + \sum_{k = 1}^{n - 1} \frac{h^{2 k - 1}}{(2 k)!} B_{2 k} (F^{(2 k - 1)} (b) - F^{(2 k - 1)} (a)) + e_{n}

B_{2 k}

is a Bernoulli number,

b = a + m h

, and

e_{n}

is the error. The summation can be split into two parts.

\sum_{k = 1}^{n} \ln (\frac{z}{k} + 1) = \sum_{k = 1}^{a - 1} \ln (\frac{z}{k} + 1) + \sum_{k = a}^{n} \ln (\frac{z}{k} + 1)

Put the second part in a form to use the Euler-Maclaurin summation formula.

\sum_{k = a}^{n} \ln (\frac{z}{k} + 1) = \sum_{k = 0}^{n - a} \ln (\frac{z}{a + k} + 1)

m = n - a, h = 1, F (a + k h) = \ln (\frac{z}{a + k h} + 1)

Define the function.

F (x) = \ln (\frac{z}{x} + 1)

Find the derivatives of the function.

F^{(1)} (x) = \frac{1}{x + z} - \frac{1}{x}

F^{(3)} (x) = \frac{2}{(x + z)^{3}} - \frac{2}{x^{3}}

:

F^{(2 k - 1)} (x) = (2 k - 2)! ((x + z)^{- 2 k + 1} - x^{- 2 k + 1})

Calculate the integral of the function.

\int_{a}^{b} F (t) dt = (z + b) \ln (z + b) - b \ln (b) - (z + a) \ln (z + a) + a \ln (a)

Substitute the three parts into the summation. The Euler-Maclaurin summation becomes:

\sum_{k = 0}^{m} F (a + k) = (z + b) \ln (z + b) - b \ln (b) - (z + a) \ln (z + a) + a \ln (a) + \frac{1}{2} (F (b) + F (a)) +

\sum_{k = 1}^{n - 1} \frac{B_{2 k}}{(2 k)!} (2 k - 2)! ((b + z)^{- 2 k + 1} - b^{- 2 k + 1} - (a + z)^{- 2 k + 1} + a^{- 2 k + 1}) + e_{n}

Change the form of the equation.

\sum_{k = a}^{b} F (k) = (z + b) \ln (z + b) - b \ln (b) - (z + a) \ln (z + a) + a \ln (a) + \frac{1}{2} (F (b) + F (a)) +

\sum_{k = 1}^{n - 1} \frac{B_{2 k}}{(2 k - 1) (2 k)} ((b + z)^{- 2 k + 1} - b^{- 2 k + 1} - (a + z)^{- 2 k + 1} + a^{- 2 k + 1}) + e_{n}

Add more terms to the summation.

\sum_{k = 1}^{b} \ln (\frac{z}{k} + 1) = \sum_{k = 1}^{a - 1} \ln (\frac{z}{k} + 1) + (z + b) \ln (z + b) - b \ln (b) - (z + a) \ln (z + a) + a \ln (a) + \frac{1}{2} (F (b) + F (a)) +

\sum_{k = 1}^{n - 1} \frac{B_{2 k}}{(2 k - 1) (2 k)} ((b + z)^{- 2 k + 1} - b^{- 2 k + 1} - (a + z)^{- 2 k + 1} + a^{- 2 k + 1}) + e_{n}

Change the form of the logarithm of the gamma function.

\ln (Γ (z)) = - \ln (z) + lim_{b \to \infty} z \ln (b) - \sum_{k = 1}^{b} \ln (\frac{z}{k} + 1)

Substitute the Euler-Maclaurin summation into the split series for the logarithm of the gamma function. The logarithm of the gamma function becomes:

\ln (Γ (z)) = - \ln (z) + lim_{b \to \infty} z \ln (b) - \sum_{k = 1}^{a - 1} \ln (\frac{z}{k} + 1) - (z + b) \ln (z + b) + b \ln (b) + (z + a) \ln (z + a) - a \ln (a) -

\frac{1}{2} (F (b) + F (a)) - \sum_{k = 1}^{n - 1} \frac{B_{2 k}}{(2 k - 1) (2 k)} ((b + z)^{- 2 k + 1} - b^{- 2 k + 1} - (a + z)^{- 2 k + 1} + a^{- 2 k + 1}) + e_{n}

Take the limit.

lim_{b \to \infty} z \ln (b) - (z + b) \ln (z + b) + b \ln (b) = - z

lim_{b \to \infty} \frac{1}{2} F (b) = 0

lim_{b \to \infty} ((b + z)^{- 2 k + 1} - b^{- 2 k + 1}) = 0

Simplify the equation.

\ln (Γ (z)) = - \ln (z) - \sum_{k = 1}^{a - 1} \ln (\frac{z}{k} + 1) - z + (z + a) \ln (z + a) - a \ln (a) - \frac{\ln (\frac{z}{a} + 1)}{2} +

\sum_{k = 1}^{n - 1} \frac{B_{2 k}}{(2 k - 1) (2 k)} ((a + z)^{- 2 k + 1} - a^{- 2 k + 1}) + e_{n}

Take the exponential function of both sides.

Γ (z) = \frac{(z + a)^{z + a}}{z a^{a} \sqrt{\frac{z}{a} + 1} Π_{k = 1}^{a - 1} (\frac{z}{k} + 1)}

\exp (- z + \sum_{k = 1}^{n - 1} \frac{B_{2 k}}{(2 k - 1) (2 k)} ((a + z)^{- 2 k + 1} - a^{- 2 k + 1}) + e_{n})

This very complex formula is an asymptotic series. It gives very good results for even large arguments with very few terms. It requires numbers that are very large. DERIVE can handle these numbers. This series requires the proper value of

a

to be used. If

a

is too large the evaluation will take too long and if

a

is too small the series will diverge before an accurate answer is calculated.

The Euler-Maclaurin summation formula requires Bernoulli numbers. This is how to generate the Bernoulli numbers.

B_{0} = 1, B_{1} = - 1 / 2, B_{2} = 1 / 6

There are special values for the Bernoulli numbers.

B_{2 n + 1} = 0 (n = 1, 2, \dots)

This is a recursive formula to generate Bernoulli numbers.

B_{2 n} = - \frac{1}{2 n + 1} + 1 / 2 - \frac{n}{6} - \sum_{k = 2}^{n - 1} \frac{B_{2 k}}{(2 k)!} Π_{j = 0}^{2 k - 2} (2 n - j) (n = 2, 3, \dots)

There are special values of the gamma function that do not need this complex formula.

The gamma function is defined by the limit:

Γ (z) = lim_{n \to \infty} \frac{n! n^{z}}{Π_{k = 0}^{n} z + k} (z \neq 0, - 1, - 2, \dots)

There are the values not covered by the limit definition.

lim_{z \to n} \frac{1}{Γ (- z)} = 0 (n = 0, 1, 2, \dots)

There are the integer arguments of the gamma function. First the value one.

Γ (1) = lim_{n \to \infty} \frac{n! n^{1}}{Π_{k = 0}^{n} 1 + k}

The product in the denominator divides into the factorial in the numerator.

Γ (1) = lim_{n \to \infty} \frac{n}{1 + n}

Now the limit is obvious.

Γ (1) = 1

The second value is two.

Γ (2) = lim_{n \to \infty} \frac{n! n^{2}}{Π_{k = 0}^{n} 2 + k}

The product in the denominator divides into the factorial in the numerator.

Γ (2) = lim_{n \to \infty} \frac{1! n^{2}}{(1 + n) (2 + n)}

Now the limit is obvious.

Γ (2) = 1!

The third value is three.

Γ (3) = lim_{n \to \infty} \frac{n! n^{3}}{Π_{k = 0}^{n} 3 + k}

The product in the denominator divides into the factorial in the numerator.

Γ (3) = lim_{n \to \infty} \frac{2! n^{3}}{(1 + n) (2 + n) (3 + n)}

Now the limit is obvious.

Γ (3) = 2!

Now I will make a wild guess for all the rest of the integer arguments.

:

Γ (1 + n) = n! (n = 1, 2, \dots)

A recurrence formula for the gamma function can be found. The definition of the gamma function is:

Γ (z) = lim_{n \to \infty} \frac{n! n^{z}}{Π_{k = 0}^{n} z + k}

Add one to the argument.

Γ (1 + z) = lim_{n \to \infty} \frac{n! n^{1 + z}}{Π_{k = 0}^{n} 1 + z + k}

Factor out

n

in the numerator and

z + n + 1

in the denominator.

Γ (1 + z) = lim_{n \to \infty} \frac{n! n n^{z}}{(z + n + 1) Π_{k = 1}^{n} z + k}

Multiply numerator and denomiator by

z

Γ (1 + z) = z lim_{n \to \infty} \frac{n! n n^{z}}{(z + n + 1) Π_{k = 0}^{n} z + k}

Take the limit of

\frac{n}{z + n + 1}

lim_{n \to \infty} \frac{n}{z + n + 1} = 1

Substitute this result.

Γ (1 + z) = z lim_{n \to \infty} \frac{n! n^{z}}{Π_{k = 0}^{n} z + k}

Substitute for the definition of the gamma function.

Γ (1 + z) = z Γ (z)

This is a recurrence formula for the gamma function. It is very useful to reduce real arguments to values near one. It does not help to reduce large complex arguments.

Gauss' multiplication formula is very useful.

Γ (n z) = (2 π)^{\frac{(1 - n)}{2}} n^{n z - \frac{1}{2}} Π_{k = 0}^{n - 1} Γ (z + \frac{k}{n})

I did not try to derive this formula. It would be interesting to see it derived. It can be used to get the duplication formula with

n = 2

Γ (2 z) = (2 π)^{- \frac{1}{2}} 2^{2 z - \frac{1}{2}} Γ (z) Γ (z + \frac{1}{2})

Γ (2 z) = π^{- \frac{1}{2}} 2^{2 z - 1} Γ (z) Γ (z + \frac{1}{2})

It can be used to get the value of

Γ (\frac{1}{2})

with

n = 2

and

z = \frac{1}{2}

Γ (1) = π^{- \frac{1}{2}} 2^{1 - 1} Γ (\frac{1}{2}) Γ (1)

1 = π^{- \frac{1}{2}} Γ (\frac{1}{2})

Γ (\frac{1}{2}) = π^{\frac{1}{2}}

It can be used to get the triplication formula with

n = 3

Γ (3 z) = (2 π)^{- 1} 3^{3 z - \frac{1}{2}} Γ (z) Γ (z + \frac{1}{3}) Γ (z + \frac{2}{3})

γ

is Euler's constant which is defined as:

γ = lim_{m \to \infty} [\sum_{k = 1}^{m} \frac{1}{k} - \ln (m)]

This limit converges very slowly. The Euler-Maclaurin summation formula is a way to help evaluate slowly converging series. The Euler-Maclaurin summation formula is:

\sum_{k = 0}^{m} F (a + k h) = \frac{1}{h} \int_{a}^{b} F (t) dt + \frac{1}{2} (F (b) + F (a)) + \sum_{k = 1}^{n - 1} \frac{h^{2 k - 1}}{(2 k)!} B_{2 k} (F^{(2 k - 1)} (b) - F^{(2 k - 1)} (a)) + e_{n}

B_{2 k}

is a Bernoulli number,

b = a + m h

, and

e_{n}

is the error.

First split the limit into a finite summation and a limit.

γ = lim_{b \to \infty} [\sum_{k = 1}^{a - 1} \frac{1}{k} + \sum_{k = a}^{b} \frac{1}{k} - \ln (b)]

The limit of the finite summation is known.

γ = \sum_{k = 1}^{a - 1} \frac{1}{k} + lim_{m \to \infty} [\sum_{k = a}^{b} \frac{1}{k} - \ln (b)]

Change the variable in the finite sum.

\sum_{k = a}^{b} \frac{1}{k} = \sum_{k = 0}^{b - a} \frac{1}{a + k}

Since

b = a + m

substitute for

b

\sum_{k = a}^{b} \frac{1}{k} = \sum_{k = 0}^{m} \frac{1}{a + k}

Next try to use the Euler-Maclaurin summation formula.

h = 1 F (a + k h) = \frac{1}{a + k}

I will put Euler's constant in a form to fit the Euler-Maclaurin summation formula.

\sum_{k = 0}^{m} F (a + k) = \int_{a}^{b} F (t) dt + \frac{1}{2} (F (b) + F (a)) + \sum_{k = 1}^{n - 1} \frac{B_{2 k}}{(2 k)!} (F^{(2 k - 1)} (b) - F^{(2 k - 1)} (a)) + e_{n}

The integral can be evaluated.

\int_{a}^{b} F (t) dt = \int_{a}^{b} \frac{1}{t} dt

It is just the logarithm.

\int_{a}^{b} F (t) dt = \ln (b) - \ln (a)

The next term is simple also.

\frac{1}{2} (F (b) + F (a)) = \frac{1}{2} (\frac{1}{b} + \frac{1}{a})

The derivative in the summation easy to generate.

F^{1} (x) = - \frac{1}{x^{2}}

F^{2} (x) = \frac{2}{x^{3}}

:

F^{n} (x) = (- 1)^{n} \frac{n!}{x^{n + 1}}

Put the derivative in the proper form by changing variables.

F^{2 k - 1} (x) = - \frac{(2 k - 1)!}{x^{2 k}}

The three different parts can be substituted in the Euler-Maclaurin summation formula.

\sum_{k = 0}^{m} F (a + k) = \ln (b) - \ln (a) + \frac{1}{2} (\frac{1}{b} + \frac{1}{a}) + \sum_{k = 1}^{n - 1} \frac{B_{2 k}}{(2 k)!} (- \frac{(2 k - 1)!}{b^{2 k}} + \frac{(2 k - 1)!}{a^{2 k}}) + e_{n}

Divide by the factorial in the denominator.

\sum_{k = 0}^{m} F (a + k) = \ln (b) - \ln (a) + \frac{1}{2} (\frac{1}{b} + \frac{1}{a}) - \sum_{k = 1}^{n - 1} \frac{B_{2 k}}{2 k} (\frac{1}{b^{2 k}} - \frac{1}{a^{2 k}}) + e_{n}

Substitute the Euler-Maclaurin summation formula in the limit part of Euler's constant.

γ = \sum_{k = 1}^{a - 1} \frac{1}{k} + lim_{m \to \infty} [\ln (b) - \ln (a) + \frac{1}{2} (\frac{1}{b} + \frac{1}{a}) - \sum_{k = 1}^{n - 1} \frac{B_{2 k}}{2 k} (\frac{1}{b^{2 k}} - \frac{1}{a^{2 k}}) + e_{n} - \ln (b)]

Evaluate the limit.

γ = \sum_{k = 1}^{a - 1} \frac{1}{k} - \ln (a) + \frac{1}{2 a} + \sum_{k = 1}^{n - 1} \frac{B_{2 k}}{2 k a^{2 k}} + e_{n}

This is a series approximation to Euler's constant which will give good results with very few terms.

The digamma function is defined as:

ψ (z) = \frac{d \ln (Γ (z))}{dz}

Substituting the logramithm of the gamma function gives a limit for the digamma function.

ψ (z) = \frac{d (- \ln (z) + lim_{n \to \infty} z \ln (n) - \sum_{k = 1}^{n} \ln (\frac{z}{k} + 1))}{dz}

Simplify the equation.

ψ (z) = - \frac{1}{z} + lim_{n \to \infty} \ln (n) - \sum_{k = 1}^{n} \frac{1}{k + z}

This limit converges very slowly so it can be helped by the Euler-Maclaurin summation formula. The Euler-Maclaurin summation formula is:

\sum_{k = 0}^{m} F (a + k h) = \frac{1}{h} \int_{a}^{b} F (t) dt + \frac{1}{2} (F (b) + F (a)) + \sum_{k = 1}^{n - 1} \frac{h^{2 k - 1}}{(2 k)!} B_{2 k} (F^{(2 k - 1)} (b) - F^{(2 k - 1)} (a)) + e_{n}

B_{2 k}

is a Bernoulli number,

b = a + m h

, and

e_{n}

is the error. The summation can be split into two parts.

\sum_{k = 1}^{n} \frac{1}{k + z} = \sum_{k = 1}^{a - 1} \frac{1}{k + z} + \sum_{k = a}^{n} \frac{1}{k + z}

Put the second part in a form to use the Euler-Maclaurin summation formula.

\sum_{k = a}^{n} \frac{1}{k + z} = \sum_{k = 0}^{n - a} \frac{1}{a + k + z}

m = n - a, h = 1, F (a + k h) = \frac{1}{a + k + z}

Define the function.

F (x) = \frac{1}{x + z}

Find the derivatives of the function.

F^{(1)} (x) = - \frac{1}{(x + z)^{2}}

F^{(3)} (x) = - \frac{6}{(x + z)^{4}}

:

F^{(2 k - 1)} (x) = - (2 k - 1)! (x + z)^{- 2 k}

Calculate the integral of the function.

\int_{a}^{b} F (t) dt = \ln (z + b) - \ln (z + a)

Substitute the three parts into the summation. The Euler-Maclaurin summation becomes:

\sum_{k = 0}^{m} F (a + k) = \ln (z + b) - \ln (z + a) + \frac{1}{2} (F (b) + F (a)) - \sum_{k = 1}^{n - 1} \frac{B_{2 k}}{(2 k)!} (2 k - 1)! ((b + z)^{- 2 k} - (a + z)^{- 2 k}) + e_{n}

Change the form of the equation.

\sum_{k = a}^{b} F (k) = \ln (z + b) - \ln (z + a) + \frac{1}{2} (F (b) + F (a)) - \sum_{k = 1}^{n - 1} \frac{B_{2 k}}{2 k} ((b + z)^{- 2 k} - (a + z)^{- 2 k}) + e_{n}

Add more terms to the summation.

\sum_{k = 1}^{b} \frac{1}{k + z} = \sum_{k = 1}^{a - 1} \frac{1}{k + z} + \ln (z + b) - \ln (z + a) + \frac{1}{2} (F (b) + F (a)) - \sum_{k = 1}^{n - 1} \frac{B_{2 k}}{2 k} ((b + z)^{- 2 k} - (a + z)^{- 2 k}) + e_{n}

Change the form of the digamma function.

ψ (z) = - \frac{1}{z} + lim_{b \to \infty} \ln (b) - \sum_{k = 1}^{b} \frac{1}{k + z}

Substitute the Euler-Maclaurin summation into the split series for the digamma function. The digamma function becomes:

ψ (z) = - \frac{1}{z} + lim_{b \to \infty} \ln (b) - \sum_{k = 1}^{a - 1} \frac{1}{k + z} - \ln (z + b) + \ln (z + a) - \frac{1}{2} (F (b) + F (a)) +

\sum_{k = 1}^{n - 1} \frac{B_{2 k}}{2 k} ((b + z)^{- 2 k} - (a + z)^{- 2 k}) + e_{n}

Take the limit.

lim_{b \to \infty} \ln (b) - \ln (z + b) = 0

lim_{b \to \infty} \frac{1}{2} F (b) = 0

lim_{b \to \infty} (b + z)^{- 2 k} = 0 k > 0

Simplify the equation.

ψ (z) = - \frac{1}{z} - \sum_{k = 1}^{a - 1} \frac{1}{k + z} + \ln (a + z) - \frac{1}{2 (a + z)} - \sum_{k = 1}^{n - 1} \frac{B_{2 k}}{2 k (a + z)^{2 k}} + e_{n}

The polygamma function is defined as:

ψ^{(n)} (z) = \frac{d^{n}}{{dz}^{n}} ψ (z)

Change the form of the digamma function definition.

ψ (z) = lim_{b \to \infty} \ln (n) - \sum_{k = 0}^{b} \frac{1}{k + z}

Substituting the limit for the digamma function gives the limit for the polygamma function.

ψ^{(n)} (z) = lim_{b \to \infty} (- 1)^{n + 1} n! \sum_{k = 0}^{b} \frac{1}{(k + z)^{n + 1}}

Take the limit.

ψ^{(n)} (z) = (- 1)^{n + 1} n! \sum_{k = 0}^{\infty} \frac{1}{(k + z)^{n + 1}} (z \neq 0, - 1, - 2, - 3, \dots)

This series converges very slowly so it can be helped by the Euler-Maclaurin summation formula. The Euler-Maclaurin summation formula is:

\sum_{k = 0}^{m} F (a + k h) = \frac{1}{h} \int_{a}^{b} F (t) dt + \frac{1}{2} (F (b) + F (a)) + \sum_{k = 1}^{n - 1} \frac{h^{2 k - 1}}{(2 k)!} B_{2 k} (F^{(2 k - 1)} (b) - F^{(2 k - 1)} (a)) + e_{n}

B_{2 k}

is a Bernoulli number,

b = a + m h

, and

e_{n}

is the error. The summation can be split into two parts.

\sum_{k = 0}^{\infty} \frac{1}{(k + z)^{n + 1}} = \sum_{k = 0}^{a - 1} \frac{1}{(k + z)^{n + 1}} + \sum_{k = a}^{\infty} \frac{1}{(k + z)^{n + 1}}

Put the second part in a form to use the Euler-Maclaurin summation formula.

\sum_{k = a}^{b} \frac{1}{(k + z)^{n + 1}} = \sum_{k = 0}^{b - a} \frac{1}{(a + k + z)^{n + 1}}

m = b - a, h = 1, F (a + k h) = \frac{1}{(a + k + z)^{n + 1}}

Define the function.

F (x) = \frac{1}{(x + z)^{n + 1}}

Find the derivatives of the function.

F^{(1)} (x) = - \frac{n + 1}{(x + z)^{n + 2}}

F^{(3)} (x) = - \frac{(n + 1) (n + 2) (n + 3)}{(x + z)^{n + 4}}

:

F^{(2 k - 1)} (x) = - Π_{j = 1}^{2 k - 1} (n + j) (x + z)^{- n - 2 k}

Calculate the integral of the function.

\int_{a}^{b} F (t) dt = - \frac{1}{n (z + b)^{n}} + \frac{1}{n (z + b)^{n}}

Substitute the three parts into the summation. The Euler-Maclaurin summation becomes:

\sum_{k = a}^{b} F (k) = - \frac{1}{n (z + b)^{n}} + \frac{1}{n (z + a)^{n}} + \frac{1}{2} (F (b) + F (a)) -

\sum_{k = 1}^{n - 1} \frac{B_{2 k}}{(2 k)!} Π_{j = 1}^{2 k - 1} (n + j) ((b + z)^{- n - 2 k} - (a + z)^{- n - 2 k}) + e_{n}

Take the limit.

\sum_{k = a}^{\infty} F (k) = \frac{1}{n (z + a)^{n}} + \frac{F (a)}{2} + \sum_{k = 1}^{n - 1} \frac{B_{2 k}}{(2 k)! (a + z)^{n + 2 k}} Π_{j = 1}^{2 k - 1} (n + j) + e_{n}

Add more terms to the summation.

\sum_{k = 0}^{\infty} F (k) = \sum_{k = 0}^{a - 1} F (k) + \frac{1}{n (z + a)^{n}} + \frac{F (a)}{2} + \sum_{k = 1}^{n - 1} \frac{B_{2 k}}{(2 k)! (a + z)^{n + 2 k}} Π_{j = 1}^{2 k - 1} (n + j) + e_{n}

Substitute the Euler-Maclaurin summation into the split series for the polygamma function. The polygamma function becomes:

ψ^{(n)} (z) = (- 1)^{n + 1} n!

[\sum_{k = 0}^{a - 1} \frac{1}{(k + z)^{n + 1}} + \frac{1}{n (a + z)^{n}} + \frac{1}{2 (a + z)^{n + 1}} + \sum_{k = 1}^{n - 1} \frac{B_{2 k}}{(2 k)! (a + z)^{n + 2 k}} Π_{j = 1}^{2 k - 1} (n + j) + e_{n}]

This is how to generate the Euler numbers.

E_{0} = 1, E_{1} = 0, E_{2} = - 1

There are special values for the Euler numbers.

E_{2 n + 1} = 0 (n = 1, 2, \dots)

This is a recursive formula to generate Euler numbers.

E_{2 n} = - 1 - (2 n)! \sum_{k = 1}^{n - 1} \frac{E_{2 k}}{(2 k)! (2 (n - k))!} (n = 1, 2, \dots)

The definition of the Riemann zeta function is:

ζ (s) = \sum_{k = 1}^{\infty} k^{- s} &rfraktur; s > 1

The definition is not the way to evaluate the Riemann zeta function except for large positive values of

s

. For small values of

s

the series converges very slowly. The Euler-Maclaurin summation formula is a way to help slowly converging series. The Euler summation expands the coefficients of the series.

The Euler-Maclaurin summation formula is:

\sum_{k = 0}^{m} F (a + k h) = \frac{1}{h} \int_{a}^{b} F (t) dt + \frac{1}{2} (F (b) + F (a)) + \sum_{k = 1}^{n - 1} \frac{h^{2 k - 1}}{(2 k)!} B_{2 k} (F^{(2 k - 1)} (b) - F^{(2 k - 1)} (a)) + e_{n}

B_{2 k}

is a Bernoulli number,

b = a + m h

, and

e_{n}

is the error.

The Euler-Maclaurin summation is an asymptotic series. The way to use the Euler-Maclaurin summation is to get the value of

a

large enough so the summation gets an accurate result before it diverges. This is the infinite series to evaluate the Riemann zeta function.

ζ (s) = \sum_{k = 1}^{\infty} k^{- s}

It can be split into two parts.

\sum_{k = 1}^{\infty} k^{- s} = \sum_{k = 1}^{a - 1} k^{- s} + \sum_{k = a}^{\infty} k^{- s}

Put it in a form to use the Euler-Maclaurin summation formula.

\sum_{k = a}^{\infty} k^{- s} = \sum_{k = 0}^{\infty} (a + k)^{- s}

m \to \infty, h = 1, F (a + k h) = (a + k h)^{- s}

Define the function.

F (x) = x^{- s}

Find the derivatives of the function.

F^{(1)} (x) = - s x^{- s - 1}

F^{(3)} (x) = - s (s + 1) (s + 2) x^{- s - 3}

:

F^{(2 k - 1)} (x) = - Π_{j = 0}^{2 k - 2} (s + j) x^{- s - (2 k - 1)}

Calculate the integral of the function.

{\int_{a}^{b} F (t) dt = \frac{t^{1 - s}}{1 - s} |}_{a}^{b}

Some of the terms in the Euler-Maclaurin summation drop out.

lim_{x \to \infty +} (x^{- s}) = 0

lim_{x \to \infty +} (\frac{x^{1 - s}}{1 - s}) = 0

lim_{x \to \infty +} (x^{- s - (2 k - 1)}) = 0 (k > 0)

Substitute the three remaining parts into the summation. The Euler-Maclaurin summation becomes:

\sum_{k = 0}^{\infty} (a + k)^{- s} = - \frac{a^{1 - s}}{1 - s} + \frac{a^{- s}}{2} + \sum_{k = 1}^{n - 1} \frac{B_{2 k}}{(2 k)!} (- (- Π_{j = 0}^{2 k - 2} (s + j) a^{- s - (2 k - 1)})) + e_{n}

Substitute the Euler-Maclaurin summation into the split series for the Riemann zeta function. The zeta function becomes:

ζ (s) = \sum_{k = 1}^{a - 1} k^{- s} + a^{1 - s} (\frac{1}{s - 1} + \frac{1}{2 a} + \sum_{k = 1}^{n - 1} \frac{B_{2 k}}{(2 k)! a^{2 k}} Π_{j = 0}^{2 k - 2} (s + j)) + e_{n}

This is the way to evaluate the Riemann zeta function for arguments with a positive real part. For arguments with negative real parts there is another formula.

ζ (s) = 2^{s} π^{s - 1} \sin (\frac{π s}{2}) Γ (1 - s) ζ (1 - s)

There are special values for the Riemann zeta function.

ζ (0) = - \frac{1}{2}

ζ (1) = \infty

ζ (- 2 n) = 0 (n = 1, 2, \dots)

ζ (1 - 2 n) = - \frac{B_{2 n}}{2 n} (n = 1, 2, \dots)

ζ (2 n) = \frac{(2 π)^{2 n}}{2 (2 n)!} | B_{2 n} | (n = 1, 2, \dots)

We need the zeta function minus one to evaluate the gamma function.

ζ (s) - 1 = \sum_{k = 2}^{a - 1} k^{- s} + a^{1 - s} (\frac{1}{s - 1} + \frac{1}{2 a} + \sum_{k = 1}^{n - 1} \frac{B_{2 k}}{(2 k)! a^{2 k}} Π_{j = 0}^{2 k - 2} (s + j)) + e_{n}

The limit is not a very good way to numerically evaluate the gamma function. Some other form would be better. A start is to invert the gamma function. That will bound the value of the function. It will go to zero rather than going to infinity.

\frac{1}{Γ (z)} = lim_{n \to \infty} \frac{Π_{k = 0}^{n} z + k}{n! n^{z}} (‖ z ‖ < \infty)

Factor

z

from the numerator and divide by the factorial in the denominator.

\frac{1}{Γ (z)} = lim_{n \to \infty} z n^{- z} Π_{k = 1}^{n} (\frac{z}{k} + 1)

Now I will do a wild guess since I know the answer.

n^{- z} = e^{- z \ln (n)}

\ln (n) = \sum_{k = 1}^{n} \frac{1}{k} - \sum_{k = 1}^{n} \frac{1}{k} + \ln (n)

Substitute for

\ln (n)

n^{- z} = e^{- z [\sum_{k = 1}^{n} \frac{1}{k} - \sum_{k = 1}^{n} \frac{1}{k} + \ln (n)]}

Convert one of the sums to a product.

n^{- z} = Π_{k = 1}^{n} e^{- \frac{z}{k}} e^{z [\sum_{k = 1}^{n} \frac{1}{k} - \ln (n)]}

Substitute for

n^{- z}

and factor out

z

\frac{1}{Γ (z)} = z lim_{n \to \infty} Π_{k = 1}^{n} e^{- \frac{z}{k}} e^{z [\sum_{k = 1}^{n} \frac{1}{k} - \ln (n)]} Π_{k = 1}^{n} (\frac{z}{k} + 1)

\frac{1}{Γ (z)} = z lim_{n \to \infty} e^{z [\sum_{k = 1}^{n} \frac{1}{k} - \ln (n)]} Π_{k = 1}^{n} [(\frac{z}{k} + 1) e^{- \frac{z}{k}}]

γ

is Euler's constant which is defined as:

γ = lim_{m \to \infty} [\sum_{k = 1}^{m} \frac{1}{k} - \ln (m)]

Substitute Euler's constant for the limit.

\frac{1}{Γ (z)} = z e^{γ z} Π_{k = 1}^{\infty} [(\frac{z}{k} + 1) e^{- \frac{z}{k}}]

Now we have an infinite product instead of a limit. This is easier to evaluate but this is not the final form yet.

The Bernoulli numbers can be generated in DERIVE. This is a simple but not very fast way to generate the Bernoulli numbers.

B_AUX(n):=-1/(n+1)+1/2-n/12-SUM(B_AUX(2*k)/(2*k)!*PRODUCT(n-j,j,0,2*k-2),k,2,n~ /2-1) B(n):=IF(n=0,1,IF(n=1,-1/2,IF(n=2,1/6,IF(MOD(n,2)=1,0,B_AUX(n))))) VECTOR(B(n),n,0,12) [1,-1/2,1/6,0,-1/30,0,1/42,0,-1/30,0,5/66,0,-691/2730]

This method suffers from fan out problems. The same values are calculated again. A better way is to store values in an array and just recall any values already calculated.

Now I will calculate Euler's constant to 64 digits.

‖ e_{n} ‖ < 10^{- 64}

‖ e_{n} ‖ < | \frac{B_{2 (n - 1)}}{2 (n - 1) a^{2 (n - 1)}} |

a = 2^{12} n = 10 ‖ e_{n} ‖ < 2.8999 10^{- 65}

γ = 0.5772156649015328606065120900824024310421593359399235988057672349

This is the DERIVE code used to generate Euler's constant. The Bernoulli numbers were stored in an array to speed up the calculation.

B_AUX(n):=-1/(n+1)+1/2-n/12-SUM(B_AUX(2*k)/(2*k)!*PRODUCT(n-j,j,0,2*k-2),k,2,n~ /2-1) B(n):=IF(n=0,1,IF(n=1,-1/2,IF(n=2,1/6,IF(MOD(n,2)=1,0,B_AUX(n))))) B_AUX0(n):=ELEMENT(VECTOR(B(2*k),k,1,12),n/2) B_AUX0(n):=ELEMENT([1/6,-1/30,1/42,-1/30,5/66,-691/2730,7/6,-3617/510,43867/79~ 8,-174611/330,854513/138,-236364091/2730],n/2) B(n):=IF(n=0,1,IF(n=1,-1/2,IF(MOD(n,2)=1,0,IF(n<26,B_AUX0(n),B_AUX(n))))) EG(a,n):=SUM(1/k,k,1,a-1)-LOG(a)+1/(2*a)+SUM(B(2*k)/(2*k*a^. (2*k)),k,1,n-1) ER(a,k):=B(2*k)/(2*k*a^(2*k)) ER(2^12,9)=2.8999*10^(-65) euler_gamma:=EG(2^12,10) euler_gamma:=0.5772156649015328606065120900824024310421593359399235988057672349

DERIVE has a very good logarithmic function. If the system you are using does not have a very good logarithmic function, then we need a method to calculate it to high precision. The Taylor series is a good way to evaluate the logarithmic function. The Taylor power series is:

f (x + h) = \sum_{k = 0}^{\infty} \frac{h^{k}}{k!} f^{(k)} (x)

We want to expand the logarithmic function with the Taylor power series.

f (x) = \ln (x)

The derivatives of the logarithmic function are:

f^{1} (x) = \frac{1}{x}

f^{2} (x) = - \frac{1}{x^{2}}

f^{3} (x) = \frac{2}{x^{3}}

f^{4} (x) = - \frac{6}{x^{4}}

:

f^{n} (x) = - \frac{(n - 1)!}{(- x)^{n}}

Expand the logarithmic function about one.

x = 1

f (1 + h) = 0 + \sum_{k = 1}^{\infty} - \frac{h^{k}}{k!} \frac{(k - 1)!}{(- 1)^{k}}

\ln (1 + h) = \sum_{k = 1}^{\infty} - \frac{(- h)^{k}}{k}

This series converges very slowly except for small values of

h

. These idenities can be used to find the logarithm of values when

h

is not small.

\ln (a b) = \ln (a) + \ln (b) \ln (\frac{a}{b}) = \ln (a) - \ln (b) \ln (a^{b}) = b \ln (a)

Since the numbers in the computer are binary the constant which would be the most useful is

\ln (2)

\ln (2) = \sum_{k = 1}^{\infty} - \frac{(- 1)^{k}}{k}

This series converges very slow for

\ln (2)

. The Euler-Maclaurin summation will speed the convergence. The series has a power of

k

in it. That will be hard to integrate. I will have to find a simplier series which contains no powers.

\ln (2) = 1 - \frac{1}{2} + \frac{1}{3} - \frac{1}{4} + \dots

A wild guess will give this new series.

\ln (2) = \sum_{k = 1}^{\infty} \frac{1}{2 k - 1} - \frac{1}{2 k}

Now put the series in a form which can be evaluated using the Euler-Maclaurin summation formula.

F (x) = \frac{1}{2 x - 1} - \frac{1}{2 x}

Substitute this function in the

\ln (2)

series.

\ln (2) = \sum_{k = 1}^{a - 1} F (k) + \sum_{k = a}^{\infty} F (k)

The second term is the right form.

\sum_{k = a}^{\infty} F (k) = \int_{a}^{\infty} F (t) dt + \frac{F (a)}{2} - \sum_{k = 1}^{\infty} \frac{B_{2 k}}{(2 k)!} F^{(2 k - 1)} (a)

Evaluate the integral.

\int_{a}^{\infty} F (t) dt = - \frac{\ln (1 - \frac{1}{2 a})}{2}

Find a general formula for the derivatives.

F^{1} (x) = - \frac{2}{(2 x - 1)^{2}} + \frac{2}{(2 x)^{2}}

F^{2} (x) = \frac{8}{(2 x - 1)^{3}} - \frac{8}{(2 x)^{3}}

:

F^{n} (x) = (- 2)^{n} n! (\frac{1}{(2 x - 1)^{n + 1}} - \frac{1}{(2 x)^{n + 1}})

Substitute these three items into the Euler-Maclaurin summation formula.

\sum_{k = a}^{\infty} F (k) = - \frac{\ln (1 - \frac{1}{2 a})}{2} + \frac{F (a)}{2} - \sum_{k = 1}^{\infty} \frac{B_{2 k}}{(2 k)!} (- 2)^{2 k - 1} (2 k - 1)! (\frac{1}{(2 a - 1)^{2 k}} - \frac{1}{(2 a)^{2 k}})

Divide the factorial.

\sum_{k = a}^{\infty} F (k) = - \frac{\ln (1 - \frac{1}{2 a})}{2} + \frac{F (a)}{2} + \sum_{k = 1}^{\infty} \frac{B_{2 k} 2^{2 k - 2}}{k} (\frac{1}{(2 a - 1)^{2 k}} - \frac{1}{(2 a)^{2 k}})

Replace the second term in the series for

\ln (2)

\ln (2) = \sum_{k = 1}^{a - 1} F (k) - \frac{\ln (1 - \frac{1}{2 a})}{2} + \frac{F (a)}{2} + \sum_{k = 1}^{\infty} \frac{B_{2 k} 2^{2 k - 2}}{k} (\frac{1}{(2 a - 1)^{2 k}} - \frac{1}{(2 a)^{2 k}})

This is the DERIVE code to calculate

\ln (2^{12})

which is required to calculate Euler's constant. DERIVE can calculate the logarithmic function very fast. This is just an example of how to calculate the logarithmic function numerically if it was not built in.

LOG(x) DIF(LOG(x),x,k) VECTOR(DIF(LOG(x),x,k),k,1,8) [1/x,-1/x^2,2/x^3,-6/x^4,24/x^5,-120/x^6,720/x^7,-5040/x^8] [1,-1,2,-6,24,-120,720,-5040] VECTOR((-1)^(k-1)*(k-1)!,k,1,8) [1,-1,2,-6,24,-120,720,-5040] LOGA(n,x):=-SUM((1-x)^k/k,k,1,n) [1,-1/2,1/3,-1/4,1/5,-1/6,1/7,-1/8] F(k):=1/(2*k-1)-1/(2*k) VECTOR(DIF(F(x),x,n),n,1,8) [1/(2*x^2)-2/(2*x-1)^2,8/(2*x-1)^3-1/x^3,3/x^4-48/(2*x-1)^4,384/(2*x-1)^5-12/x~ ^5,60/x^6-3840/(2*x-1)^6,46080/(2*x-1)^7-360/x^7,2520/x^8-645120/(2*x-1)^8,103~ 21920/(2*x-1)^9-20160/x^9] VECTOR((-2)^n*n!*(-1/(2*x)^(n+1)+1/(2*x-1)^(n+1)),n,1,8) [1/(2*x^2)-2/(2*x-1)^2,8/(2*x-1)^3-1/x^3,3/x^4-48/(2*x-1)^4,384/(2*x-1)^5-12/x~ ^5,60/x^6-3840/(2*x-1)^6,46080/(2*x-1)^7-360/x^7,2520/x^8-645120/(2*x-1)^8,103~ 21920/(2*x-1)^9-20160/x^9] INT(F(t),t,a,b)=-LN(2*a-1)/2+LN(a)/2+LN(2*b-1)/2-LN(b)/2 INT(F(t),t,a,inf)=(LN(a)-LN(a-1/2))/2 B_AUX(n):=-1/(n+1)+1/2-n/12-SUM(B_AUX(2*k)/(2*k)!*PRODUCT(n-j,j,0,2*k-2),k,2,n~ /2-1) B(n):=IF(n=0,1,IF(n=1,-1/2,IF(n=2,1/6,IF(MOD(n,2)=1,0,B_AUX(n))))) B_AUX0(n):=ELEMENT(VECTOR(B(2*k),k,1,12),n/2) B_AUX0(n):=ELEMENT([1/6,-1/30,1/42,-1/30,5/66,-691/2730,7/6,-3617/510,43867/79~ 8,-174611/330,854513/138,-236364091/2730],n/2) B(n):=IF(n=0,1,IF(n=1,-1/2,IF(MOD(n,2)=1,0,IF(n<26,B_AUX0(n),B_AUX(n))))) LOG2(a,n):=SUM(F(k),k,1,a-1)-LN(1-1/(2*a))/2+F(a)/2+SUM(B(2*k)*2^(2*k-2)/k*(1/~ (2*a-1)^(2*k)-1/(2*a)^(2*k)),k,1,n) ER(a,k):=B(2*k)*2^(2*k-2)/k*(1/(2*a-1)^(2*k)-1/(2*a)^(2*k)) ER(2^10,11)=9.02239*10^(-67) LOG2(2^10,11)-LN(2)=1.20366*10^(-71) LOG2(2^10,11)=0.69314718055994530941723212145817656807550013436025525412068000~ 94 LOGA(n,1-1/(2*a))=-SUM((1/(2*a))^k/k,k,1,n) LNA(a,n):=-SUM((1/(2*a))^k/k,k,1,n) ERA(a,k):=(1/(2*a))^k/k ERA(2^10,24)=1.40562*10^(-81) LNA(2^10,24)-LN(1-1/(2*a)) LNA(2^10,24)-LN(1-1/(2*2^10))=6.59196*10^(-85) LNA(2^10,24)=-4.88400498108874464984963559896910983460775989412689222147149656~ 3*10^(-4) LOG(2^12)=12*LN(2) LOG(2^12)=8.317766166719343713006785457498118816906001612323063049448160113

\ln (2^{12}) = 12 \ln (2)

\ln (2) = 0.6931471805599453094172321214581765680755001343602552541206800094

\ln (2^{12}) = 8.317766166719343713006785457498118816906001612323063049448160113

The estimated error in calculating

\ln (2)

9.02239 10^{- 67}

. The actual error is

1.20366 10^{- 71}

. The calculation used 2048 terms in the Taylor series and 11 terms in the Euler-Maclaurin summation. The actual power series had 1024 terms since

F (k)

contains two Taylor series terms.

Now that Euler's constant is calculated we can get back to calculating the gamma function.

\frac{1}{Γ (z)} = z e^{γ z} Π_{k = 1}^{\infty} [(\frac{z}{k} + 1) e^{- \frac{z}{k}}]

The range of

z

can be limited to

1 \leq z < 2

by using the recurrence formula.

z > 2

then

Γ (z) = (z - 1) Γ (z - 1)

z < 1

then

Γ (z) = \frac{Γ (z + 1)}{z}

1 \leq z < 2

then

0 > Γ (z) \leq 1

. Since the range of values the gamma function is limited it can be inverted.

Γ (z) = z^{- 1} e^{- γ z} Π_{k = 1}^{\infty} [\frac{e^{\frac{z}{k}}}{(\frac{z}{k} + 1)}]

Take the logramithm to convert the infinite product into an infinite series.

\ln (Γ (z)) = - \ln (z) - γ z + \sum_{k = 1}^{\infty} [\frac{z}{k} - \ln (\frac{z}{k} + 1)]

The digamma function is defined as:

ψ (z) = \frac{d \ln (Γ (z))}{dz}

Substituting the logramithm of the gamma function gives a series for the digamma function.

ψ (z) = - \frac{1}{z} - γ + z \sum_{k = 1}^{\infty} \frac{1}{k (k + z)} (z \neq 0, - 1, - 2, - 3, \dots)

The polygamma function is defined as:

ψ^{(n)} (z) = \frac{d^{n}}{{dz}^{n}} ψ (z)

Substituting the series for the digamma function gives the series for the polygamma function.

ψ^{(n)} (z) = (- 1)^{n + 1} n! \sum_{k = 0}^{\infty} \frac{1}{(k + z)^{n + 1}} (z \neq 0, - 1, - 2, - 3, \dots)

The Taylor power series is:

f (x + h) = \sum_{k = 0}^{\infty} \frac{h^{k}}{k!} f^{(k)} (x)

The logramithm of the gamma function can be evaluated using the Taylor series.

f (x) = \ln (Γ (x))

The derivative is the digamma function.

f^{(1)} (x) = ψ (x)

The higher order derivatives are the polygamma functions.

f^{(2)} (x) = ψ^{(1)} (x)

:

f^{(n)} (x) = ψ^{(n - 1)} (x)

Expand the Taylor series about one.

x = 1

f (1 + h) = f (1) + f^{(1)} (1) h + \sum_{k = 2}^{\infty} \frac{h^{k}}{k!} f^{(k)} (1)

Substitute in the Taylor series the logarithm of the gamma function.

f (1 + h) = \ln (Γ (1)) + ψ (1) h + \sum_{k = 2}^{\infty} \frac{h^{k}}{k!} ψ^{(k - 1)} (1)

Simplify the series with known values.

f (1 + h) = - γ h + \sum_{k = 2}^{\infty} \frac{h^{k}}{k!} (- 1)^{k} (k - 1)! \sum_{j = 0}^{\infty} \frac{1}{(j + 1)^{k}}

Divide by the factorial.

\ln (Γ (1 + h)) = - γ h + \sum_{k = 2}^{\infty} \frac{(- h)^{k}}{k} \sum_{j = 1}^{\infty} \frac{1}{j^{k}}

The definition of the Riemann zeta function is:

ζ (s) = \sum_{k = 1}^{\infty} k^{- s} &rfraktur; s > 1

Substitute the Riemann zeta function in the series for the logarithm of the gamma function.

\ln (Γ (1 + h)) = - γ h + \sum_{k = 2}^{\infty} \frac{(- h)^{k}}{k} ζ (k)

Subtract one from the argument of the Riemann zeta function so the series converges faster.

\ln (Γ (1 + h)) = - γ h + \sum_{k = 2}^{\infty} \frac{(- h)^{k}}{k} (ζ (k) - 1) + \sum_{k = 2}^{\infty} \frac{(- h)^{k}}{k}

The Taylor power series for the logarithmic function is:

\ln (1 + h) = \sum_{k = 1}^{\infty} - \frac{(- h)^{k}}{k}

Substitute the logarithm for its series.

\ln (Γ (1 + h)) = - γ h + \sum_{k = 2}^{\infty} \frac{(- h)^{k}}{k} (ζ (k) - 1) - \ln (1 + h) + h

Change variables and simplify.

\ln (Γ (1 + z)) = - \ln (1 + z) + z (1 - γ) + \sum_{n = 2}^{\infty} \frac{(- z)^{n} (ζ (n) - 1)}{n}

This is now a power series to evaluate the gamma function.

The DERIVE code to generate the zeta function is:

B_AUX(n):=-1/(n+1)+1/2-n/12-SUM(B_AUX(2*k)/(2*k)!*PRODUCT(n-j,j,0,2*k-2),k,2,n~ /2-1) B(n):=IF(n=0,1,IF(n=1,-1/2,IF(n=2,1/6,IF(MOD(n,2)=1,0,B_AUX(n))))) B_AUX0(n):=ELEMENT(VECTOR(B(2*k),k,1,12),n/2) B_AUX0(n):=ELEMENT([1/6,-1/30,1/42,-1/30,5/66,-691/2730,7/6,-3617/510,43867/79~ 8,-174611/330,854513/138,-236364091/2730],n/2) B(n):=IF(n=0,1,IF(n=1,-1/2,IF(MOD(n,2)=1,0,IF(n<26,B_AUX0(n),B_AUX(n))))) F(k):=k^(-s) F(x) DIF(F(x),x) DIF(F(x),x,k) VECTOR(DIF(F(x),x,k),k,1,4) [-s*x^(-s-1),s*x^(-s-2)*(s+1),-s*x^(-s-3)*(s+1)*(s+2),s*x^(-s-4)*(s+1)*(s+2)*(~ s+3)] (-1)^n*PRODUCT(s+j,j,0,n-1)*x^(-s-n) VECTOR((-1)^n*PRODUCT(s+j,j,0,n-1)*x^(-s-n),n,1,4) [-s*x^(-s-1),s*x^(-s-2)*(s+1),-s*x^(-s-3)*(s+1)*(s+2),s*x^(-s-4)*(s+1)*(s+2)*(~ s+3)] INT(F(t),t,a,inf)=a^(1-s)/(s-1) Z1(s,a,m):=SUM(k^(-s),k,2,a-1)+a^(1-s)*(1/(s-1)+1/(2*a)+SUM(B(2*k)/((2*k)!*a^(~ 2*k))*PRODUCT(s+j,j,0,2*k-2),k,1,m)) ER(s,a,k):=a^(1-s)*(B(2*k)/((2*k)!*a^(2*k))*PRODUCT(s+j,j,0,2*k-2)) ER(2,2^10,11)=3.58872*10^(-66) Z1(2,2^10,11)+1-ZETA(2)=4.78534*10^(-71) Z1(2,2^10,11)=0.64493406684822643647241516664602518921894990120679843773555822~ 93

The worst case error is when

z = 2

. That is the lowest value required to calculate the gamma function. The predicted error is

3.58872 10^{- 66}

. The actual error is

4.78534 10^{- 71}

. The calculation used 1024 terms in the series and 11 terms in the Euler-Maclaurin summation. For larger values of

z

fewer terms in the series are required.

We now have a way to generate the zeta function, so let us get back to the gamma function. The next value we need to calculate is

\ln (1 + z)

for

- \frac{1}{2} < z \leq \frac{1}{2}

. The Taylor power series is:

\ln (1 + x) = - \sum_{k = 1}^{\infty} \frac{(- x)^{k}}{k}

The worst case convergence occurs when

x = - \frac{1}{2}

. The estimated error for 207 terms is

2.34866 10^{- 65}

. The actual error is

2.32629 10^{- 65}

. That is too many terms. The logarithmic idenities can easily reduce the size of

x

and thus reduce the number of terms required. If

x < 2^{- \frac{1}{4}}

then

\ln (x) = \ln (x 2^{\frac{1}{2}}) - \ln (2^{\frac{1}{2}})

. If

x > 2^{\frac{1}{4}}

then

\ln (x) = \ln (x 2^{- \frac{1}{2}}) + \ln (2^{\frac{1}{2}})

. The worst case convergence occurs when

x = 2^{\frac{1}{4}}

. The estimated error for 87 terms is

- 1.42569 10^{- 65}

. The actual error is

2.24664 10^{- 66}

. That saves time calculating. The two new constants needed now are

\ln (2^{\frac{1}{2}})

and

2^{\frac{1}{2}}

\ln (2^{\frac{1}{2}}) = \frac{\ln (2)}{2}

\ln (2^{\frac{1}{2}}) = 0.3465735902799726547086160607290882840377500671801276270603400047

2^{\frac{1}{2}} = e^{\frac{\ln (2)}{2}}

e^{x} = \sum_{k = 0}^{\infty} \frac{x^{k}}{k!}

The estimated error for 40 terms is

4.78821 10^{- 66}

. The actual error is

- 3.40707 10^{- 66}

. The result is:

e^{\frac{\ln (2)}{2}} = 1.414213562373095048801688724209698078569671875376948073176679737

The logarithmic idenities can be used again to reduce the size of

x

. If

x < 2^{- \frac{1}{8}}

then

\ln (x) = \ln (x 2^{\frac{1}{4}}) - \ln (2^{\frac{1}{4}})

. If

x > 2^{\frac{1}{8}}

then

\ln (x) = \ln (x 2^{- \frac{1}{4}}) + \ln (2^{\frac{1}{4}})

. The worst case convergence occurs when

x = 2^{\frac{1}{8}}

. The estimated error for 61 terms is

- 3.73665 10^{- 66}

. The actual error is

3.05538 10^{- 67}

. That saves some more time. The two new constants needed now are

\ln (2^{\frac{1}{4}})

and

2^{\frac{1}{4}}

\ln (2^{\frac{1}{4}}) = \frac{\ln (2)}{4}

\ln (2^{\frac{1}{4}}) = 0.1732867951399863273543080303645441420188750335900638135301700023

2^{\frac{1}{4}} = e^{\frac{\ln (2)}{4}}

The estimated error for 35 terms is

2.20069 10^{- 67}

. The actual error is

- 1.06429 10^{- 69}

. The result is:

e^{\frac{\ln (2)}{4}} = 1.189207115002721066717499970560475915292972092463817413019002224

Finally we can calculate the gamma function using the power series.

\ln (Γ (1 + z)) = - \ln (1 + z) + z (1 - γ) + \sum_{n = 2}^{\infty} \frac{(- z)^{n} (ζ (n) - 1)}{n}

The worst case value is

z = - \frac{1}{2}

. The estimated error with

z = - \frac{1}{2}

and 104 terms is

2.33737 10^{- 65}

. The actual error is

- 1.07061 10^{- 65}

. The time consuming part of the calculations is the calculation of the zeta function. The 104 values of the zeta function should be saved so the next value of the gamma function will take less time. This is the DERIVE code used to calculate the zeta function. 199 values of the zeta function are saved in an array.

Z1(k):=ELEMENT([1915102941003759390681506705367659/296945539001034692893587104~ 1107579,404980503851378643563268289083258/2004289373531106300586436804239795,3~

There are several pages of numbers left out. There are a total of 199 rational numbers.

00000000000000000000000000000000000000000000000/463119544355440997411193475012~ 723059245072018560683104122996917*10^(-61)],k-1) euler_gamma:=0.5772156649015328606065120900824024310421593359399235988057672349 eg1:=1-euler_gamma eg1:=2702377911746907625973774294833889/6391859128644204715319046100816835 LG(z,n):=-LN(1+z)+z*eg1+SUM((-z)^k*Z1(k)/k,k,2,n) EG(-1/2,104)=2.33737*10^(-65) LG(-1/2,104)-LN(GAMMA(1-1/2))=-1.07061*10^(-65)

The algorithms are programmed in LISP. This is how to generate the Bernoulli numbers.

B_{0} = 1, B_{1} = - 1 / 2, B_{2} = 1 / 6

B_{2 n + 1} = 0 (n = 1, 2, \dots)

B_{2 n} = - \frac{1}{2 n + 1} + 1 / 2 - \frac{n}{6} - \sum_{k = 2}^{n - 1} \frac{B_{2 k}}{(2 k)!} Π_{j = 0}^{2 k - 2} (2 n - j) (n = 2, 3, \dots)

This is how I calculate Bernoulii numbers in muLISP. The output is compatible with DERIVE.

$ (precision 213) ;1024 digits $ (setq *print-point* nil) $ (setq *exact-mode* T) $ (defun Bernoulli (n l total term k j i) $ ((zerop n) 1) $ ((eql n 1) -1/2) $ ((eql n 2) 1/6) $ ((oddp n) 0) $ (setq l (* 2 (add1 (length B)))) $ (loop $ ((<= n l) (nth (- (/ n 2) 2) B) ) $ (setq k 2 j (incq l 2) i -1) $ (setq term (/ j k)) $ (setq total (- 1/2 (/ 1 (+ 1 l)) (/ l 12))) $ (loop $ ((eql k (- l 2))) $ (setq term (* term (decq j))) $ (setq term (* term (decq j))) $ (setq term (/ term (incq k))) $ (setq term (/ term (incq k))) $ (decq total (* (nth (incq i) B) term)) ) $ (setq B (append B (list total))) ) ) $ (open-output-file 'bernoull) $ (setq B nil n 2) $ (setq *auto-newline* nil) $ (princ ``['') $ (time T) $ (loop $ (prin1 (Bernoulli n)) $ ((eql n 512)) $ (princ ``, '')(terpri) $ (incq n 2) ) $ (princ ``]'') $ (terpri) $ (print (time)) $ (system)

This is how I calculate one minus Euler's constant in LISP.

γ = \sum_{k = 1}^{a - 1} \frac{1}{k} - \ln (a) + \frac{1}{2 a} + \sum_{k = 1}^{n - 1} \frac{B_{2 k}}{2 k a^{2 k}} + e_{n}

$ (load 'irratnal) $ (precision 213) $ (setq *print-point* nil) $ (setq *exact-mode* T) $ (defun Bernoulli (n l total term k j i) $ ((zerop n) 1) $ ((eql n 1) -1/2) $ ((eql n 2) 1/6) $ ((oddp n) 0) $ (setq l (* 2 (add1 (length *Bernoulli*)))) $ (loop $ ((<= n l) (nth (- (/ n 2) 2) *Bernoulli*) ) $ (setq total 1/2 k 2 j (incq l 2) i -1) $ (decq total (+ (/ 1 (+ 1 l)) (/ (setq term (/ j k)) 6) ) ) $ (loop $ ((eql k (- l 2))) $ (setq term (* term (decq j))) $ (setq term (* term (decq j))) $ (setq term (/ term (incq k))) $ (setq term (/ term (incq k))) $ (decq total (* (nth (incq i) *Bernoulli*) term)) ) $ (setq *Bernoulli* (append *Bernoulli* (list total))) ) ) $ (defun Euler-gamma (n m k bit total s) $ ((null *Euler-gamma*) $ (setq k 1 n (shift 1 13) *Euler-gamma* 0) $ (loop $ (incq *Euler-gamma* (/ 1 k)) $ ((>= (incq k) n)) ) $ (incq *Euler-gamma* (/ 1/2 n)) $ (print-time) $ (setq *exact-mode* nil) $ (decq *Euler-gamma* (log n)) $ (setq *exact-mode* T) $ (print-time) $ (setq k 0 bit 0 m 404 s (shift 1 3416) total 0) $ (loop $ ((> (incq k 2) m)) $ (incq total $ (round (* s (/ (Bernoulli k) (shift k (incq bit 26))))) ) ) $ (print-time) $ (setq total (/ total s)) $ (setq *exact-mode* nil) $ (incq *Euler-gamma* total) $ (setq *exact-mode* T) $ (print-time) $ *Euler-gamma* ) $ *Euler-gamma* )

$ (defun print-time (temp) $ (setq temp (time)) $ (print (list temp (- temp start-time))) $ (setq start-time temp) ) $ (open-output-file 'eulerg) $ (setq *Bernoulli* nil *Euler-gamma* nil start-time 0) $ (setq *auto-newline* nil) $ (time T) $ (print (setq temp (- 1 (Euler-gamma)))) $ (print-time) $ (setq *exact-mode* nil) $ (setq *print-point* 1024) $ (print temp) $ (setq *exact-mode* T) $ (setq *print-point* nil) $ (print-time) $ (system)

This is gamma function power series again.

\ln (Γ (1 + z)) = - \ln (1 + z) + z (1 - γ) + \sum_{n = 2}^{\infty} \frac{(- z)^{n} (ζ (n) - 1)}{n}

The worst case value is

z = - \frac{1}{2}

. The estimated error with

z = - \frac{1}{2}

and 1697 terms is

1.18716 10^{- 1025}

. The time consuming part of the calculations is the calculation of the zeta function. The 1697 values of the zeta function would take several hours on my computer. There must be a better way to calculate the gamma function to very high precision.

Let's try something. Using the duplication formula with

z = \frac{1}{2} + δ

Γ (1 + 2 δ) = π^{- \frac{1}{2}} 2^{2 δ} Γ (\frac{1}{2} + δ) Γ (1 + δ)

Γ (\frac{1}{2} + δ) = \frac{Γ (1 + 2 δ)}{Γ (1 + δ)} π^{\frac{1}{2}} 2^{- 2 δ}

This formula will reduce the time to compute the gamma function if

\frac{1}{2} - δ > 2 δ

. The cross over point is when

\frac{1}{2} - δ = 2 δ

\frac{1}{2} = 3 δ δ = \frac{1}{6}

Γ (\frac{1}{2} - δ) = \frac{Γ (1 - 2 δ)}{Γ (1 - δ)} π^{\frac{1}{2}} 2^{2 δ}

Γ (\frac{3}{2} - δ) = (\frac{1}{2} - δ) \frac{Γ (1 - 2 δ)}{Γ (1 - δ)} π^{\frac{1}{2}} 2^{2 δ}

The new worst case in the power series is now

z = - 2 δ = - \frac{2}{6} = - \frac{1}{3}

. The estimated error with

z = - \frac{1}{3}

and 1313 terms is

1.47619 10^{- 1025}

. 1313 values of the zeta function need to be calculated. That will still take too long. Even though it takes several hours to run here is the muLISP program that generates 1696 values of the zeta function to 1024 decimal digits.

ζ (s) - 1 = \sum_{k = 2}^{a - 1} k^{- s} + a^{1 - s} (\frac{1}{s - 1} + \frac{1}{2 a} + \sum_{k = 1}^{n - 1} \frac{B_{2 k}}{(2 k)! a^{2 k}} Π_{j = 0}^{2 k - 2} (s + j)) + e_{n}

$ ; This is how to generate the zeta function minus one to 1024 digits. $ $ (load 'irratnal) $ (precision 216) $ (setq *print-point* nil) $ (setq *exact-mode* T)

$ (defun Bernoulli (n l total term k j i) $ ((zerop n) 1) $ ((eql n 1) -1/2) $ ((eql n 2) 1/6) $ ((oddp n) 0) $ (setq l (* 2 (add1 (length *Bernoulli*)))) $ (loop $ ((<= n l) (nth (- (/ n 2) 2) *Bernoulli*) ) $ (setq total 1/2 k 2 j (incq l 2) i -1) $ (decq total (+ (/ 1 (+ 1 l)) (/ (setq term (/ j k)) 6) ) ) $ (loop $ ((eql k (- l 2))) $ (setq term (* term (decq j))) $ (setq term (* term (decq j))) $ (setq term (/ term (incq k))) $ (setq term (/ term (incq k))) $ (decq total (* (nth (incq i) *Bernoulli*) term)) $ ) $ (setq *Bernoulli* (append *Bernoulli* (list total))) $ ) $ ) $ $ (defun Zeta1 (s a) $ (precision 4) $ (setq *exact-mode* nil) $ (setq a (ceiling (expt 2 (/ 3416 s)))) $ (precision 216) $ (setq *exact-mode* T) $ ((< a 256) (Zeta3 s a)) $ ((> s 300) (Zeta2 s 256)) $ ((> s 64) (Zeta2 s 512)) $ ((> s 1) (Zeta2 s 1024)) $ )

$ (defun Zeta2 (s a k total r temp term j test) $ (print (list s a)) $ (setq k 2 total 0 r (shift 1 3416)) ; r is used to set the precision. $ (loop $ (incq total (round (/ r (expt k s)))) $ ((>= (incq k) a)) $ ) $ (incq total (round (/ r 2 (expt a s)))) $ (incq total (round (/ r (expt a (sub1 s)) (sub1 s) ))) $ (print-time) $ (setq term (/ r (expt a (sub1 s)))) $ (setq j s) $ (setq k 0 temp 0 term (* term j) a (* a a)) $ (loop $ (setq term (/ term (incq k))) $ (setq term (/ term (incq k))) ; 1/k! $ ((> k 2048)) $ (setq term (/ term a)) ; 1/a^(2k) $ (incq temp (setq test (round (* (Bernoulli k) term))) ) $ ((< (abs test) 16)) $ (setq term (* term (incq j))) $ (setq term (* term (incq j))) ; s + j $ ) $ (print (append (list k) (print-time-aux))) $ (incq total temp) $ ) $ $ (defun Zeta3 (s a k total r) $ (print (list s a)) $ (setq k 2 total 0 r (shift 1 3416)) ; r is used to set the precision. $ (loop $ (incq total (round (/ r (expt k s)))) $ ((>= (incq k) a)) $ ) $ (incq total (round (/ r 2 (expt a s)))) $ (incq total (round (/ r (expt a (sub1 s)) (sub1 s) ))) $ ) $ $ (defun print-time-aux (temp) $ (setq temp (time)) $ (setq temp (list temp (- temp start-time))) $ (setq start-time (car temp)) $ temp $ ) $ $ (defun print-time () $ (print (print-time-aux)) $ )

$ (open-output-file 'zeta) $ (setq *Bernoulli* nil start-time 0) $ (setq *auto-newline* nil) $ ; (setq *gcstat* T) $ (time T) $ (precision 216) $ (setq *exact-mode* T) $ (setq s 1697) $ (loop $ (print (Zeta1 s a)) $ (print-time) $ ((< (decq s) 2)) $ ) $ (system)

To get results in less than a hour, the following muLISP program only calculates 257 values of the gamma function to 154 digits.

\ln (Γ (1 + z)) = - \ln (1 + z) + z (1 - γ) + \sum_{n = 2}^{\infty} \frac{(- z)^{n} (ζ (n) - 1)}{n}

$ ; This is how to generate the gamma function to 154 digits. $ $ (load 'irratnal) $ $ (defun Bernoulli (n l total term k j i) $ ((zerop n) 1) $ ((eql n 1) -1/2) $ ((eql n 2) 1/6) $ ((oddp n) 0) $ (setq l (* 2 (add1 (length *Bernoulli*)))) $ (loop $ ((<= n l) (nth (- (/ n 2) 2) *Bernoulli*) ) $ (setq total 1/2 k 2 j (incq l 2) i -1) $ (decq total (+ (/ 1 (+ 1 l)) (/ (setq term (/ j k)) 6) ) ) $ (loop $ ((eql k (- l 2))) $ (setq term (* term (decq j))) $ (setq term (* term (decq j))) $ (setq term (/ term (incq k))) $ (setq term (/ term (incq k))) $ (decq total (* (nth (incq i) *Bernoulli*) term)) $ ) $ (setq *Bernoulli* (append *Bernoulli* (list total))) $ ) $ )

$ (defun Euler-gamma (a k total term r test) $ ((null *Euler-gamma*) $ (setq k 1 a (shift 1 13) *Euler-gamma* 0) $ (setq r (shift 1 (* 16 (precision)))) $ ; r is used to set the precision. $ (loop $ (incq *Euler-gamma* (/ 1 k)) $ ((>= (incq k) a)) $ ) $ (incq *Euler-gamma* (/ 1/2 a)) $ (setq *exact-mode* nil) $ (setq total (ln a)) $ (setq *exact-mode* T) $ (decq *Euler-gamma* total) $ (print-time) $ (setq k 2 total 0 term r a (* a a)) $ (loop $ ((> k 2048)) $ (setq term (/ term a)) ; 1/a^(2k) $ (incq total (setq test (round (* (Bernoulli k) (/ term k)) ))) $ ((< (abs test) 16)) $ (incq k 2) $ ) $ (print (append (list k) (print-time-aux))) $ (incq *Euler-gamma* (/ total r)) $ ) $ *Euler-gamma* $ ) $ $ (defun Zeta1 (s a temp) $ (setq temp (precision 4)) $ (setq *exact-mode* nil) $ (setq a (ceiling (expt 2 (/ (* 16 temp) s)))) $ (precision temp) $ (setq *exact-mode* T) $ ((< a 256) (Zeta3 s a)) $ ((> s 300) (Zeta2 s 256)) $ ((> s 64) (Zeta2 s 512)) $ ((> s 1) (Zeta2 s 1024)) $ )

$ (defun Zeta2 (s a k total r temp term j test) $ (print (list s a)) $ (setq k 2 total 0) $ (setq r (shift 1 (* 16 (precision)))) $ ; r is used to set the precision. $ (loop $ (incq total (round (/ r (expt k s)))) $ ((>= (incq k) a)) $ ) $ (incq total (round (/ r 2 (expt a s)))) $ (incq total (round (/ r (expt a (sub1 s)) (sub1 s) ))) $ (print-time) $ (setq term (/ r (expt a (sub1 s)))) $ (setq j s) $ (setq k 0 temp 0 term (* term j) a (* a a)) $ (loop $ (setq term (/ term (incq k))) $ (setq term (/ term (incq k))) ; 1/k! $ ((> k 2048)) $ (setq term (/ term a)) ; 1/a^(2k) $ (incq temp (setq test (round (* (Bernoulli k) term))) ) $ ((< (abs test) 16)) $ (setq term (* term (incq j))) $ (setq term (* term (incq j))) ; s + j $ ) $ (print (append (list k) (print-time-aux))) $ (incq total temp) $ ) $ $ (defun Zeta3 (s a k total r) $ (print (list s a)) $ (setq k 2 total 0) $ (setq r (shift 1 (* 16 (precision)))) $ ; r is used to set the precision. $ (loop $ (incq total (round (/ r (expt k s)))) $ ((>= (incq k) a)) $ ) $ (incq total (round (/ r 2 (expt a s)))) $ (incq total (round (/ r (expt a (sub1 s)) (sub1 s) ))) $ )

$ (defun lngamma (z n total term a test) $ (setq n 2 total 0 a (- 0 z) term (* z z)) $ (loop $ ((> n 1697)) $ (incq total (setq test (* (nth (- n 2) *Zeta1*) (/ term n) ))) $ ((< (abs test) 16)) $ (incq n) $ (setq term (* term a)) $ ) $ (print (append (list n) (print-time-aux))) $ (setq total (/ total (shift 1 (* 16 (precision))))) $ (setq *exact-mode* nil) $ (setq temp (ln (add1 z))) $ (setq *exact-mode* T) $ (decq total temp) $ (incq total (* z EG1)) $ ) $ $ (defun gamma (z temp) $ (setq temp (lngamma (sub1 z))) $ (setq *exact-mode* nil) $ (setq temp (exp temp)) $ (setq *exact-mode* T) $ temp $ ) $ $ (defun print-time-aux (temp) $ (setq temp (time)) $ (setq temp (list temp (- temp start-time))) $ (setq start-time (car temp)) $ temp $ ) $ $ (defun print-time () $ (print (print-time-aux)) $ )

$ (open-output-file 'gamma) $ (setq *Bernoulli* nil *Euler-gamma* nil start-time 0) $ (setq *auto-newline* nil) $ ; (setq *gcstat* T) $ (time T) $ (precision 32) $ (setq temp (precision 4)) $ (print (setq pp (ceiling (* 16 temp (/ (ln 2) (ln 10)) )))) $ (precision temp) $ (setq *print-point* nil) $ (setq *exact-mode* T) $ (setq EG1 (- 1 (Euler-gamma))) $ (print-time) $ (setq *print-point* pp) $ (print EG1) $ (setq *print-point* nil) $ (setq n 2 term 1/4 *Zeta1* nil) $ (loop $ ((> n 1697)) $ (setq test (* (setq temp (zeta1 n)) (/ term n) )) $ (setq *Zeta1* (append *Zeta1* (list temp))) $ ((< (abs test) 16)) $ (incq n) $ (setq term (/ term 2)) $ ) $ (print (append (list n) (print-time-aux))) $ (setq z 1/2) $ (loop $ ((> z 3/2)) $ (print z) $ (print (setq temp (gamma z))) $ (print-time) $ (setq *print-point* pp) $ (print temp) $ (setq *print-point* nil) $ (incq z 1/256) $ ) $ (system)

A definition of the gamma function is:

\frac{1}{Γ (z)} = z e^{γ z} Π_{k = 1}^{\infty} [(\frac{z}{k} + 1) e^{- \frac{z}{k}}]

The range of

z

can be limited to

1 \leq z < 2

by using the recurrence formula.

z > 2

then

Γ (z) = (z - 1) Γ (z - 1)

z < 1

then

Γ (z) = \frac{Γ (z + 1)}{z}

1 \leq z < 2

then

0 > Γ (z) \leq 1

. Since the range of values the gamma function is limited it can be inverted.

Γ (z) = z^{- 1} e^{- γ z} Π_{k = 1}^{\infty} [\frac{e^{\frac{z}{k}}}{(\frac{z}{k} + 1)}]

Take the logramithm to convert the infinite product into an infinite series.

\ln (Γ (z)) = - \ln (z) - γ z + \sum_{k = 1}^{\infty} [\frac{z}{k} - \ln (\frac{z}{k} + 1)]

This series converges very slowly so it can be helped by the Euler-Maclaurin summation formula. The Euler-Maclaurin summation formula is:

\sum_{k = 0}^{m} F (a + k h) = \frac{1}{h} \int_{a}^{b} F (t) dt + \frac{1}{2} (F (b) + F (a)) + \sum_{k = 1}^{n - 1} \frac{h^{2 k - 1}}{(2 k)!} B_{2 k} (F^{(2 k - 1)} (b) - F^{(2 k - 1)} (a)) + e_{n}

B_{2 k}

is a Bernoulli number,

b = a + m h

, and

e_{n}

is the error. The way to use the Euler-Maclaurin summation is to get the value of

a

large enough so the summation converges fast. The summation can be split into two parts.

\sum_{k = 1}^{\infty} [\frac{z}{k} - \ln (\frac{z}{k} + 1)] = \sum_{k = 1}^{a - 1} [\frac{z}{k} - \ln (\frac{z}{k} + 1)] + \sum_{k = a}^{\infty} [\frac{z}{k} - \ln (\frac{z}{k} + 1)]

Put the second part in a form to use the Euler-Maclaurin summation formula.

\sum_{k = a}^{\infty} [\frac{z}{k} - \ln (\frac{z}{k} + 1)] = \sum_{k = 0}^{\infty} [\frac{z}{a + k} - \ln (\frac{z}{a + k} + 1)]

m \to \infty, h = 1, F (a + k h) = [\frac{z}{a + k h} - \ln (\frac{z}{a + k h} + 1)]

Define the function.

F (x) = [\frac{z}{x} - \ln (\frac{z}{x} + 1)]

Find the derivatives of the function.

F^{(1)} (x) = - \frac{1}{x + z} + \frac{1}{x} - \frac{z}{x^{2}}

F^{(3)} (x) = - \frac{2}{(x + z)^{3}} + \frac{2}{x^{3}} - \frac{6 z}{x^{4}}

:

F^{(2 k - 1)} (x) = (- 1)^{2 k - 1} (2 k - 1)! (\frac{(x + z)^{- 2 k + 1} - x^{- 2 k + 1}}{2 k - 1} + z x^{- 2 k})

Calculate the integral of the function.

{\int_{a}^{b} F (t) dt = - (z + t) \ln (\frac{z}{t} + 1) |}_{a}^{b}

Some of the terms in the Euler-Maclaurin summation drop out.

lim_{x \to \infty +} (- (z + x) \ln (\frac{z}{x} + 1)) = z

lim_{x \to \infty +} F (x) = 0

lim_{x \to \infty +} (\frac{(x + z)^{- 2 k + 1} - x^{- 2 k + 1}}{2 k - 1} + z x^{- 2 k}) = 0 (k > 0)

Substitute the three remaining parts into the summation. The Euler-Maclaurin summation becomes:

\sum_{k = 0}^{\infty} [\frac{z}{a + k} - \ln (\frac{z}{a + k} + 1)] = (z + a) \ln (\frac{z}{a} + 1) - z +

\frac{[\frac{z}{x} - \ln (\frac{z}{x} + 1)]}{2} + \sum_{k = 1}^{n - 1} \frac{B_{2 k}}{(2 k)!} [(2 k - 1)! (\frac{(a + z)^{- 2 k + 1} - a^{- 2 k + 1}}{2 k - 1} + z a^{- 2 k})] + e_{n}

Substitute the Euler-Maclaurin summation into the split series for the logarithm of the gamma function. The logarithm of the gamma function becomes:

\ln (Γ (z)) = - \ln (z) - γ z + \sum_{k = 1}^{a - 1} [\frac{z}{k} - \ln (\frac{z}{k} + 1)] + (z + a) \ln (\frac{z}{a} + 1) - z +

\frac{[\frac{z}{a} - \ln (\frac{z}{a} + 1)]}{2} + \sum_{k = 1}^{n - 1} \frac{B_{2 k}}{2 k} (\frac{(a + z)^{- 2 k + 1} - a^{- 2 k + 1}}{2 k - 1} + z a^{- 2 k}) + e_{n}

Take the exponential function of both sides.

Γ (z) = z^{- 1} Π_{k = 1}^{a - 1} [\frac{e^{\frac{z}{k}}}{(\frac{z}{k} + 1)}] {(\frac{z}{a} + 1)}^{z + a} {(\frac{z}{a} + 1)}^{- \frac{1}{2}}

e^{\frac{z}{2 a} - γ z - z + \sum_{k = 1}^{n - 1} \frac{B_{2 k}}{2 k} (\frac{(a + z)^{- 2 k + 1} - a^{- 2 k + 1}}{2 k - 1} + z a^{- 2 k}) + e_{n}}

File translated from T_EX by T_TM, version 3.01.
On 5 Jul 2001, 00:33.