# Gamma Function Golf

Given a real number t in (-10^9,13) (not including -10^9 or 13) as input, output Γ(t), also known as the Gamma function, which is defined as follows:

You may not use a built-in Gamma function to solve this task, nor may you use built-in numeric or symbolic integration functions. Your output should be accurate to 6 significant figures or within 10^-6 of the actual value, whichever is less restrictive for the given value. Python's built-in Gamma function will be used for determining the actual value. You may assume Γ(t) is defined - that is, t is either a positive real number or a non-integer negative real number - and that |Γ(t)| ≤ 10^9. Here is a reference program that you may use to get the actual values, using Python's built-in Gamma function.

### Examples

1 -> 1.000000
-2.5 -> -0.945309
3.14159265 -> 2.288038
-2.71828182846 -> -0.952682
12 -> 39916800.000000
0.5 -> 1.772454
8.675309 -> 20248.386956
-10.1 -> -0.000002

## Rules

• This is , so shortest answer (in bytes) wins.
• Standard loopholes are forbidden.
• Input and output may be performed in whatever manner is considered standard for your language.
• You may write a full program, a function, or anything that is normally considered a valid answer for your language

The Stack Snippet at the bottom of this post generates the leaderboard from the answers a) as a list of shortest solution per language and b) as an overall leaderboard.

## Language Name, N bytes

where N is the size of your submission. If you improve your score, you can keep old scores in the headline, by striking them through. For instance:

## Ruby, <s>104</s> <s>101</s> 96 bytes

If there you want to include multiple numbers in your header (e.g. because your score is the sum of two files or you want to list interpreter flag penalties separately), make sure that the actual score is the last number in the header:

## Perl, 43 + 2 (-p flag) = 45 bytes

You can also make the language name a link which will then show up in the snippet:

## [><>](http://esolangs.org/wiki/Fish), 121 bytes

• Please provide explicit bounds for t such that |gamma(t)|<10^9 – flawr Nov 14 '15 at 23:49
• The link is not a reference implementation,... – sergiol Mar 9 '18 at 0:03
• @sergiol Reworded it – user45941 Mar 9 '18 at 0:04

# Pyth, 21 bytes

As with my TI-BASIC answer, I haven't been able to test this with the full 8^10 iterations, but everything seems good with smaller cases.

cu*Gc^hc1HQhcQHS^8T1Q

Explanation:

[implicit: Q=input]
^8T         8**10
S^8T         [1,2,3,...,8**10]
*Gc^hc1HQhcQH             lambda G,H:G*(1+1/H)**Q/(1+Q/H)
1        Base case
u*Gc^hc1HQhcQHS^8T1        Reduce with base case 1
c                   Q       Divide by Q

Try it here with 2000 iterations instead of 8^10.

# C++14, 8685 81 bytes

[](auto t){auto v=1.;for(int x=1;x<1e9;++x)v*=pow(1+1./x,t)/(1+t/x);return v/t;};

I didn't spend much time on this one. I just looked on approximation that seemed the easiest to implement (in the manner of bytes). It will take some time to compute the value (since the loop is over all positive integers), but time limitation is not specified in the challenge. It is anonymous function (lambda), which takes any argument (convertible to T on which pow(double, T) and operator/(T,int) can be called) and returns double.

Ungolfed with usage

#include <iostream>
int main()
{
auto r = [](auto t)
{
auto v = 1.;
for (int x = 1; x < 1e9; ++x)
v *= pow(1 + 1. / x, t) / (1 + t / x);
return v / t;
};
std::cout << r(-2.71828182846); // outputs -0.952682
}
• @Mego Of course it is! Thanks. – Zereges Nov 14 '15 at 23:38
• So what value o you get for -10^9 and for 10^9? I first want to know how well your stuff works before you get my upvote. – flawr Nov 14 '15 at 23:40
• @Mego Microsoft compiler does not need neither of theese includes. – Zereges Nov 14 '15 at 23:46
• @Mego Microsoft (R) C/C++ Optimizing Compiler Version 19.00.23026 for x86 – Zereges Nov 14 '15 at 23:49
• @flawr Mine program outputs 0 for gamma(-10e9) but OP stated, that only parameters, for which gamma function is defined, can be considered. gamma(10e9) returns inf, while Python's built-in Gamma function will be used for determining the actual value says OverflowError: math range error – Zereges Nov 14 '15 at 23:59

## Minkolang 0.12, 3534 25 bytes

n$zl8;dz;z$:r[i1+dz+$:*]N This does halt with an error (on trying to divide by 0), but that is allowed as per Meta consensus. Add a . at the end for a program that halts normally. Try all test cases at once. (The loop iterates only 1e4 times so it would finish sooner rather than later.) ### Explanation Zereges used one of the alternative, infinite product definitions. As it turns out, the other is much more amenable to implement in Minkolang. This is a limit as n goes to infinity, which means that I can calculate both n! and (t+n) as I go. So I take out 1/t (because 0!=1) and n^t because that one cannot be calculated sequentially without knowing the ending value of n. As it happens, because n is the limit, I can use it twice. Once as a factor in the calculation and once as the number of times to run the loop. An sequential infinite product has to start with something, usually 1. In this case, it's n^t/t. In the body of the loop, I calculate k/(t+k) and multiply this with the product so far. At the end, the whole product has been calculated and output. This is essentially what my program does, with n high enough that the answer is precise enough. n Take number from input$z                          Store it in the register (this is t; retrieved with z)
l8;                       10^8 (this is n, the limit)
d                      n,n
z;                    n,n^t
z$: n,n^t/t r Reverse stack -> n^t/t,n [ For loop that runs n times i1+ k d k,k z+ k,t+k$:       k/(t+k)
*      Multiply
]N    Close for loop and output as integer

As there is no ., it wraps around and starts over. However, n now produces -1 because the input is empty, which eventually leads to attempting to divide by 0, which halts the program.

# Julia, 141 bytes

z->(z-=1;a=90;c(k)=(k=big(k);(-1)^(k-1)/factorial(k-1)*(a-k)^(k-.5)*exp(a-k));(z+a)^(z+.5)*exp(-z-a)*(√(2π)+sum([c(k)/(z+k)for k=1:a-1])))

This creates an unnamed lambda function that accepts a real number and returns a real number. It uses Spounge's approximation to compute Gamma.

Ungolfed:

function Γ(z::Real)
# Spounge's approxmation is for Γ(z+1), so subtract 1
z -= 1

# Choose a number for the constant a, which determines the
# bound on the error
a = 90

# Define a function for the sequence c_k
function c(k::Integer)
# Convert k to a BigInt
k = big(k)
return (-1)^(k-1) / factorial(k-1) * (a-k)^(k-1/2) * exp(a-k)
end

# Compute the approximation
return (z+a)^(z+1/2) * exp(-z-a) * (√(2π) + sum([c(k)/(z+k) for k=1:a-1]))
end
• Very, very late golf, but z->(z-=1;a=90;c(k)=(k=big(k);(-1)^~-k/factorial(k-1)*(a-k)^(k-.5)*exp(a-k));(z+a)^(z+.5)*exp(-z-a)*(√(2π)+sum(c(k)/(z+k)for k=1:a-1))) should work for 137 bytes (at least in Julia 0.6) – Mr. Xcoder Mar 6 '18 at 19:07

# Japt, 45 bytes

Japt is a shortened version of JavaScript. Interpreter

$for(V=X=1;X<1e9;)$V*=(1+1/X pU /(1+U/X++;V/U

Of course, 1e9 = 1,000,000,000 iterations takes forever, so for testing, try replacing the 9 with a 6. (1e6 is accurate to ~5 significant figures. Using 1e8 on an input of 12 is enough to get the first six.)

Test-case results: (using 1e7 precision)

1:  1
-2.5: -0.9453083...
pi:  2.2880370...
-e: -0.9526812...
12:  39916536.5...
0.5:  1.7724538...
8.675309:  20248.319...
-10.1: -0.0000022...

### How it works

// Implicit: U = input number
$for( // Ordinary for loop. V=X=1; // Set V and X to 1. X<1e9;)$ //  Repeat while X is less than 1e9.
V*=      // Multiply V by:
(1+1/X   //  1 plus (1 over X),
pU /     //  to the power of U, divided by
(1+U/X++ //  1 plus (U over X). Increment X by 1.
;V/U     // Output the result of (V over U).

# TI-BASIC, 35 bytes

Input Z
1
For(I,1,ᴇ9
Ans(1+I⁻¹)^Z/(1+Z/I
End
Ans/Z

This uses the same algorithm as Zereges.

Caveat: I haven't actually tested this with the full 1e9 iterations; based on the time taken for smaller values, I expect the runtime to be on the order of months. However, it appears to converge, and there should be no problems with rounding errors. TI stores numbers as decimal floats with 14 digits of precision.

• You didn't test it?! – TanMath Nov 16 '15 at 2:25
• @TanMath I would, but I need my calculator for a final exam next month. – lirtosiast Nov 16 '15 at 2:48

# Python 3, 746878 73 bytes

Thanks @Mego and @xnor

This is a translation of the C++ answer by Zereges. Basically, this is an alternate definition of the gamma function, hence more accurate (and what is great is that is uses less bytes!)

I am sorry for all the mistakes!

def g(z,v=1):
for i in range(1,10**9):v*=(1+1/i)**z/(1+z/i)
return v/z
• The +1 on the range doesn't matter when you are dealing with billions. Also, you should specify that this is Python 3 - you'd need from __future__ import division for float division and a few terabytes of RAM to deal with the fact that range returns a list in Python 2. Plus, you can replace the 1.0s with 1s and shave off 4 bytes. – user45941 Nov 15 '15 at 5:21
• @TanMath: ^ is xor, didn't you mean ** as for exponentiation? – jermenkoo Nov 15 '15 at 5:34
• int(1e9) is just 10**9, and you don't need the parens around (1+1/i)**z. – xnor Nov 15 '15 at 8:32

# Python, 348448407390 389 bytes

Special thanks to @Mego!

A crossed out 448 is (almost) still a 448! :p

This is based on the Lanzcos approximation. Golfed from here

from cmath import*
C=[0.9999999999998099,676.5203681218851,-1259.1392167224028,771.3234287776531,-17‌6.6150291621406,12.507343278686905,-0.13857109526572012,9.984369578019572e-6,1.5‌​056327351493116e-7]
def g(z):
z-=1;if z.real<0.5:return pi/(sin(pi*z)*gamma(1-z))
else:
x=C[0]
for i in range(1,9):x+=C[i]/(z+i)
t=z+7.5;return sqrt(2*pi)*t**(z+0.5)*exp(-t)*x
• Please golf your submission by, at a minimum, removing whitespace (the spaces before and after -= and in import * for example) and using a one-character function name. Also note that you only need to support real input. – lirtosiast Nov 15 '15 at 2:47
• @ThomasKwa I have edited it. My original version did not work, here is a newer one. – TanMath Nov 15 '15 at 2:55
• @Mego edited... – TanMath Nov 15 '15 at 3:10
• This causes a recursion error - remove the z-=1; in the first line of gamma to fix it. You should also rename gamma to g for bytes saves and to avoid naming conflicts with cmath.gamma. Also drop the extraneous leading zeroes. – user45941 Nov 15 '15 at 3:40

# APL (Dyalog Unicode), 13 bytes

1e8∘(*÷⊣!+)÷⊢

Try it online!

The above code uses the ! function, but it is used as generalized binomial coefficient, not factorial. A version without the ! symbol altogether:

# APL (Dyalog Unicode), 20 bytes

1e8∘(*×(×/⊢÷+)∘⍳⍨)÷⊢

Try it online!

In both cases, 1e7 is used instead of 1e8 to show that the actual values calculated are close enough to the true values.

### How it works: the math

Uses the formula used in TanMath's Python answer and therefore Zereges' C++ answer.

Basically, the solutions by TanMath/Zereges compute the following, which I developed further to suit better to APL:

\begin{align} \Gamma(x) &\approx \frac1x \prod_{i=1}^{n}{\Bigl(1+\frac{1}{i}\Bigr)^x \div \Bigl(1+\frac{x}{i}\Bigr)} \\ &= \frac1x \prod_{i=1}^{n}{\Bigl(\frac{i+1}{i}\Bigr)^x \times \Bigl(\frac{i}{x+i}\Bigr)} \\ &= \frac1x \Bigl(\prod_{i=1}^{n}{\frac{i+1}{i}}\Bigr)^x \times \prod_{i=1}^{n}{\Bigl(\frac{i}{x+i}\Bigr)} \\ &\approx \frac{n^x}x \times \prod_{i=1}^{n}{\Bigl(\frac{i}{x+i}\Bigr)} \tag{1}\label{eq1}\\ &= \frac{n^x}x \div \binom{x+n}{n} \tag{2}\label{eq2}\\ \end{align}

The second solution uses equation $$\\eqref{eq1}\$$ while the first uses equation $$\\eqref{eq2}\$$, both using $$\n=10^8\$$ in order to get enough precision.

### How it works: the code

1e8∘(*÷⊣!+)÷⊢  ⍝ Input←x; n←1e8
1e8∘(    +)    ⍝ n+x
⊣!      ⍝ binom(n+x, n)
*÷        ⍝ n^x ÷ binom(n+x, n)
÷⊢  ⍝ Divide by x

1e8∘(*×(×/⊢÷+)∘⍳⍨)÷⊢  ⍝ Input←x; n←1e8
1e8∘(  (     )∘⍳⍨)    ⍝ Generate 1..n; for each i in 1..n,
⊢÷+         ⍝ i÷(x+i)
×/            ⍝ Product over all i
*×               ⍝ Multiply n^x
÷⊢  ⍝ Divide by x

# Julia, 41 bytes

x->prod([(1+1/i)^x/(1+x/i)for i=1:1E7])/x

This is a translation of Zereges' C++ answer. While my other Julia answer finishes instantaneously, this is rather slow. It computes the test cases in a couple seconds each on my computer.

Ungolfed:

function f(x::Real)
prod([(1 + 1/i)^x / (1 + x/i) for i = 1:1E7]) / x
end

# Prolog, 114 bytes

This is a translation of Zereges' C++ answer.

q(F,F,V,Z):-X is V/Z,write(X).
q(F,T,V,Z):-W is(1+1/F)**Z/(1+Z/F)*V,I is F+1,q(I,T,W,Z).
p(N):-q(1.0,1e9,1,N),!.

Try it out online here
Run it with a query of the form:

p(12).

Running it with 1e9 recursions takes about 15 minutes.
If you decrease it to 1e6 it takes about 1 second which makes for easier (but less acurate) testing.
Running it in an interpreter on you computer/laptop is most likely faster for most people as well.

# Mathematica, 40 bytes

NProduct[(1+1/n)^#/(1+#/n),{n,1,∞}]/#&