function y = fpdf_scalar(x, v1, v2, outProto)
%FPDF   F probability density function for a single element, optimized for the GPU.

%   Algorithm must be kept in-sync with toolbox/stats/stats/fpdf.m

%   Copyright 2019 The MathWorks, Inc.

% Initialize Y to zero
y = zeros("like", outProto);

% Special case for invalid parameters
if (v1 <= 0) || (v2 <= 0) || isnan(x) || isnan(v1) || isnan(v2)
    y = NaN("like", outProto);
    return
end

const = 1e8;
if (x > 0 && x < Inf) && (v1 > 0 && v1 < Inf) && (v2 > 0 && v2 < Inf)
    % Regular case with valid parameters and positive x
    ratio = v1 .* x ./ v2;
    logtemp = -log(x) - 0.5 * v1 .* log1p(1./ratio) - 0.5 * v2 .* log1p(ratio) - betaln(v1/2, v2/2);
    y = exp(logtemp);
    
elseif ( (v1>const) && (v1<Inf) ) || ( (v2>const) && (v2<Inf) ) && (x > 0 && x < Inf)
    % Special case - if either v1 or v2 is too large then return NaN
    y = NaN("like", outProto);
    
elseif (x==Inf) && (v1>0) && (v2>0)
    % Special case - at x==Inf, result should be 0, regardless of any valid v1 or v2
    % (y is already zero)
    
elseif (x==0) && (v1>0) && (v2>0)
    % Special case - at x==0, result depends only on the value of v1
    if (v1==2)
        y = ones("like", outProto);
    elseif (v1<2)
        y = Inf("like", outProto);
    elseif (v1>2)
        % Leave at zero
    end
    
elseif (v2 > 0 && v2 < Inf) && (v1 == Inf) && (x > 0 && x < Inf)
    % If v1 = Inf AND v2 is not Inf AND both v1, v2 > 0 AND x is not NaN
    % then fpdf( t, v1, v2 ) = (v2/t^2) * chi2pdf( v2/t, v2 )
    y = ( v2 ./ (x.^2) ) .* chi2pdf_scalar( v2 ./ x , v2, outProto );
    
elseif (v1 > 0 && v1 < Inf) && (v2 == Inf) && (x > 0 && x < Inf)
    % If v1 is not Inf AND v2 = Inf AND both v1, v2 > 0 AND x is not NaN
    % then fpdf( t, v1, v2 ) = v1 * chi2pdf( v1*t, v1 )
    y = v1 .* chi2pdf_scalar( v1 .* x, v1, outProto );
    
elseif (v1 == Inf) && (v2 == Inf) && (x > 0 && x < Inf)
    % If v1 is Inf AND v2 = Inf AND x is not NaN
    % then set fpdf( t, v1, v2 ) = 0 if t not equal to 1 and Inf if t = 1
    if (x==1)
        y = Inf("like", outProto);
    else
        y = zeros("like", outProto);
    end
    
end

end

function y = chi2pdf_scalar(x, v, outProto)
y = gampdf_scalar(x, v/2, 2, outProto);
end
