picomath

-module(gamma).
-export([gamma/1, log_gamma/1]).

% Visit http://www.johndcook.com/stand_alone_code.html for the source of this code and more like it.

% Note that the functions Gamma and LogGamma are mutually dependent.

gamma(X) when X > 0 ->

    % Split the function domain into three intervals:
    % (0, 0.001), [0.001, 12), and (12, infinity)

    if

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % First interval: (0, 0.001)
        %
        % For small x, 1/Gamma(x) has power series x + gamma x^2  - ...
        % So in this range, 1/Gamma(x) = x + gamma x^2 with error on the order of x^3.
        % The relative error over this interval is less than 6e-7.

        X < 0.001 ->
            Gamma = 0.577215664901532860606512090, % Euler's gamma constant
            1.0/(X*(1.0 + Gamma*X));

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Second interval: [0.001, 12)

        X < 12.0 ->
            % The algorithm directly approximates gamma over (1,2) and uses
            % reduction identities to reduce other arguments to this interval.
            
            Arg_was_less_than_one = (X < 1.0),

            % Add or subtract integers as necessary to bring y into (1,2)
            % Will correct for this below
            [Y, N] = if
                Arg_was_less_than_one ->
                    [X + 1.0, 0];
                true ->
                    N1 = trunc(X) - 1,
                    [X - N1, N1]
            end,

            % numerator coefficients for approximation over the interval (1,2)
            P =
            [
                -1.71618513886549492533811E+0,
                 2.47656508055759199108314E+1,
                -3.79804256470945635097577E+2,
                 6.29331155312818442661052E+2,
                 8.66966202790413211295064E+2,
                -3.14512729688483675254357E+4,
                -3.61444134186911729807069E+4,
                 6.64561438202405440627855E+4
            ],

            % denominator coefficients for approximation over the interval (1,2)
            Q =
            [
                -3.08402300119738975254353E+1,
                 3.15350626979604161529144E+2,
                -1.01515636749021914166146E+3,
                -3.10777167157231109440444E+3,
                 2.25381184209801510330112E+4,
                 4.75584627752788110767815E+3,
                -1.34659959864969306392456E+5,
                -1.15132259675553483497211E+5
            ],

            Z = Y - 1,
            [Num, Den] = gamma_iter(Z, 0.0, 1.0, P, Q),
            Result = Num/Den + 1.0,

            % Apply correction if argument was not initially in (1,2)
            if
                Arg_was_less_than_one ->
                    % Use identity gamma(z) = gamma(z+1)/z
                    % The variable "result" now holds gamma of the original y + 1
                    % Thus we use y-1 to get back the orginal y.
                    Result / (Y - 1.0);
                true ->
                    % Use the identity gamma(z+n) = z*(z+1)* ... *(z+n-1)*gamma(z)
                    gamma_z_n(Result, Y, N)
            end;

        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        % Third interval: [12, infinity)

        X < 171.624 ->
            math:exp(log_gamma(X));

        true ->
            999.0e299

    end.

gamma_iter(_Z, Num, Den, [], []) ->
    [Num, Den];
gamma_iter(Z, Num, Den, [P|Ps], [Q|Qs]) ->
    NewNum = (Num + P) * Z,
    NewDen = Den * Z + Q,
    gamma_iter(Z, NewNum, NewDen, Ps, Qs).

gamma_z_n(Result, _Y, 0) -> Result;
gamma_z_n(Result, Y, N) ->
    gamma_z_n(Result * Y, Y + 1, N - 1).
    

log_gamma(X) when X > 0 ->

    if
        X < 12.0 ->
            math:log(abs(gamma(X)));

        true ->

            % Abramowitz and Stegun 6.1.41
            % Asymptotic series should be good to at least 11 or 12 figures
            % For error analysis, see Whittiker and Watson
            % A Course in Modern Analysis (1927), page 252

            C = lists:reverse(
            [
                 1.0/12.0,
                -1.0/360.0,
                 1.0/1260.0,
                -1.0/1680.0,
                 1.0/1188.0,
                -691.0/360360.0,
                 1.0/156.0,
                -3617.0/122400.0
            ]),
            Z = 1.0/(X*X),
            Sum = log_gamma_iter(Z, 0, C),
            Series = Sum/X,

            HalfLogTwoPi = 0.91893853320467274178032973640562,
            LogGamma = (X - 0.5)*math:log(X) - X + HalfLogTwoPi + Series,    
            LogGamma
    end.

log_gamma_iter(_Z, Sum, []) ->
    Sum;
log_gamma_iter(Z, Sum, [C|Cs]) ->
    S = (Sum * Z) + C,
    log_gamma_iter(Z, S, Cs).