Factorial, Gamma Function, and More

Cover image
🏛 A Project Cauchy Op-ed

Named after French mathematician Augustin-Louis Cauchy, Project Cauchy column is where I invite some of the HFI Programming club members to provide neat proofs or explanations about some number theory puzzles. I highly suggest that you read these articles with a pencil and paper so you can sketch things out and scribble solutions to exercises as you come across them. This week, Travor is demonstrateing us one of the most commonly used extension of the factorial function to complex numbers.

It is well-known that factorial is defined by the following recursive relation,


with 0!=10!=1, but however it is possible to generalize this operation to complex numbers. Let's begin our journey of generalizations!

First generalization: integral representation

It can be easily shown that the following integral satisfies

0eλtdt=1λ\int_0^\infty e^{-\lambda t}\mathrm dt=\frac1\lambda

Differentiation with respect to ss on both side for n1n-1 times gives

0tn1eλtdt=(n1)!λn\int_0^\infty t^{n-1}e^{-\lambda t}\mathrm dt={(n-1)!\over\lambda^n}

Setting λ=1\lambda=1 gives us the first generalization of factorial: the Gamma function.

(s1)!=Γ(s)=0ts1etdt(s-1)!=\Gamma(s)=\int_0^\infty t^{s-1}e^{-t}\mathrm dt

In fact, by integration by parts we can show that the Gamma function satisfies the recursive relationship:


Furthermore, we have the relation

0ts1eλtdt=Γ(s)λs\int_0^\infty t^{s-1}e^{-\lambda t}\mathrm dt={\Gamma(s)\over\lambda^s}

Digression: the Beta function

To convenience our derivation process, let's define the Euler integral of the first kind: the Beta function B(x,y)B(x,y):

B(x,y)=01τx1(1τ)y1dτ(2)B(x,y)=\int_0^1\tau^{x-1}(1-\tau)^{y-1}\mathrm d\tau\tag2

which can be seen as a convolution between two power functions:

B(x,y;t)=0tτx1(tτ)y1dτB(x,y;t)=\int_0^t\tau^{x-1}(t-\tau)^{y-1}\mathrm d\tau

If we were to perform Laplace transform on both side, we get

L{B(x,y;t)}(λ)=Γ(x)Γ(y)λx+y\mathscr L\{B(x,y;t)\}(\lambda)={\Gamma(x)\Gamma(y)\over\lambda^{x+y}}

If we juxtapose this fact with the relation

L{tx+y1}(λ)=Γ(x+y)λx+y\mathscr L\{t^{x+y-1}\}(\lambda)={\Gamma(x+y)\over\lambda^{x+y}}

then we see that


which will be extremely useful for us to generalize factorial even further.

Γ(s)\Gamma(s) evaluated at complex numbers

Via taking absolute variable, we know that the Gamma integral converges whenever (s)>0\Re(s)>0

Γ(s)0t(s)1etdt|\Gamma(s)|\le\int_0^\infty t^{\Re(s)-1}e^{-t}\mathrm dt

which means that this improper integral is converges absolutely on the right half plane. Hence, a stronger definition is needed for us to expand it to the entire complex plane.

Γ(s)\Gamma(s) as a limit

Let's consider a sequence of functions:

fn(t,s)={ts1(1tn)n0xn0x>nf_n(t,s)= \begin{cases} t^{s-1}\left(1-\frac tn\right)^n & 0\le x\le n \\ 0 & x>n \end{cases}

Then by the exponential limit we see that

fn(t,s)ts1(1tn)nf_n(t,s)\to t^{s-1}\left(1-\frac tn\right)^n

in a pointwise sense. However, it is possible to show that this sequence converges uniformly for t[0,+)t\in[0,+\infty). Let's set T>0T>0 such that for all t>Tt>T we have


Then, let's consider the interval [0,T][0,T]:

ts1etfn(t,s)=t(s)1et(1tn)n=t(s)1etebn(t)\begin{aligned} |t^{s-1}e^{-t}-f_n(t,s)| &=t^{\Re(s)-1}\left|e^{-t}-\left(1-\frac tn\right)^n\right| \\ &=t^{\Re(s)-1}\left|e^{-t}-e^{-b_n(t)}\right| \end{aligned}

where we define bn(t)b_n(t) as

bn(t)=nlog(1tn)b_n(t)=-n\log\left(1-\frac tn\right)

By the uniform continuity of ete^{-t} on [0,T][0,T], all we need is to prove that bn(t)b_n(t) converges uniformly to tt. In fact, for all n>tn>t we can use the Taylor expansion of logarithm to obtain

bn(t)t=nlog(1tn)t=n[tn+O(1n2)]t=O(1n)\begin{aligned} |b_n(t)-t| &=\left|-n\log\left(1-\frac tn\right)-t\right| \\ &=\left|-n\left[\frac tn+\mathcal O\left(1\over n^2\right)\right]-t\right| \\ &=\mathcal O\left(\frac1n\right) \end{aligned}

As a result, we conclude that fn(s,t)f_n(s,t) converges uniformly to ts1ett^{s-1}e^{-t}, which allows us to interchange the limit operation and integral to obtain

limn0nts1(1tn)ndt=Γ(s)(3)\lim_{n\to\infty}\int_0^n t^{s-1}\left(1-\frac tn\right)^n\mathrm dt=\Gamma(s)\tag3

In the following procedure, we are going to expand the left hand side limit in a subtle sense so that the right hand side can be analytically continued to the left half plane.

Expanding the integral sequence

Performing a change of variables on (3) gives

Γn(s)0nts1(1tn)ndtt=nr=ns01rs1(1r)ndr\Gamma_n(s)\triangleq\underbrace{\int_0^nt^{s-1}\left(1-\frac tn\right)^n\mathrm dt}_{t=nr}=n^s\int_0^1r^{s-1}(1-r)^n\mathrm dr

In fact, the right hand side can be expressed by Beta function as in (3), which yields


Continuous application of (1) gives us

Γn(s)=Γ(s)nsn!Γ(s+n+1)=nsn!s(s+1)(s+2)(s+n)=nssk=1nks+k\begin{aligned} \Gamma_n(s) &={\Gamma(s)n^sn!\over\Gamma(s+n+1)} \\ &={n^sn!\over s(s+1)(s+2)\cdots(s+n)} \\ &={n^s\over s}\prod_{k=1}^n{k\over s+k} \\ \end{aligned}

which transforms Gamma function into a product representation, however this expression still looks ugly, why not go deeper?

Weierstrass product for Γ(s)\Gamma(s)

First, let's turn this equation up side down to obtain

1Γn(s)=snsk=1n(1+sk)=seslognk=1n(1+sk)\begin{aligned} {1\over\Gamma_n(s)} &=sn^{-s}\prod_{k=1}^n\left(1+\frac sk\right) \\ &=se^{-s\log n}\prod_{k=1}^n\left(1+\frac sk\right) \end{aligned}

Then, employing the fact that Hnk=1n1k=logn+γ+O(1/n)H_n\triangleq\sum_{k=1}^n\frac1k=\log n+\gamma+\mathcal O(1/n) we can replace the logarithm with harmonic numbers:

1Γn(s)=ses[Hnγ+O(1/n)]k=1n(1+sk)=sesγ+O(1/n)k=1n(1+sk)es/k\begin{aligned} {1\over\Gamma_n(s)} &=se^{-s[H_n-\gamma+\mathcal O(1/n)]}\prod_{k=1}^n\left(1+\frac sk\right) \\ &=se^{s\gamma+\mathcal O(1/n)}\prod_{k=1}^n\left(1+\frac sk\right)e^{-s/k} \end{aligned}

Now, if we take logarithm on the product, we have

logk=1n(1+sk)es/ks2k=1n1k2<s2π26\begin{aligned} \log\left|\prod_{k=1}^n\left(1+\frac sk\right)e^{-s/k}\right| &\le|s^2|\sum_{k=1}^n{1\over k^2}<{|s|^2\pi^2\over6} \end{aligned}

As a result, the product converges absolutely for all sCs\in\mathbb C, giving us the Weierstrass product representation of Gamma function:

1Γ(s)=seγsk=1(1+sk)es/k{1\over\Gamma(s)}=se^{\gamma s}\prod_{k=1}^\infty\left(1+\frac sk\right)e^{-s/k}

which allows us to analytically continue Γ(s)\Gamma(s) to the entire complex planes except for nonpositive integers:

Γ(s)=eγssk=1(1+sk)1es/k(4)\Gamma(s)={e^{-\gamma s}\over s}\prod_{k=1}^\infty\left(1+\frac sk\right)^{-1}e^{s/k}\tag4

Remark: (4) also reveals Γ(s)\Gamma(s) is non-zero for all sCs\in\mathbb C

Now, we successfully expanded factorial to the entire complex plane except at negative integers, but Gamma function has some other brilliant properties. Let's have a look at some of them:

Alternative Definition for Euler-Mascheroni constant

From the last article, we know that Euler-Mascheroni constant is defined by

γ=limn(Hnlogn)=11{x}x2dx\gamma=\lim_{n\to\infty}(H_n-\log n)=1-\int_1^\infty{\{x\}\over x^2}\mathrm dx

However, it is possible for us to create a new definition for γ\gamma by using Γ(s)\Gamma(s).

Remark: I hypothesize this explains why they name the function Γ\Gamma and the constant γ\gamma since they are highly correlated.

To begin with, we take logarithm on (4) to get

logΓ(s)=logsγs+k=1{sklog(1+sk)}\log\Gamma(s)=-\log s-\gamma s+\sum_{k=1}^\infty\left\{\frac sk-\log\left(1+\frac sk\right)\right\}

If we were to define Digamma function ψ(s)\psi(s) as the logarithmic derivative of Γ(s)\Gamma(s), then

ψ(s)=γ1s+k=1{1k1s+k}\psi(s)=-\gamma-\frac1s+\sum_{k=1}^\infty\left\{\frac1k-{1\over s+k}\right\}

If we were to move the 1s\frac1s term into the summation, we deduce the standard definition for Digamma function.

ψ(s)=γ+m=0{1m+11m+s}(5)\psi(s)=-\gamma+\sum_{m=0}^\infty\left\{{1\over m+1}-{1\over m+s}\right\}\tag5

Plugging s=1s=1 gives ψ(1)=Γ(1)/Γ(1)=γ\psi(1)=\Gamma'(1)/\Gamma(1)=-\gamma. Because Γ(1)=1\Gamma(1)=1, we also know that Γ(1)=γ\Gamma'(1)=-\gamma. Combining this with the original integral definition for Γ(s)\Gamma(s) gives this elegant integral identity to represent Euler-Mascheroni constant:

0etlog(t)dt=γ(6)\int_0^\infty e^{-t}\log(t)\mathrm dt=-\gamma\tag6

Integral representation for Digamma function

While deriving (5) in the previous section, we introduce the Digamma function, so why don't we do some calculus on that

ψ(s)=γ+m=001(xmxm+s1)dt=γ+01(1xs1)m=0xmdt=γ+011xs11xdt\begin{aligned} \psi(s) &=-\gamma+\sum_{m=0}^\infty\int_0^1(x^m-x^{m+s-1})\mathrm dt \\ &=-\gamma+\int_0^1(1-x^{s-1})\sum_{m=0}^\infty x^m\mathrm dt \\ &=-\gamma+\int_0^1{1-x^{s-1}\over1-x}\mathrm dt \end{aligned}

Well, the constant lying outside is not beautiful, so why not use (6):

γ=0etlog(t)dtx=et=10log(logx)dx=01loglog(1x)dx\begin{aligned} -\gamma &=\underbrace{\int_0^\infty e^{-t}\log(t)\mathrm dt}_{x=e^{-t}} \\ &=-\int_1^0\log(-\log x)\mathrm dx \\ &=\int_0^1\log\log\left(\frac1x\right)\mathrm dx \end{aligned}

Combining all these gives us the ultimate integral definition for ψ(s)\psi(s):

ψ(s)=01{loglog(1x)+xs11x1}dx(7)\psi(s)=\int_0^1\left\{\log\log\left(\frac1x\right)+{x^{s-1}-1\over x-1}\right\}\mathrm dx\tag7


In this blog, we first use the technique of differentiation under integral to deduce an integral representation for factorial that introduces the concept of Gamma function Γ(s)\Gamma(s). Then, by using an identity connecting B(x,y)B(x,y) and Γ(s)\Gamma(s), we obtain a product formula that turns Γ(s)\Gamma(s) into a meromorphic function on C\mathbb C. Using this newly obtained product formula, we are also able to discover some new identities. In fact, Gamma function is a function that often appears in the field of analytic number theory, and we will begin future investigation based on the following identity (which you could try prove it yourself):

n=11ns=1Γ(s)0xs1ex1dx\sum_{n=1}^\infty{1\over n^s}={1\over\Gamma(s)}\int_0^\infty{x^{s-1}\over e^x-1}\mathrm dx
◀ Extending FGSM to Other NormsAdversarial Camouflage ▶