picomath

generic
    type Float_Type is digits <>;
package Gamma is

    function Gamma(x: Float_Type'Base) return Float_Type'Base;
    function Log_Gamma(x: Float_Type'Base) return Float_Type'Base;

end Gamma;

with Ada.Numerics.Generic_Elementary_Functions;

package body Gamma is

    package Elementary_Functions is new Ada.Numerics.Generic_Elementary_Functions(Float_Type'Base);
    use Elementary_Functions;

    function Positive_Infinity return Float_Type'Base is

        result: Float_Type'Base := Float_Type'Base'Last;

    begin
        return Float_Type'Base'Succ(result);
    exception
        when Constraint_Error =>
            return result;
    end Positive_Infinity;

    function Gamma(x: Float_Type'Base) return Float_Type'Base is

        gamma: constant := 0.577215664901532860606512090; -- Euler's gamma constant

        -- numerator coefficients for approximation over the interval (1,2)
        p: constant array(1..8) of Float_Type'Base := (
            -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: constant array(1..8) of Float_Type'Base := (
            -3.08402300119738975254353E+1,
             3.15350626979604161529144E+2,
            -1.01515636749021914166146E+3,
            -3.10777167157231109440444E+3,
             2.25381184209801510330112E+4,
             4.75584627752788110767815E+3,
            -1.34659959864969306392456E+5,
            -1.15132259675553483497211E+5
        );

        y, num, den, z, result: Float_Type'Base;
        n: Natural;
        arg_was_less_than_one: Boolean;

    begin
        if x <= 0.0 then
            raise Ada.Numerics.Argument_Error;
        end if;

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

        ---------------------------------------------------------------------------
        -- 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.

        if x < 0.001 then
            return 1.0/(x*(1.0 + gamma*x));

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

        elsif x < 12.0 then
            -- The algorithm directly approximates gamma over (1,2) and uses
            -- reduction identities to reduce other arguments to this interval.
            
            y := x;
            n := 0;
            arg_was_less_than_one := (y < 1.0);

            -- Add or subtract integers as necessary to bring y into (1,2)
            -- Will correct for this below
            if arg_was_less_than_one then
                y := y + 1.0;
            else
                n := Natural(Float_Type'Base'Floor(y)) - 1;  -- will use n later
                y := y - Float_Type'Base(n);
            end if;

            num := 0.0;
            den := 1.0;

            z := y - 1.0;
            for i in 1..8 loop
                num := (num + p(i))*z;
                den := den*z + q(i);
            end loop;
            result := num/den + 1.0;

            -- Apply correction if argument was not initially in (1,2)
            if arg_was_less_than_one then
                -- 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 := result / (y-1.0);
            else
                -- Use the identity gamma(z+n) = z*(z+1)* ... *(z+n-1)*gamma(z)
                for i in 1..n loop
                    result := result * y;
                    y := y + 1.0;
                end loop;
            end if;

            return result;

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

        elsif x <= 171.624 then
            return Exp(Log_Gamma(x));

        else
            -- Correct answer too large to display. 
            return Positive_Infinity; -- float infinity

        end if;
    end Gamma;

    function Log_Gamma(x: Float_Type'Base) return Float_Type'Base is

        c: constant array(1..8) of Float_Type'Base := (
             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
        );
        halfLogTwoPi: constant := 0.91893853320467274178032973640562;

        z, sum, series: Float_Type'Base;

    begin
        if x <= 0.0 then
            raise Ada.Numerics.Argument_Error;
        end if;

        if x < 12.0 then
            return Log(Abs(Gamma(x)));
        end if;

        -- 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

        z := 1.0/(x*x);
        sum := c(8);
        for i in reverse 1..7 loop
            sum := sum * z;
            sum := sum + c(i);
        end loop;
        series := sum/x;

        return (x - 0.5)*Log(x) - x + halfLogTwoPi + series;
    end Log_Gamma;

end Gamma;