function p = hygecdf_scalar(x, m, k, n, outProto)
% HYGECDFPDF Hypergeometric cumulative distribution 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)
    p = NaN("like", outProto);
    return
end

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

% Handle values of X for which P is zero by inspection.
if (m - k - n + x + 1 <= 0 || x < 0)
    p = zeros("like", outProto);
    return
end

% Handle values of X for which P is one by inspection.
if (x >= n || x >= k)
    p = ones("like", outProto);
    return
end


% Compute p when xx >= 0.

% For accuracy, for x values that are larger than the mean, compute the
% upper tail 1-p instead of the lower tail p.  This will tend to
% increase the number of significant bits in the result.
if (x <= k*n/m)
    p = F(floor(x),m,k,n,outProto);
else
    p = 1 - F(n-floor(x)-1,m,m-k,n,outProto);
end
end


% -----------------------------------------
% Hypergeometric cdf without error checking
function p=F(x,m,k,n,outProto)
dens = hygepdf_scalar(x,m,k,n,outProto);
% Compute hygecdf(x,m,k,n)/hygepdf(x,m,k,n) with a series whose terms can
% be computed recursively.
%
% For a vector x = [a b c] we want to calculate sum(cumprod(x,'reverse')),
% i.e. (a*b*c) + (b*c) + (c). Note that the first terms are small compared
% to the last so it is important to add from beginning to end for minimum
% rounding error. We therefore expand to c*(1+b*(1+a)), which maintains
% numerical stability without requiring two loops over the data.
ratio = zeros("like", outProto);
for y=0:x-1
    thisTerm = ((y+1) * (m-k-n+y+1)) / ((n-y) * (k-y));
    ratio = thisTerm*(1 + ratio);
end
ratio = 1 + ratio;

p = ratio*dens;

% Make sure that round-off errors never make P greater than 1.
if p > 1
    p = ones("like", outProto);
end
end
