function y = binocdf_scalar(x, n, p, outType)
% BINOCDF Binomial cumulative probability density function for one element, optimised for GPU.

%   Copyright 2019 The MathWorks, Inc.

% Deal with special cases for a single element

smallp = 1e-4;  % to maintain accuracy, don't switch tails below here
bign = 1e5;     % to conserve memory, sum over smaller range above here

if (x >= n)
    % Y = 1 if X >= N
    y = ones("like", outType);
    
elseif (x < 0)
    % Y = 0 if X < 0
    y = zeros("like", outType);
    
elseif isnan(x)
    % Y = NaN if X is NaN
    y = NaN("like", outType);
    
elseif (n < 0 || p < 0 || p > 1 || round(n) ~= n)
    % Return NaN if any arguments are outside of their respective limits.
    y = NaN("like", outType);
    
elseif (p == 0)
    % assign 1 to p==0 indices.
    y = ones("like", outType);
    
elseif (p == 1)
    % assign 0 to p==1 indices
    y = zeros("like", outType);
    
else
    % Sum the PDF directly
    xx = floor(x);
    if (n<bign) && ((x <= n*p) || (p<smallp))
        % For the simple case of small(ish) x, sum the PDF from 0 to n
        y = sumfrom0(xx, n, p, outType);
    elseif (n<bign)
        % Get more accuracy by summing the PDF from the other direction in
        % the upper tail
        y = 1 - sumfrom0(n-xx-1, n, 1-p, outType);
    else
        % Do the full sum from A-to-B
        y = sumA2B(xx, n, p, outType);
    end
    
end

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

end


function y = sumfrom0(xx, n, p, outType)
% sum PDF series starting from 0

% In the original this is cumsum( binopdf(0:xx,n,p) )
y = zeros("like", outType);
for ii=0:xx
    y = y + binopdf_scalar(ii, n, p, outType);
end
end


function y = sumA2B(x, N, p, outType)
% sum series from A to B, with these selected to avoid negligible values

% See how far we need to look into each tail
mu = N*p;
std = sqrt(N*p*(1-p));
t1=40;
t2=10;
while binopdf_scalar(floor(mu-t1*std), N, p, outType) > eps(0)
    t1 = 1.5*t1;
end
while binopdf_scalar(ceil(mu+t2*std), N, p, outType) > eps(1)
    t2 = 1.5*t2;
end

% find the limits for the sum
a = max(0, floor(mu-t1*std));  % lower limit of sum
b = ceil(mu+t2*std);           % upper limit of sum

% outside limits, set to known values
y = zeros("like", outType);
if x<a
    % Leave at zero
elseif x>b
    y = ones("like", outType);
else
    % Sum from a to x
    for ii=a:x
        y = y + binopdf_scalar(ii, N, p, outType);
    end
end

end
