function y = gampdf_scalar(x, a, b, outProto)
% GAMPDF Gamma probability density function for one element, optimized for GPU.

%   Copyright 2019 The MathWorks, Inc.

% Initialize Y to zero
y = zeros("like", outProto);

% Return NaN for out of range parameters.
if (a < 0 || b <= 0 || isnan(a) || isnan(b) || isnan(x))
    y = NaN("like", outProto);
end

% Scale
z = x./b;

% Special cases
if (z==0 && a==1 && b>0) %#ok<BDSCI>
    y = cast(1/b, "like", outProto);
elseif (z==0 && a<1 && b>0) %#ok<BDSCI>
    y = Inf("like", outProto);
elseif (z>0 && z<Inf && a>0 && a<Inf && b>0) %#ok<BDSCI>
    % Normal cases
    a = a-1;
    if a<0
        y = f(z,a+1) .* exp(log(a+1)-log(z)) ./ b;
    else
        y = f(z,a) ./ b;
    end
end

end


function y = f(z,a)
% Compute gampdf without error checking for z>0 and a>0.
if a<=realmin*z
    % z term dominates
    y = exp(-z);
elseif z<realmin*a
    % Normal expansion through logs
    y = exp( a.*log(z) -z -gammaln(a+1) );
else
    % Loader's saddle point expansion
    lnsr2pi = 0.9189385332046727; % log(sqrt(2*pi))
    y = exp(-lnsr2pi -0.5*log(a) - stirlerr_scalar(a) - binodeviance_scalar(a,z));
end
end
