function delta = stirlerr_scalar(n)
%STIRLERR Error of Stirling-De Moivre approximation to n factorial,
%optimised for GPU. All inputs must be scalar.

%   Copyright 2019 The MathWorks, Inc.

delta = 0*n; % To get the right type
nn = n*n;

% Define S0=1/12 S1=1/360 S2=1/1260 S3=1/1680 S4=1/1188
S0 = 8.333333333333333e-02;
S1 = 2.777777777777778e-03;
S2 = 7.936507936507937e-04;
S3 = 5.952380952380952e-04;
S4 = 8.417508417508418e-04;

if n<=15
    n1 = n;
    n2 = 2*n1;
    if (n2==round(n2))
        delta = delta + sfe(n2+1); % Use PLUS to preserve type
    else
        lnsr2pi = 0.9189385332046728;
        delta = gammaln(n1+1)-(n1+0.5)*log(n1)+n1-lnsr2pi;
    end
elseif (n>15 && n<=35)
    delta = (S0-(S1-(S2-(S3-S4/nn)/nn)/nn)/nn)/n;
elseif (n>35 && n<=80)
    delta = (S0-(S1-(S2-S3/nn)/nn)/nn)/n;
elseif (n>80 && n<=500)
    delta = (S0-(S1-S2/nn)/nn)/n;
elseif (n>500)
    delta = (S0-S1/nn)/n;
end

end

function out = sfe(in)
% We can't use indexing or switchyards in GPU arrayfun, so a long list of
% "if" is the best we can do.

out = 0;
if in==1
    out = 0;
elseif in==2
    out = 1.534264097200273e-01;
elseif in==3
    out = 8.106146679532726e-02;
elseif in==4
    out = 5.481412105191765e-02;
elseif in==5
    out = 4.134069595540929e-02;
elseif in==6
    out = 3.316287351993629e-02;
elseif in==7
    out = 2.767792568499834e-02;
elseif in==8
    out = 2.374616365629750e-02;
elseif in==9
    out = 2.079067210376509e-02;
elseif in==10
    out = 1.848845053267319e-02;
elseif in==11
    out = 1.664469118982119e-02;
elseif in==12
    out = 1.513497322191738e-02;
elseif in==13
    out = 1.387612882307075e-02;
elseif in==14
    out = 1.281046524292023e-02;
elseif in==15
    out = 1.189670994589177e-02;
elseif in==16
    out = 1.110455975820868e-02;
elseif in==17
    out = 1.041126526197210e-02;
elseif in==18
    out = 9.799416126158803e-03;
elseif in==19
    out = 9.255462182712733e-03;
elseif in==20
    out = 8.768700134139385e-03;
elseif in==21
    out = 8.330563433362871e-03;
elseif in==22
    out = 7.934114564314021e-03;
elseif in==23
    out = 7.573675487951841e-03;
elseif in==24
    out = 7.244554301320383e-03;
elseif in==25
    out = 6.942840107209530e-03;
elseif in==26
    out = 6.665247032707682e-03;
elseif in==27
    out = 6.408994188004207e-03;
elseif in==28
    out = 6.171712263039458e-03;
elseif in==29
    out = 5.951370112758848e-03;
elseif in==30
    out = 5.746216513010116e-03;
elseif in==31
    out = 5.554733551962801e-03;
end
end
