function gamma( x ) {
	var result;
	var lzgh;
	var zgh;
	var hp;

	if ( x === NINF || isnan( x ) ) {
		return NaN;
	}
	if ( x === 0.0 ) {
		if ( isNegativeZero( x ) ) {
			return NINF;
		}
		return PINF;
	}
	result = 1.0;
	if ( x < 0.0 ) {
		if ( isInteger( x ) ) {
			return NaN;
		}
		if ( x <= -20.0 ) {
			result = gamma( -x ) * sinpx( x );
			if ( abs(result) < 1.0 && MAX_VALUE * abs(result) < PI ) {
				return ( signum( result ) === 1 ) ? NINF : PINF;
			}
			result = -PI / result;
			return result;
		}
		// Shift x to > 1:
		while ( x < 0.0 ) {
			result /= x;
			x += 1.0;
		}
	}
	if ( floor( x ) === x && x < MAX_FACTORIAL ) {
		result *= FACTORIALS[ floor( x ) - 1 ];
	}
	else if ( x < SQRT_EPSILON ) {
		if ( x < 1.0 / MAX_VALUE ) {
			result = PINF;
		}
		result *= ( 1.0 / x ) - EULERGAMMA;
	}
	else {
		result *= gammaLanczosSum( x );
		zgh = ( x + LANCZOS_G - HALF );
		lzgh = ln( zgh );
		if ( x * lzgh > MAX_LN ) {
			// We're going to overflow unless this is done with care:
			if ( lzgh * x / 2.0 > MAX_LN ) {
				return ( signum( result ) === 1 ) ? PINF : NINF;
			}
			hp = pow( zgh, ( x / 2.0 ) - 0.25 );
			result *= hp / exp( zgh );
			if ( MAX_VALUE / hp < result ) {
				return ( signum( result ) === 1 ) ? PINF : NINF;
			}
			result *= hp;
		} else {
			result *= pow( zgh, x - HALF ) / exp( zgh );
		}
	}
	return result;
}
Evaluates the gamma function.

@param {number} x - input value
@returns {number} function value

@example
var v = gamma( 4.0 );
// returns 6.0

@example
var v = gamma( -1.5 );
// returns ~2.363

@example
var v = gamma( -0.5 );
// returns ~-3.545

@example
var v = gamma( 0.5 );
// returns ~1.772

@example
var v = gamma( 0.0 );
// returns Infinity

@example
var v = gamma( -0.0 );
// returns -Infinity

@example
var v = gamma( NaN );
// returns NaN