function b = gammainc(x,a,tail)
%MATLAB Code Generation Library Function

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

narginchk(2,3);
coder.internal.assert(isreal(x) && isreal(a) && ...
    isfloat(x) && isfloat(a) && ...
    ~issparse(x) && ~issparse(a), ...
    'MATLAB:gammainc:inputType');
coder.internal.errorIf( ...
    coder.internal.vAllOrAny('any',a,@(t)t < zeros('like',t)), ...
    'MATLAB:gammainc:NegativeArg');
if nargin < 3
    upper = false;
    scaled = false;
else
    coder.internal.prefer_const(tail);
    isUpperStr = strcmp(tail,'upper');
    isScaledUpperStr = strcmp(tail,'scaledupper');
    isLowerStr = strcmp(tail,'lower');
    isScaledLowerStr = strcmp(tail,'scaledlower');
    coder.internal.assert(isUpperStr || isScaledUpperStr || ...
        isLowerStr || isScaledLowerStr, ...
        'MATLAB:gammainc:InvalidTailArg');
    upper = isUpperStr || isScaledUpperStr;
    scaled = isScaledUpperStr || isScaledLowerStr;
end
zero = complex(coder.internal.scalarEg(x,a));
b = coder.internal.scalexpAlloc(zero,x,a);
SCALARX = coder.const(eml_is_const(isscalar(x)) && isscalar(x));
if SCALARX
    x1 = double(x(1));
end
SCALARA = coder.const(eml_is_const(isscalar(a)) && isscalar(a));
if SCALARA
    a1 = double(a(1));
    la1 = log(a1);
    lgap11 = gammaln(a1 + 1);
end
if SCALARA && SCALARX
    b(1) = scalar_gammainc(x1,a1,la1,lgap11,upper,scaled);
    return
elseif isempty(b)
    return
end
n = coder.internal.indexInt(eml_numel(b));
if coder.internal.useParforConst(mfilename,n)
    parfor k = 1:n
        if SCALARX
            xk = x1;
        else
            xk = double(x(k));
        end
        if SCALARA
            ak = a1;
            la = la1;
            lgap1 = lgap11;
        else
            ak = double(a(k));
            la = log(ak);
            lgap1 = gammaln(ak + 1);
        end
        b(k) = scalar_gammainc(xk,ak,la,lgap1,upper,scaled);
    end
else
    for k = 1:n
        if SCALARX
            xk = x1;
        else
            xk = double(x(k));
        end
        if SCALARA
            ak = a1;
            la = la1;
            lgap1 = lgap11;
        else
            ak = double(a(k));
            la = log(ak);
            lgap1 = gammaln(ak + 1);
        end
        b(k) = scalar_gammainc(xk,ak,la,lgap1,upper,scaled);
    end
end

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

function b = scalar_gammainc(x,a,la,lgap1,upper,scaled)
b = complex(coder.internal.scalarEg(x,a));
if ~(x < 0) % includes NaN
    b(1) = eml_gammainc(x,a,la,lgap1,upper,scaled);
else
    [yr,yi] = gammaincNegX(x,a,lgap1,upper,scaled);
    b(1) = complex(yr,yi);
end

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

function [gami_r,gami_i] = gammaincNegX(x,a,lgap1,upper,scaled)
% Compute the incomplete gamma function, lower or upper tail, x < 0.
% Based on Knuesel, L. (1986) 'Computation of the Chi-square and Poisson
% Distribution', SIAM J. Sci. Stat. Comput., 7(3):1022-1036.
coder.internal.prefer_const(upper,scaled);
if upper
    tail = 1;
else
    tail = 0;
end
tol = eps;
% Check for trivial and edge cases
if x == 0
    if scaled
        if upper
            gami_r = coder.internal.inf;
        else
            gami_r = 0;
        end
        gami_i = 0;
    else
        gami_r = tail;
        gami_i = 0;
    end
    return
elseif ~(x < 0) || x == -coder.internal.inf || ...
        ~(a > 0) || a == coder.internal.inf
    gami_r = coder.internal.nan;
    gami_i = coder.internal.nan;
    return
end
% Compute log(p(a,x)*x/a) = log(exp(-x)*x^a/gamma(a+1))
% pax = exp(-x + a*log(x) - lgap1);
pax_modulus = exp(-x + a*log(-x) - lgap1);
pax_angle = a;
if a == 1
    pax_r = -pax_modulus;
    pax_i = 0;
else
    [cpax,spax] = coder.internal.scalar.cospiAndSinpi(pax_angle);
    pax_r = pax_modulus*cpax;
    pax_i = pax_modulus*spax;
end
% Compute the lower tail using the alternating series.
Iax = 1;
term = 1;
i = 1;
while true
    b = a + i; % from scratch to avoid accumulating roundoff
    term = term*x/b;
    Iax = Iax + term;
    if abs(term) <= tol*abs(Iax)
        break % <= is important, in the case when Iax and term overflow
    end
    i = i + 1;
end
% For x = -r, I(a,x) = ...
Iax = abs(Iax);
if scaled
    % This has the limit 1 as x->0, for any finite a>0. ...
    gami_r = Iax;
    gami_i = 0;
    if upper
        % When a is an integer, this becomes Inf + 0*i, the same ...
        if a == 1
            paxinv_r = -1/pax_modulus ;
            paxinv_i =  0.0;
        else
            [cpax,spax] = coder.internal.scalar.cospiAndSinpi(pax_angle);
            paxinv_r = cpax/pax_modulus;
            paxinv_i = -spax/pax_modulus;
        end
        gami_r = paxinv_r - gami_r;
        gami_i = paxinv_i - gami_i;
    end
else  % gami = pax * Iax;
    gami_r = Iax*pax_r;
    if pax_i == 0
        gami_i = 0;
    else
        gami_i = Iax*pax_i;
    end
    if upper
        % gami = 1 - pax * Iax;
        gami_r = 1 - gami_r;
        gami_i = 0 - gami_i;
    end
end

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