Repository navigation
Gamma and Related Functions
Gamma functions were not originally part of this library in the pre-releases, but have been in it since the first 1.0 release.
The functions in this library that have been implemented and which pertain to the gamma function are:
Complex gamma(Complex z) the gamma function, Γ(z)
Complex factorial(Complex z) the factorial function, z!
Complex loggamma(Complex z) the log-gamma function, lnΓ(z)
Complex digamma(Complex z) the digamma function, ψ(z) (sometimes written ψ₀(z))
Complex polygamma(int m, Complex z) the polygamma function, ψ(m,z) (sometimes written ψ_m(z), where _m indicates a subscript m)
Complex loggammaReflector(Complex z) lnΓ(z)+Γ(1-z)
These functions are all fundamentally tied to the gamma function, one of the most well-known non-elementary mathematical functions, possibly next to the Riemann Zeta Function. For that reason, I will start by explaining the Gamma function.
The gamma function, Γ(z), is an extension of the factorial. The factorial of an integer is defined as the product between that number and all positive integers before it, with the special case that 0! is 1 for reasons we'll get to in a bit.
There is no agreed upon consensus on whether or not it is acceptable to perform the factorial on a non-integer. Some say you can, other's say you can't. In any case, all mathematicians agree that you can perform the Gamma function on any number (except non-positive integers, because that either gives ±∞ or undefined). The gamma function is essentially just the factorial, but no one says you can't use it on integers. There's a bit of nuance to that statement, but we'll circle back to that in a moment.
The first question that comes to mind is "how does it even make sense to perform a factorial on a non-integer?". That's a really good question. Mathematics is built upon patterns. Mathematics is really just the study of things that follow patterns. And sure enough, factorials follow a basic pattern. Observe:
1! = 1
2! = 1*2
3! = 1*2*3
4! = 1*2*3*4
5! = 1*2*3*4*5
...
16000! = 1*2*3*4*5*...*15999*16000
If we want, though, we can rewrite this as such:
1! = 1
2! = 1!*2
3! = 2!*3
4! = 3!*4
5! = 4!*5
...
16000! = 15999!*16000
From this, we can define the factorial of all integers recursively by saying:
1! = 1
n! = (n-1)!*n
In fact, this very fact explains precisely why 0! must be equal to 1.
1! = 1
1! = (1-1)!*1 = 0!*1
1 = 0!*1
1 = 0!
0! = 1
We can also use this to define the factorial of all negative integers. Though the result might be a bit unsatisfying:
0! = (-1)!*0
1 = (-1)!*0
1/0 = (-1)!
(-1)! = ∞ (or undefined)
(-1)! = (-2)!*-1
-(-1)! = (-2)!
(-2)! = -∞ (or undefined)
(-2)! = (-3)!*-2
-(-2)!/2 = (-3)!
(-3)! = ∞/2 (or undefined)
Thus, the factorial of all negative integers is either ±∞ or undefined.
But how do we find the factorial of non-integers? Well, that's a good question. To make a long story short, we just do. But if that's not a satisfying explanation for you, consider the following. There's another pattern that the factorial creates. Consider the factorial of a really large number. Say, googol (1 followed by 100 0s; 10^100).
googol! = 1*2*3*4*5*...*googol
(googol+1)! = googol! * (googol+1) ≈ googol! * googol
(googol+2)! = (googol+1)! * (googol+2) ≈ (googol+1)! * googol ≈ googol! * googol * googol = googol! * googol²
(googol+3)! = (googol+2)! * (googol+3) ≈ (googol+2)! * googol ≈ googol! * googol² * googol = googol! * googol³
(googol+n)! ≈ googol! * googol^n, so long as n is much smaller than googol.
Note, the expression above is slightly off. Just slightly. As the number gets bigger and bigger, saying "multiplying by x+1 is basically the same as multiplying by x" gets a lower and lower percent error. Using googol, the percent error is on the order of 10^-100. If we used googolplex instead of googol, the percent error would be on the order of 10^(-10^100), and if we used googolplexian (1 followed by googolplex 0s), the percent error would be impractical to write. And in the limiting case as we go off to ∞, the error goes to 0. If we follow this pattern, we can safely say that
(x+n)! ≈ x! * x^n
for all sufficiently large values of x.
If we extend that pattern to hold regardless of if n is an integer, we can define the factorial for non-integers. For instance, (googol+1/2)! ≈ googol! * √(googol). Again, this isn't exact, but it's pretty close. If we could somehow do enough multiplication to compute the factorial of googol, then using this expression, we could get a value for (googol+0.5)! that's accurate to the first 99 digits. And if we calculated the factorial of googolplex, we could use (googolplex+1/2)! ≈ googolplex! * √(googolplex) to obtain an approximation for that factorial accurate to the first googol digits, which is very impressive.
So this lets us (approximately) calculate the factorial of non-integers...but only for large non-integers. How can we use this to calculate the factorial of, say, 1/2? Well, observe:
If n! = (n-1)!*n, that means (n-1)! = n!/n, which means n! = (n+1)!/(n+1)
(googol-1/2)! = (googol+1/2)!/(googol+1/2)
(googol-3/2)! = (googol-1/2)!/(googol-1/2) = (googol+1/2)!/((googol+1/2)(googol-1/2))
(googol-5/2)! = (googol-3/2)!/(googol-3/2) = (googol+1/2)!/((googol+1/2)(googol-1/2)(googol-3/2))
...
(1/2)! = (3/2)!/(3/2) = (googol+1/2)!/((googol+1/2)(googol-1/2)(googol-3/2)(googol-5/2)...(7/2)(5/2)(3/2))
(googol+1/2)! ≈ googol! * √(googol)
(1/2)! ≈ googol! * √(googol)/((googol+1/2)(googol-1/2)(googol-3/2)(googol-5/2)...(7/2)(5/2)(3/2))
(1/2)! ≈ 1*2*3*4*5*...*googol * √(googol)/((googol+1/2)(googol-1/2)(googol-3/2)(googol-5/2)...(7/2)(5/2)(3/2))
(1/2)! ≈ 1*2*3*4*5*...*googol * √(googol)/((3/2)(5/2)(7/2)...(googol-5/2)(googol-3/2)(googol-1/2)(googol+1/2))
Note, googol was arbitrary. (1/2)! is actually the limiting case as this large number approaches ∞.
(1/2)! = the limit as W->∞ of 1*2*3*4*5*...*W * √(W) / ((3/2)(5/2)(7/2)...(W-3/2)(W-1/2)(W+1/2))
It should also be noted that this can be done for any number.
(1/3)! = lim W->∞ 1*2*3*4*5*...*W * W^(1/3) / ((4/3)(7/3)(10/3)...(W-5/3)(W-2/3)(W+1/3))
(1/4)! = lim W->∞ 1*2*3*4*5*...*W * W^(1/4) / ((5/4)(9/4)(13/4)...(W-7/4)(W-3/4)(W+1/4))
And, in general,
n! = lim W->∞ 1*2*3*4*5*...*W * W^n / ((n+1)(n+2)(n+3)...(n+W-2)(n+W-1)(n+W))
This even works for integers, surprisingly:
3! = lim W->∞ 1*2*3*4*5*6*7*8*...*W * W^3 / (4*5*6*7*8*...*(W-3)(W-2)(W-1)W(w+1)(W+2)(W+3))
If you cancel out the numerator and denominator, we get
3! = lim W->∞ 1*2*3 *W^3 / ((W+1)(W+2)(W+3)) = 1*2*3 = 6
And if we try it for half integers...well, the answer isn't quite straightforward, but if you either wait long enough for the number to converge to something, or are just really good at math, you'll find
(-1/2)! = √(π)
(1/2)! = √(π)/2
(3/2)! = 3√(π)/4
Weird, you wouldn't expect π to come into play would you? Let alone the square root of it. Also, if you're curious, there's nothing special about (1/3)!, (2/3)!, (1/4)!, (3/4)!, (1/5)!, etc. The only factorials that evaluate to something involving familiar numbers are integers and half-integers. For some reason...
And also, this pattern works for complex numbers, hence it is completely possible to extend this definition to complex numbers.
Now, before, I said there was a bit of nuance to the gamma function. You see, technically, the gamma function Γ(z) is not the same as z!. Actually, Γ(z) equals (z-1)!. Why? I have no clue. Mathematicians say it's more "natural" this way. Because now, it's the product of all numbers up to z, not including z. It would actually be more natural to define Γ(z) as (z-1/2)! (for reasons I won't go into, defining it that way causes it to obey SEVERAL symmetries). In any case, mathematicians heard our plea, that we wanted to extend the factorial to non-integers, not just this new fancy gamma function. And they invented the PI function, Π(z). A function which no one uses because they just find it easier to call it z!.
In any case, the gamma function obeys several formulas. Just as z! = (z-1)!*z, Γ(z) = Γ(z-1)*(z-1). Or, alternatively,
Γ(z)*z = Γ(z+1)
We also have what is known as the reflection formula:
Γ(z)*Γ(1-z) = π/sin(πz).
We also have what is known as the duplication formula:
Γ(z)*Γ(z+1/2) = √(π)*Γ(2z)/2^(2z-1)
Finally, the formula that causes the gamma function to come up so often:
Γ(z) = ∫ e^-t*t^(z-1) dt from 0 to ∞ (it's a definite integral, but I couldn't put in the bounds without using LaTex, something not supported by this wiki).
A lot of definite integrals, when appropriate substitution is performed, evaluate to this integral with some value for z. This causes many patterns to be explained using the gamma function, and causes it to be used so frequently.
I've already told you what the gamma function does. Hopefully you understood my explanation. If not, I do apologize. Likewise, the factorial function is just the gamma function of the input plus 1.
The loggamma function is essentially just the logarithm of the gamma function. There are a few good reasons for this function to be here, actually. For starters, if we take the logarithm of both sides, we find that the identity
Γ(z+1) = Γ(z)*z
turns into
lnΓ(z+1) = lnΓ(z)+ln(z)
Alternatively, this means
lnΓ(x) = ln(1)+ln(2)+ln(3)+...+ln(x-1) for integer values of x.
However, in order for this relation to hold up, we have to make a slight adjustment. lnΓ(z) can't just be the principal value of the logarithm. Sometimes, it has to be the logarithm plus some multiple of 2πi. For instance, the gamma function of -1.5 is about 2.36327. The logarithm of that is 0.860047. However, lnΓ(-1.5) is 0.860047-2πi. Why is this? It's complicated. The point is, the loggamma function has fewer discontinuities than naively taking the principal value of the logarithm of the gamma function.
Something also worth noting, however, is that the gamma function grows pretty big pretty fast. By the time we get to z=172, the output is already too big to represent with double precision floating points. The gamma function grows at a rate of about z^z. However, the loggamma function grows at a rate of about zln(z) (linearithmic). As such, it doesn't grow nearly as fast, and you don't reach an overflow until the input is about equal to 10^305. For reference, the maximum double is 1.7*10^308.
The digamma function is the derivative of the loggamma function. While the loggamma function obeys the relation:
lnΓ(z+1) = lnΓ(z)+ln(z)
The digamma function obeys the relation:
ψ(z+1) = ψ(z)+1/z
If you take the derivative of the first expression, you'll find it equals the second expression, so hopefully this should all make sense.
The digamma function comes up quite a bit, though not nearly as much as the gamma function. Something of note is that, if z is an integer, ψ(z) = -γ + 1/1 + 1/2 + 1/3 + 1/4 + ... + 1/(z-1), where γ is the Euler-Mascheroni constant. In this library, you can obtain this constant by saying Mafs.GAMMA, in the same way you can obtain π by saying Math.PI or e by saying Math.E. For those of you familiar with the harmonic series, you can obtain the harmonic series from 1 to n by asking for ψ(n+1)+γ. In this library, you'd do that by saying Cpx2.digamma(new Complex(n+1)).re+Mafs.GAMMA. For large integers, this is much more efficient than naively adding up all the reciprocals from 1 to n.
The digamma function grows at a rate of about ln(z). Because the harmonic series never converges, when you plug in ∞, you get ∞.
The trigamma function is the derivative of the digamma function. It's often denoted as ψ₁(z). It doesn't have it's own function name in this library (yet???), but it can be calculated via polygamma(1,z). The reason I'm explaining the trigamma function, though, is because it's a good introduction into the concept of the polygamma function.
Just as the digamma function obeys the relation:
ψ(z+1) = ψ(z)+1/z
The trigamma function obeys the relation:
ψ₁(z+1) = ψ₁(z)-1/z²
The trigamma function doesn't come up as much as the digamma function, though a lot of situations which use the digamma function could be rearranged to another problem which uses the trigamma function. For instance, the mean amount of time needed in the infamous coupon collector's problem to finish can be calculated via the digamma function. However, the standard deviation for that time can be calculated with the trigamma function. In addition, for integer values of z, ψ₁(z) = π²/6 - 1/1²-1/2²-1/3²-...-1/(z-1)². For those of you familiar with the Basel problem, this π²/6 term might seem familiar to you. Unlike with the digamma function, which goes off to ∞ as x goes to ∞, the trigamma function goes off to 0. Since the sum of the squares of the reciprocals of all positive integers is π²/6, that term must be π²/6 for said relation to hold. In addition, if someone wants to compute the Basel sum 1/1²+1/2²+1/3²+...+1/n², they can simply compute π²/6-ψ₁(n+1). In this library, that would be Math.PI*Math.PI/6-polygamma(1,new Complex(n+1)).re. For large values of n, this is more efficient than naively adding up the reciprocal squares up to n.
The trigamma function shrinks at a rate of 1/z. When you plug in ∞, you get 0.
The polygamma function is the m-th derivative of the digamma function, often denoted ψ_m(z) (_m meaning subscript m), but for simplicity we'll use ψ(m,z). In this library, m must be an integer. While there is a way of defining it for non-integers, it's extremely complicated, and due to certain ambiguities, it'd likely be more convenient to use the Hurwitz Zeta function.
As the digamma function obeys the relation:
ψ(z+1) = ψ(z)+1/z
The polygamma function obeys the relation:
ψ(m,z+1) = ψ(m,z)+(-1)^m*m!/z^(m+1)
The m-th derivative of the first expression is equal to the second expression. When m is 1, this just yields the trigamma function. When m is 0, this just yields the digamma function. When m is -1, this just yields the loggamma function. When m is -2, this is supposed to yield some variation of the logarithm of the Barnes G function (or the logarithm of the K function). However, implementing such a function would require calculation of the polylogarithms, which have yet to be implemented. As such, using m=-2 will yield NaN until said functions are implemented. Likewise, using an m less than -2 will likely always yield NaN since the definition for the polygamma function for negative parameters is very difficult to work with.
When m is greater than or equal to 1, this function will always converge to 0 as the input goes to ∞. Otherwise, it'll go off to ∞.
The polygamma function is a useful tool for computing sums of negative integer powers. For instance, 1/1³+1/2³+1/3³+...+1/n³ = ψ₂(n+1)/2-ζ(3).
As was mentioned before, the gamma function obeys the relation:
Γ(z)Γ(1-z) = π/sin(πz)
As such, it would be easy to take the logarithm of both sides and assume the loggamma function obeys the relation:
lnΓ(z)+lnΓ(1-z) = ln(π) - ln(sin(πz))
However, this is not true, it actually obeys the relation:
lnΓ(z)+lnΓ(1-z) = ln(π) - ln(sin(πz)) + 2πi*N, where N is some integer.
Calculating this integer takes a lot of work, surprisingly. I spent hours analyzing the function, trying to figure out what it must be in order for the relation
lnΓ(z+1) = lnΓ(z)+ln(z)
to be preserved, while also making sure the function was as smooth as possible. I initially was going to make the end result a private function, just like many other weird functions. The reason I made this function to begin with is because all the algorithms used to compute gamma and related functions fail when the input has a negative real part, so you have to use a reflection formula to calculate for those values. However, this function works for all inputs, it doesn't have any special cases, so I figured there wasn't really any harm in making it public. But how could you, the programmer, possibly use this? Well, for starters, if you were to try to plot out the 3D function z = Im(ln(sin(x+yi))), you'd find that it constantly jumps up and down by π. However, if you plot out the function z = -Im(loggammaReflector((x+yi)/π)), you'd find it to be mostly smooth. There's an unavoidable discontinuity when crossing the x-axis, but other than that, smooth! If you ever need a smooth version of the log sine function, this should come in handy. It should be noted that this function is similar to the gudermannian function. The gudermannian is the integral of the hyperbolic secant, and the inverse gudermannian is the integral of the secant. Likewise, the logarithm of the sine is the integral of the cotangent, and the logarithm of the hyperbolic sine is the integral of the hyperbolic cotangent. If you're comfortable with shifting by π/2, the logarithm of the cosine is the integral of negative tangent, and the logarithm of the hyperbolic cosine is the integral of the hyperbolic tangent.
To obtain a semi-continuous mapping of the logarithm of the sine, simply compute ln(π)-loggammaReflector(z/π). In code, that's Cpx.sub(Math.log(Math.PI), Cpx2.loggammaReflector(z.div(Math.PI))).