How to make a function that computes the factorial for numbers with decimals?

Question 1

I might have found an existing solution... It's an implementation of Lanczos method, I found it at the swedish wikipedia (http://sv.wikipedia.org/wiki/Gammafunktionen). It was written in python and says to be correct up to 15 decimals. I ported it to js, cross checked some random values against (http://www.efunda.com/math/gamma/findgamma.cfm).

http://jsfiddle.net/Fzy9C/

var g = 7;
var C = [0.99999999999980993, 676.5203681218851, -1259.1392167224028,771.32342877765313, -176.61502916214059, 12.507343278686905, -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7];

function gamma(z) {

    if (z < 0.5) return Math.PI / (Math.sin(Math.PI * z) * gamma(1 - z));
    else {
        z -= 1;

        var x = C[0];
        for (var i = 1; i < g + 2; i++)
        x += C[i] / (z + i);

        var t = z + g + 0.5;
        return Math.sqrt(2 * Math.PI) * Math.pow(t, (z + 0.5)) * Math.exp(-t) * x;
    }
}

(and ofcourse it does not support imaginary numbers, since js does not)

Question 2

As an alternative to the other answers here, here's a much simpler approximation for the gamma function, proposed in 2007 by Gergő Nemes. (See the wikipedia page on Stirling's approximation).

enter image description here

This can be implemented directly in JavaScript in a single line:

function gamma(z) {
  return Math.sqrt(2 * Math.PI / z) * Math.pow((1 / Math.E) * (z + 1 / (12 * z - 1 / (10 * z))), z);
}

You can see this in action on jsFiddle.

This is accurate to 8 digits for z > 8, but it is still accurate to a handful of digits for smaller z. It is not quite as accurate as Lanczos approximation, but it is simpler and also slightly faster.

Note that the gamma function and the factorial function are slightly different. The factorial function can be defined in terms of the gamma function thus:

function factorial(n) {
  return gamma(n + 1);
}

Question 3

This is not a trivial problem. There is not a simple closed-form formula for the gamma function. That said, there are some numerical approximations that should suit your needs.

The following answer will be using a technique called Lanczos approximation. The formula is as follows:

enter image description here

where g is an arbitrarily chosen constant that controls how accurate the approximation will be. For larger g, the approximation will be more accurate. A_g(z) is defined thus:

enter image description here

The hardest part is finding A_g(z), since p_n is also defined with a complicated formula dependent on g.

I can't take too much credit for the following code, since I am just writing a port of the Python program on the wikipedia page.

function gamma(n) {  // accurate to about 15 decimal places
  //some magic constants 
  var g = 7, // g represents the precision desired, p is the values of p[i] to plug into Lanczos' formula
      p = [0.99999999999980993, 676.5203681218851, -1259.1392167224028, 771.32342877765313, -176.61502916214059, 12.507343278686905, -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7];
  if(n < 0.5) {
    return Math.PI / Math.sin(n * Math.PI) / gamma(1 - n);
  }
  else {
    n--;
    var x = p[0];
    for(var i = 1; i < g + 2; i++) {
      x += p[i] / (n + i);
    }
    var t = n + g + 0.5;
    return Math.sqrt(2 * Math.PI) * Math.pow(t, (n + 0.5)) * Math.exp(-t) * x;
  }
}

and of course, by definition of the gamma function:

function factorial(n) {
  return gamma(n + 1);
}

You can see this in action on jsFiddle.

Question 4

Just to complete @apelsinapa answer to correct the calculation for an integer (we didn't get an integer solution when inputing an integer number).

@apelsinapa's great solution:

var g = 7;
var C = [0.99999999999980993, 676.5203681218851, -1259.1392167224028,771.32342877765313, -176.61502916214059, 12.507343278686905, -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7];

function gamma(z) {
    if (z < 0.5) return Math.PI / (Math.sin(Math.PI * z) * gamma(1 - z));
    else {
        z -= 1;

        var x = C[0];
        for (var i = 1; i < g + 2; i++)
        x += C[i] / (z + i);

        var t = z + g + 0.5;
        return Math.sqrt(2 * Math.PI) * Math.pow(t, (z + 0.5)) * Math.exp(-t) * x;
    }
}

And to get a correct answer for integer:

 function factorialOfNumber(number) {
    if (number % 1 != 0 || number<0){
        return gamma(number + 1);
    }
    else {
        if(number == 0) {
           return 1;
        }
        for(var i = number; --i; ) {
           number *= i;
        }
        return number;
    }
}

Question 5

Here's a version I wrote a few years ago ... a bit messy but tested :)

var
    M_GAMMA = [76.18009172947146, -86.50532032941677, 24.01409824083091, -1.231739572450155, 0.1208650973866179e-2, -0.5395239384953e-5],
    M_GAMMAS = 6;

function gammaFn(x) // Modified to JavaScript from "Numerical Recipies in C" by P. Mainwaring
{
    var i = 0, n = ++x, tmp = x + 5.5, ser = 1.000000000190015;
    for (tmp -= (x + 0.5) * Math.log(tmp); i < M_GAMMAS; ++i) ser += M_GAMMA [i] / ++n;
    return Math.log(2.5066282746310005 * ser / x) - tmp;
}

function fact(x) { return x > 1 ? Math.exp(gammaFn(x)) : 1 }

function combin(n, k) { return (Math.exp(gammaFn(n) - gammaFn(n - k) - gammaFn(k)) + 0.5) | 0 } // Ms Excel COMBIN() n! / k!(n - k)!

n = 49; k = 6; alert(fact(n) + ' ' + fact(k) + ' ' + combin(n, k)); // Lottery odds! (13,983,816)

The gamma and gammaLn functions are then:

function gammaLn(x) { return gammaFn(--x) }

function gamma(x) { return Math.exp(gammaLn(x)) }

:-)

Question 6

If you're just looking for the function to compute factorials of real numbers then you only need this bit of code from Lanczos approximation:

function = factorial(z) {

  var g = 7;
  var C = [0.99999999999980993, 676.5203681218851, -1259.1392167224028, 771.32342877765313, -176.61502916214059, 12.507343278686905, -0.13857109526572012, 9.9843695780195716e-6, 1.5056327351493116e-7];
  var x = C[0];

  for (var i = 1; i < g + 2; i++)
  x += C[i] / (z + i);

  var t = z + g + 0.5;
  return Math.sqrt(2 * Math.PI) * Math.pow(t, (z + 0.5)) * Math.exp(-t) * x;
}

Works great for negative numbers in addition to decimals.