23
\$\begingroup\$

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:

gamma function definition

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

Leaderboard

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.

To make sure that your answer shows up, please start your answer with a headline, using the following Markdown template:

## 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

<style>body { text-align: left !important} #answer-list { padding: 10px; width: 290px; float: left; } #language-list { padding: 10px; width: 290px; float: left; } table thead { font-weight: bold; } table td { padding: 5px; }</style><script src="https://ajax.googleapis.com/ajax/libs/jquery/2.1.1/jquery.min.js"></script> <link rel="stylesheet" type="text/css" href="//cdn.sstatic.net/codegolf/all.css?v=83c949450c8b"> <div id="language-list"> <h2>Shortest Solution by Language</h2> <table class="language-list"> <thead> <tr><td>Language</td><td>User</td><td>Score</td></tr> </thead> <tbody id="languages"> </tbody> </table> </div> <div id="answer-list"> <h2>Leaderboard</h2> <table class="answer-list"> <thead> <tr><td></td><td>Author</td><td>Language</td><td>Size</td></tr> </thead> <tbody id="answers"> </tbody> </table> </div> <table style="display: none"> <tbody id="answer-template"> <tr><td>{{PLACE}}</td><td>{{NAME}}</td><td>{{LANGUAGE}}</td><td>{{SIZE}}</td><td><a href="{{LINK}}">Link</a></td></tr> </tbody> </table> <table style="display: none"> <tbody id="language-template"> <tr><td>{{LANGUAGE}}</td><td>{{NAME}}</td><td>{{SIZE}}</td><td><a href="{{LINK}}">Link</a></td></tr> </tbody> </table><script>var QUESTION_ID = 63887; var ANSWER_FILTER = "!t)IWYnsLAZle2tQ3KqrVveCRJfxcRLe"; var COMMENT_FILTER = "!)Q2B_A2kjfAiU78X(md6BoYk"; var OVERRIDE_USER = 45941; var answers = [], answers_hash, answer_ids, answer_page = 1, more_answers = true, comment_page; function answersUrl(index) { return "https://api.stackexchange.com/2.2/questions/" + QUESTION_ID + "/answers?page=" + index + "&pagesize=100&order=desc&sort=creation&site=codegolf&filter=" + ANSWER_FILTER; } function commentUrl(index, answers) { return "https://api.stackexchange.com/2.2/answers/" + answers.join(';') + "/comments?page=" + index + "&pagesize=100&order=desc&sort=creation&site=codegolf&filter=" + COMMENT_FILTER; } function getAnswers() { jQuery.ajax({ url: answersUrl(answer_page++), method: "get", dataType: "jsonp", crossDomain: true, success: function (data) { answers.push.apply(answers, data.items); answers_hash = []; answer_ids = []; data.items.forEach(function(a) { a.comments = []; var id = +a.share_link.match(/\d+/); answer_ids.push(id); answers_hash[id] = a; }); if (!data.has_more) more_answers = false; comment_page = 1; getComments(); } }); } function getComments() { jQuery.ajax({ url: commentUrl(comment_page++, answer_ids), method: "get", dataType: "jsonp", crossDomain: true, success: function (data) { data.items.forEach(function(c) { if (c.owner.user_id === OVERRIDE_USER) answers_hash[c.post_id].comments.push(c); }); if (data.has_more) getComments(); else if (more_answers) getAnswers(); else process(); } }); } getAnswers(); var SCORE_REG = /<h\d>\s*([^\n,<]*(?:<(?:[^\n>]*>[^\n<]*<\/[^\n>]*>)[^\n,<]*)*),.*?(\d+)(?=[^\n\d<>]*(?:<(?:s>[^\n<>]*<\/s>|[^\n<>]+>)[^\n\d<>]*)*<\/h\d>)/; var OVERRIDE_REG = /^Override\s*header:\s*/i; function getAuthorName(a) { return a.owner.display_name; } function process() { var valid = []; answers.forEach(function(a) { var body = a.body; a.comments.forEach(function(c) { if(OVERRIDE_REG.test(c.body)) body = '<h1>' + c.body.replace(OVERRIDE_REG, '') + '</h1>'; }); var match = body.match(SCORE_REG); if (match) valid.push({ user: getAuthorName(a), size: +match[2], language: match[1], link: a.share_link, }); else console.log(body); }); valid.sort(function (a, b) { var aB = a.size, bB = b.size; return aB - bB }); var languages = {}; var place = 1; var lastSize = null; var lastPlace = 1; valid.forEach(function (a) { if (a.size != lastSize) lastPlace = place; lastSize = a.size; ++place; var answer = jQuery("#answer-template").html(); answer = answer.replace("{{PLACE}}", lastPlace + ".") .replace("{{NAME}}", a.user) .replace("{{LANGUAGE}}", a.language) .replace("{{SIZE}}", a.size) .replace("{{LINK}}", a.link); answer = jQuery(answer); jQuery("#answers").append(answer); var lang = a.language; lang = jQuery('<a>'+lang+'</a>').text(); languages[lang] = languages[lang] || {lang: a.language, lang_raw: lang.toLowerCase(), user: a.user, size: a.size, link: a.link}; }); var langs = []; for (var lang in languages) if (languages.hasOwnProperty(lang)) langs.push(languages[lang]); langs.sort(function (a, b) { if (a.lang_raw > b.lang_raw) return 1; if (a.lang_raw < b.lang_raw) return -1; return 0; }); for (var i = 0; i < langs.length; ++i) { var language = jQuery("#language-template").html(); var lang = langs[i]; language = language.replace("{{LANGUAGE}}", lang.lang) .replace("{{NAME}}", lang.user) .replace("{{SIZE}}", lang.size) .replace("{{LINK}}", lang.link); language = jQuery(language); jQuery("#languages").append(language); } }</script>

\$\endgroup\$
3
  • 2
    \$\begingroup\$ Please provide explicit bounds for t such that |gamma(t)|<10^9 \$\endgroup\$
    – flawr
    Nov 14, 2015 at 23:49
  • \$\begingroup\$ The link is not a reference implementation,... \$\endgroup\$
    – sergiol
    Mar 9, 2018 at 0:03
  • \$\begingroup\$ @sergiol Reworded it \$\endgroup\$
    – user45941
    Mar 9, 2018 at 0:04

12 Answers 12

9
\$\begingroup\$

C++14, 86 85 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 an 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
    }
\$\endgroup\$
9
  • \$\begingroup\$ @Mego Of course it is! Thanks. \$\endgroup\$
    – Zereges
    Nov 14, 2015 at 23:38
  • \$\begingroup\$ 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. \$\endgroup\$
    – flawr
    Nov 14, 2015 at 23:40
  • \$\begingroup\$ @Mego Microsoft compiler does not need neither of theese includes. \$\endgroup\$
    – Zereges
    Nov 14, 2015 at 23:46
  • \$\begingroup\$ @Mego Microsoft (R) C/C++ Optimizing Compiler Version 19.00.23026 for x86 \$\endgroup\$
    – Zereges
    Nov 14, 2015 at 23:49
  • \$\begingroup\$ @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 \$\endgroup\$
    – Zereges
    Nov 14, 2015 at 23:59
6
\$\begingroup\$

Minkolang 0.12, 35 34 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.

Euler's alternative formulation of the gamma function

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.

exploded version of the infinite product

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.

\$\endgroup\$
4
\$\begingroup\$

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
\$\endgroup\$
1
  • \$\begingroup\$ 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) \$\endgroup\$
    – Mr. Xcoder
    Mar 6, 2018 at 19:07
3
\$\begingroup\$

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.

\$\endgroup\$
3
\$\begingroup\$

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).
\$\endgroup\$
3
\$\begingroup\$

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.

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

Python 3, 74 68 78 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
\$\endgroup\$
3
  • 1
    \$\begingroup\$ 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. \$\endgroup\$
    – user45941
    Nov 15, 2015 at 5:21
  • 2
    \$\begingroup\$ @TanMath: ^ is xor, didn't you mean ** as for exponentiation? \$\endgroup\$
    – jermenkoo
    Nov 15, 2015 at 5:34
  • 3
    \$\begingroup\$ int(1e9) is just 10**9, and you don't need the parens around (1+1/i)**z. \$\endgroup\$
    – xnor
    Nov 15, 2015 at 8:32
3
\$\begingroup\$

Python, 348 448 407 390 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
\$\endgroup\$
4
  • 1
    \$\begingroup\$ 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. \$\endgroup\$
    – lirtosiast
    Nov 15, 2015 at 2:47
  • \$\begingroup\$ @ThomasKwa I have edited it. My original version did not work, here is a newer one. \$\endgroup\$
    – TanMath
    Nov 15, 2015 at 2:55
  • \$\begingroup\$ @Mego edited... \$\endgroup\$
    – TanMath
    Nov 15, 2015 at 3:10
  • \$\begingroup\$ 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. \$\endgroup\$
    – user45941
    Nov 15, 2015 at 3:40
2
+200
\$\begingroup\$

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
\$\endgroup\$
1
\$\begingroup\$

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
\$\endgroup\$
1
\$\begingroup\$

Mathematica, 40 bytes

NProduct[(1+1/n)^#/(1+#/n),{n,1,∞}]/#&
\$\endgroup\$
0
\$\begingroup\$

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.

\$\endgroup\$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.