function y = hygepdf_scalar(x, m, k, n, outProto)
% UNIFPDF Uniform probability density function for one element, optimised for GPU.

%   Copyright 2019 The MathWorks, Inc.

% NaN in, NaN out.
if isnan(x) || isnan(m) || isnan(k) || isnan(n)
    y = NaN("like", outProto);
    return
end

% Return NaN for values of the parameters outside their respective limits.
if (m < 0 || k < 0 || n < 0 || n > m || k > m ...
        || round(m) ~= m || round(k) ~= k || round(n) ~= n)
    y = NaN("like", outProto);
    return
end

% Values of X for which Y is zero by inspection.
if (m - k - n + x + 1 <= 0 || x < 0 || round(x) ~= x || x > n || x > k)
    y = zeros("like", outProto);
    return
end

% At n=0 or n==m must have pdf=1
if (n==0 || n==m)
    y = ones("like", outProto);
    return
end

% Remaining integer values of x are within the correct range
lnsqr2pi = 0.9189385332046727;

% Get the biggest term for m-choose-n.
% Use binomial pdf with p=n/m. Since p=n/m, binomial deviance is zero.
p = n./m;
mn = -lnsqr2pi -0.5*log(n.*(1-n./m)) ...
    +stirlerr_scalar(m) -stirlerr_scalar(n) -stirlerr_scalar(m-n);

if x==0
    kx = cast(k.*log(1-p), "like", outProto);
elseif x==k
    kx = cast(x.*log(p), "like", outProto);
else
    kx = -lnsqr2pi -0.5*log(x.*(1-x./k)) ...
        +stirlerr_scalar(k) -stirlerr_scalar(x) -stirlerr_scalar(k-x) ...
        -binodeviance_scalar(x,k.*p) ...
        -binodeviance_scalar(k-x,k.*(1-p));
end

if n==x
    mknx = cast((m-k).*log(1-p), "like", outProto);
elseif (n-x)==(m-k)
    mknx = cast((n-x).*log(p), "like", outProto);
else
    mknx = -lnsqr2pi -0.5*log((n-x).*(1-(n-x)./(m-k))) ...
        +stirlerr_scalar(m-k) -stirlerr_scalar(n-x) -stirlerr_scalar(m-k-n+x) ...
        -binodeviance_scalar(n-x,(m-k).*p) ...
        -binodeviance_scalar(m-k-n+x,(m-k).*(1-p));
end

y = exp(kx + mknx - mn);
end
