function x = gevinv(p,k,sigma,mu)
%GEVINV Inverse of the generalized extreme value cumulative distribution function, optimised for GPU.
%   X = GEVINV(P,K,SIGMA,MU)
%
%   See also GEVINV, GPUARRAY.

%   See toolbox/stats/stats/gevinv for implementation details.

%   Copyright 2021 The MathWorks, Inc.

if nargin < 1
    error(message('stats:gevinv:TooFewInputs'));
end
if nargin < 2, k = 0;     end
if nargin < 3, sigma = 1; end
if nargin < 4, mu = 0;    end

err = internal.stats.statsizechk(4,p,k,sigma,mu);
if err > 0
    error(message('stats:gevinv:InputSizeMismatch'));
end

% Return NaN for out of range parameters.
sigma(sigma <= 0) = NaN;

% Work out the expected output type
outProto = cast(1,"like",internal.stats.dominantType(p,k,sigma,mu));

% Process each element independently using arrayfun. This causes the
% function to be compiled into an optimised GPU kernel.
x = arrayfun(@gevinv_scalar, p, k, sigma, mu, outProto);
end

function x = gevinv_scalar(p, k, sigma, mu, outProto)
% Process one element of the GEVPDF output
if (0<p) && (p<1)
    % P ok
    if abs(k) < eps
        % Find the k==0 cases and fill them in.
        z = cast(-log(-log(p)),"like",outProto);
    else
        % Find the k~=0 cases and fill them in.
        z = cast(expm1(-k*log(-log(p)))./k,"like",outProto); % ((-log(u)).^(-k) - 1) ./ k
    end
else
    % When k~=0, the support is 1 + k.*(x-mu)/sigma > 0
    if p==0 && k<=0
        z = -inf("like",outProto);
    elseif p==1 && k<0
        z = cast(-1/k,"like",outProto);
    elseif p==0 && k>0
        z = cast(-1/k,"like",outProto);
    elseif p==1 && k>=0
        z = inf("like",outProto);
    else
        z = nan("like",outProto);
    end
end

x = mu + sigma*z;
end