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

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

eml_invariant(nargin >= 3,'MATLAB:minrhs');
eml_invariant(isfloat(x) && isfloat(a) && isfloat(b), ...
    'Coder:toolbox:betaincinv_2');
eml_invariant(isreal(x) && isreal(a) && isreal(b), ...
    'Coder:toolbox:betaincinv_3');
if nargin < 4
    upper = false;
else
    coder.internal.prefer_const(tail);
    upper = strcmp(tail, 'upper');
    eml_invariant(upper || strcmp(tail,'lower'), ...
        'Coder:toolbox:betaincinv_4');
end
y = coder.internal.scalexpAlloc( ...
    complex(coder.internal.scalarEg(x,a,b)),x,a,b);
check_betaincinv_inputs(x,a,b);
% Process scalar inputs.
SCALARX = coder.const(eml_is_const(isscalar(x)) && isscalar(x));
SCALARA = coder.const(eml_is_const(isscalar(a)) && isscalar(a));
SCALARB = coder.const(eml_is_const(isscalar(b)) && isscalar(b));
SCALARAB = coder.const(SCALARA && SCALARB);
if SCALARX
    x1 = double(x(1));
end
if SCALARA
    a1 = double(a(1));
    glna1 = gammaln(a1);
end
if SCALARB
    b1 = double(b(1));
    glnb1 = gammaln(b1);
end
if SCALARAB
    lbetascalar = glna1 + glnb1 - gammaln(a1 + b1);
end
% Determine whether to use PARFOR.
n = coder.internal.indexInt(numel(y));
if coder.internal.useParforConst(mfilename,n)
    parfor k = 1:n
        if SCALARX
            xk = x1;
        else
            xk = double(x(k));
        end
        if SCALARA
            ak = a1;
            glna = glna1;
        else
            ak = double(a(k));
            glna = gammaln(ak);
        end
        if SCALARB
            bk = b1;
            glnb = glnb1;
        else
            bk = double(b(k));
            glnb = gammaln(bk);
        end
        if SCALARAB
            lbeta = lbetascalar;
        else
            lbeta = glna + glnb - gammaln(ak + bk);
        end
        y(k) = eml_betaincinv(xk,ak,bk,lbeta,upper);
    end
else
    for k = 1:n
        if SCALARX
            xk = x1;
        else
            xk = double(x(k));
        end
        if SCALARA
            ak = a1;
            glna = glna1;
        else
            ak = double(a(k));
            glna = gammaln(ak);
        end
        if SCALARB
            bk = b1;
            glnb = glnb1;
        else
            bk = double(b(k));
            glnb = gammaln(bk);
        end
        if SCALARAB
            lbeta = lbetascalar;
        else
            lbeta = glna + glnb - gammaln(ak + bk);
        end
        y(k) = eml_betaincinv(xk,ak,bk,lbeta,upper);
    end
end

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

function rval = eml_betaincinv(p,a,b,lbeta,upper)
% Compute the inverse incomplete beta function, lower or upper tail.
% Uses Newton's method to find a root of p = betainc(x,a,b).  Starting
% values are based on Griffiths & Hill, Algorithm AS 109, Appl. Statist.
% (1977) Vol.26, and incorporates the suggested changes in AS R83
% (vol.40(1), p.236, 1991).
coder.internal.prefer_const(upper);
if upper
    tail = 1;
else
    tail = 0;
end
pIn = p;
ftol = 100*eps;
ztol = eps;
% Check for trivial and limit cases.
if ~(a >= 0) % including a == NaN
    rval = coder.internal.nan;
    return
elseif ~(b >= 0) % including b == NaN
    rval = coder.internal.nan;
    return
elseif ~((0 < p) && (p < 1)) % including p == NaN
    if p == 0
        rval = tail;
        return
    elseif p == 1
        rval = 1 - tail;
        return
    else
        rval = coder.internal.nan;
        return
    end
elseif a == 0
    if b > 0
        if (p < 1)
            rval = tail;
        else
            rval = 1 - tail;
        end
        return
    else % b == 0)
        if p == 0
            rval = tail;
        elseif p == 1
            rval = 1 - tail;
        else
            rval = coder.internal.nan;
        end
        return
    end
elseif b == 0 % && a > 0
    if p == 0
        rval = tail;
    else
        rval = 1 - tail;
    end
    return
elseif a == coder.internal.inf
    if b < coder.internal.inf
        if p == 0
            rval = tail;
        else
            rval = 1 - tail;
        end
        return
    else % b == inf
        if p == 0
            rval = tail;
        elseif p == 1
            rval = 1 - tail;
        else
            rval = coder.internal.nan;
        end
        return
    end
elseif b == coder.internal.inf % && (a < inf)
    if p < 1
        rval = tail;
    else
        rval = 1 - tail;
    end
    return
end
% Work with the smaller of the two tails, ...
if p > 0.5
    p = 1 - p; % no round-off here, p is > 0.5
    % tail = 1 - tail;
    upper = ~upper;
end
% Use the lower tail for starting value calculation, keeping away from 1
if upper
    pLower = 1 - p;
    if pLower == 1
        pLower = 1 - ftol;
    end
else
    pLower = p;
end
% Starting approximation.
r = sqrt(-2*log(pLower));
y = r - (2.30753 + 0.27061*r)/(1 + (0.99229 + 0.04481*r)*r);
if (a > 1) && (b > 1)
    r = (y*y - 3)/6;
    s = 1 / (2*a - 1);
    t = 1 / (2*b - 1);
    h = 2 / (s+t);
    w = y*sqrt(h+r)/h - (t-s)*(r + 5/6 - 2/(3*h));
    z = a/(a + b*exp(2*w));
else
    r = 2*b;
    t = 1/(9*b); t = (1 - t + y*sqrt(t));
    t = r*t*t*t;
    if t > 0
        t = (4*a + r - 2)/t;
        if t > 1
            z = 1 - 2/(t+1);
        else
            z = exp((log(pLower*a) + lbeta)/a);
        end
    else
        z = 1 - exp((log((1-pLower)*b) + lbeta)/b);
    end
end
% Keep starting value away from boundaries.
z = eml_min(eml_max(z,0.0001),0.9999);
am1 = a - 1;
bm1 = b - 1;
% Newton's method: new x = x - f/f'
% Halley's rational method: new x = x - 1/(f'/f - (f"/f')/2)
% This function: f  = betainc(x,a,b) - p,
%              f' = x^(a-1)*(1-x)^(b-1)/beta(a,b),
%              f" = ((a-1)/x - (b-1)/(1-x))*f'
oldf = coder.internal.inf;
oldz = coder.internal.inf;
maxIter = 1000;
% Impose a lower limit of realmin on the convergence tolerance for the function value
if p > 1.0021e-294
    ftol1 = ftol*p;
else
    ftol1 = ftol*1.0021e-294;
end
if upper
    sgn = -1;
else
    sgn = 1;
end
zlo = 0;
zhi = 1;
i = 0;
while i < maxIter
    f = sgn*(eml_betainc(z,a,b,lbeta,upper) - p);
    % Halve the step if it overshot the root badly.
    if f*oldf < 0
        if abs(oldf) <= abs(f)
            z = .5*(z+oldz);
            f = sgn*(eml_betainc(z,a,b,lbeta,upper) - p);
        end
    end
    % Update the brackets around the root.
    if (f > 0)
        zhi = z;
    else
        zlo = z;
    end
    % If (betainc(z,a,b) - p) is small enough, exit.  ...
    if abs(f) < ftol1
        break
    end
    if abs(z-oldz) < ztol*z+realmin
        break;
    end
    % Use Newton's method.  Halley's method seems to improve some cases
    % a little, but worsen others a lot.
    if i < 200
        oldz = z;
        oldf = f;
        logz = log(z);
        u = 1-z;
        if u ~= 1 % log(1-z)
            log1mz = log(u)*(-z/(u-1));
        else
            log1mz = -z;
        end
        df = exp(am1*logz + bm1*log1mz - lbeta); % this may suffer from subtraction cancellation
        z = z - f/df;
        % z = z - f/(f - 0.5*f*(am1/z - bm1/(1-z)));
        % Make sure we don't step out of the brackets.
        if z <= zlo
            if zlo == 0
                if abs(double(upper)-p) < abs(eml_betainc(realmin,a,b,lbeta,upper)-p)
                    z = 0;
                    break
                else
                    zlo = realmin;
                end
            end
            z = 0.99*zlo + 0.01*zhi;
        elseif z >= zhi
            if zhi == 1
                if abs(1-double(upper)-p) < abs(eml_betainc(1-0.5*eps,a,b,lbeta,upper)-p)
                    z = 1;
                    break
                else
                    zhi = 1 - 0.5*eps;
                end
            end
            z = 0.01*zlo + 0.99*zhi;
        end
    else
        % Reasons for Newton's method not converging: ...
        z = 0.5*(zlo + zhi);
    end
    i = i + 1;
end
if i >= maxIter
    if ~coder.setwarningflag
        coder.internal.warning('Coder:MATLAB:betaincinv_FailedToConverge', ...
            coder.internal.flt2str(pIn),coder.internal.flt2str(a),coder.internal.flt2str(b));
    end
end
rval = z;

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

function check_betaincinv_inputs(x,a,b)
%MATLAB Code Generation Private Function

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

for k = 1:numel(a)
    if ~(a(k) >= 0) % including a == NaN
        coder.internal.error('Coder:MATLAB:betaincinv_PositiveZ');
    end
end
for k = 1:numel(b)
    if ~(b(k) >= 0) % including b == NaN
        coder.internal.error('Coder:MATLAB:betaincinv_PositiveW');
    end
end
for k = 1:numel(x)
    if ~((0 <= x(k)) && (x(k) <= 1)) % including x == NaN
        coder.internal.error('Coder:MATLAB:betaincinv_YOutOfRange');
    end
end

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