function x = gamma(x)
%MATLAB Code Generation Library Function

%   Copyright 1984-2018 The MathWorks, Inc.
%#codegen

eml_invariant(nargin == 1,'MATLAB:minrhs');
eml_invariant(isfloat(x),'Coder:toolbox:gamma_2');
eml_invariant(isreal(x) && ~issparse(x),'MATLAB:gamma:notFullReal');
x = coder.internal.applyScalarFunctionInPlace(mfilename,@scalar_gamma,x);

%--------------------------------------------------------------------------

function x = scalar_gamma(x)
gam = cast([1.0; 1.0; 2.0; 6.0; 24.0; 120.0; 720.0; 5040.0; 40320.0;
    362880.0; 3628800.0; 39916800.0; 479001600.0; 6227020800.0;
    87178291200.0; 1307674368000.0; 20922789888000.0; 355687428096000.0;
    6402373705728000.0; 121645100408832000.0; 2432902008176640000.0;
    51090942171709440000.0; 1124000727777607680000.0], ...
    'like',real(x));
p = cast([-1.71618513886549492533811e+0; 2.47656508055759199108314e+1;
    -3.79804256470945635097577e+2; 6.29331155312818442661052e+2;
    8.66966202790413211295064e+2; -3.14512729688483675254357e+4;
    -3.61444134186911729807069e+4; 6.64561438202405440627855e+4], ...
    'like',real(x));
q = cast([-3.08402300119738975254353e+1; 3.15350626979604161529144e+2;
    -1.01515636749021914166146e+3; -3.10777167157231109440444e+3;
    2.25381184209801510330112e+4; 4.75584627752788110767815e+3;
    -1.34659959864969306392456e+5; -1.15132259675553483497211e+5], ...
    'like',real(x));
c = cast([-1.910444077728e-03; 8.4171387781295e-04;
    -5.952379913043012e-04; 7.93650793500350248e-04;
    -2.777777777777681622553e-03; 8.333333333333333331554247e-02;
    5.7083835261e-03], ...
    'like',real(x));
if (x >= 1) && (x <= 23) && (x == floor(x))
    x = cast(gam(x),'like',x);
elseif (x < 1) && (x == floor(x)) %x is a non-positive integer
    x = coder.internal.inf('like',x);
elseif isnan(x)
    x = coder.internal.nan('like',x);
elseif isinf(x)
    x = coder.internal.inf('like',x);
else
    fact = ones('like',x);
    n = zeros('like',x);
    parity = false;
    if x <= 0  % Catch negative x
        y = -x;
        yint = floor(y);
        parity = yint ~= floor(y/2)*2;
        r = y - yint;
        fact = -pi/sinpi(r);
        x = y + 1;
    end
    if x < 12
        xkold = x;
        if x < 1
            z = x;
            x = x + 1;
        else
            n = floor(x) - 1;
            x = x - n;
            z = x - 1;
        end
        xnum = 0*z;
        xden = xnum + 1;
        for i = 1:8
            xnum = (xnum + p(i))*z;
            xden = xden*z + q(i);
        end
        r = xnum/xden + 1;
        if xkold < x
            r = r/xkold;
        elseif xkold > x
            for j = 1:n
                r = r*x;
                x = x + 1;
            end
        end
    else
        y = x;
        ysq = y*y;
        sum = cast(c(7),'like',x);
        for i = 1:6
            sum = sum/ysq + c(i);
        end
        spi = cast(0.9189385332046727417803297,'like',real(x));
        sum = sum/y - y + spi;
        sum = sum + (y - 0.5)*eml_log(y);
        r = eml_exp(sum);
    end
    if parity
        r = -r;
    end
    if fact ~= 1
        r = fact/r;
    end
    x = r;
end

%--------------------------------------------------------------------------
