Fast High Precision Gamma Calculator

Introduction

What is Gamma Function?

In mathematics, the gamma function is an extension of the factorial function, with its argument shifted down by 1, to real and complex numbers (look at Wikipedia).

Simple Explanation: What is the factorial of non-Integer number? The answer is gamma function!

In this tip, I show you how to calculate gamma function. Simple, fast and accurate!

Background

First of all, you need to know some basic rules about factorials:

• Γ(x) = (x - 1)!
• x! = (x -1)! * x
• x! = (x+1)! / x

To calculate gamma function, we use taylor series.

Note: taylor series does not converge fast enough for all real numbers. So, we expand this series at the range of numbers between -1 and 1. stepping size is 0.05 and we calculate 9 terms of series to achieve high precision.

This is the taylor series of x! (look at Wolframalpha).

where n0 is expansion point and Γ^(k) is the k^th derivative of gamma function.

Note: As you choose n0 closer to x, you will get more accurate results. Thus, choose the closest n0 to x.
Note: Γ^(k) (1 + n0)/ k! is a constant so we precalculate all these values and keep them in arrays.

Note: To calculate Γ(x), for x between -1 and 1, use taylor series directly.

Note: To calculate Γ(x), for x more than 1 or less than -1, first use factorial rules to shift x between -1 and 1, then multiply result by Γ(x) using series.

Using the Code

Store all precalculated Gamma derivatives inside arrays.

Initialize arrays once. Because these arrays are constant, so mark them as readonly.

You will have 41 arrays (because -20 < x / 0.05 < 20), each array has a length of 10.

First index is n0 and rest of the arrays are constant part of 9 terms of the series.

Since step size is 0.05, we should choose the array that is closest to the x.

For example, when x = 0.5784, then we have x/0.05=11.568 therefore we choose 11^th or 12^th array.

static class GammaDerivative
{
    private static readonly double[]
        GammaDerivativesRange20N =
        {
            0D, Double.NaN, 0D, 0D, 0D, 0D, 0D, 0D, 0D, 0D
        };
    private static readonly double[]
        GammaDerivativesRange19N = {
            -0.95, 19.4700853112555D, -399.0947906775D,
            7999.22620441159D, -159999.192130127D,
            3199999.2329081D, -63999999.259529D, 1279999999.29272D, -25599999999.3247D,
            511999999999.356D
        };
        
    private static readonly double[]
        GammaDerivativesRange18N = {
            -0.90, 9.51350769866873D, -99.166472874637D,
            999.336741961287D, -9999.32866773128D,
            99999.3941978608D, -999999.440326772D, 9999999.48975095D, -99999999.5347559D,
            999999999.576533D
        };
        
    private static readonly double[]
        GammaDerivativesRange17N = {
            -0.85, 6.22027287404988D, -43.6724944505672D,
            295.725305699639D, -1974.74565909015D,
            13168.2411259365D, -87791.0668614162D, 585276.261230467D, -3901843.90519612D,
            26012294.5899988D
        };
        
    private static readonly double[]
        GammaDerivativesRange16N = {
            -0.80, 4.5908437119988D, -24.281155551781D,
            124.506711961476D, -624.523806406947D,
            3124.61123747756D, -15624.6683906425D, 78124.723097606D, -390624.768282098D,
            1953124.80661177D
        };
        
    private static readonly double[]
        GammaDerivativesRange15N = {
            -0.75, 3.62560990822191D, -15.3270974171567D,
            63.5726995550479D, -255.593913429678D,
            1023.68469107349D, -4095.74053268546D, 16383.7921939606D, -65535.8329260497D,
            262143.866126271D
        };
        
    private static readonly double[]
        GammaDerivativesRange14N = {
            -0.70, 2.99156898768759D, -10.4780428417585D, 36.6662258367875D, 123.107763832927D,
            411.265075256399D, -1371.53707456811D, 4584.60422680214D, -653.188004194767D,
            1905.18212638032D
        };
        
    private static readonly double[]
        GammaDerivativesRange13N = {
            -0.65, 2.54614697721229D, -7.56478311429081D,
            23.0015368028935D, -66.3366289787721D,
            190.18514185276D, -543.82748921185D, 1574.64401643923D, -222.01988683161D, 555.081610140393D
        };
        
    private static readonly double[]
        GammaDerivativesRange12N = {
            -0.60, 2.21815954375769D, -5.68155957280492D, 15.3452802672229D, -38.7987670035651D,
            97.4812519180141D, -244.009046282463D, 610.258055159647D, -1525.81152787317D,
            3814.64906436589D
        };
        
    private static readonly double[]
        GammaDerivativesRange11N = {
            -0.55, 1.96813640060238D, -4.39590870900005D,
            10.7313082744161D, -24.1547234259443D,
            54.046946033126D, -120.320539267501D, 267.543235620895D, -594.651805969081D,
            1321.52636874769D
        };
        
    private static readonly double[]
        GammaDerivativesRange10N = {
            -0.50, 1.77245385090552D, -3.48023090691326D,
            7.79008872120313D, -15.7947670515358D,
            31.8788248211608D, -63.9126957469214D, 127.942612101477D, -255.961228074064D,
            511.974131361039D
        };
        
    private static readonly double[]
        GammaDerivativesRange9N = {
            -0.45, 1.61612426873358D, -2.80552587449537D,
            5.82967843679375D, -10.7451701862446D,
            19.7681401306886D, -36.054376432714D, 65.6387780703025D, -119.396066949852D,
            217.118739506589D
        };
        
    private static readonly double[]
        GammaDerivativesRange8N = {
            -0.40, 1.48919224881282D, -2.29427819170183D,
            4.47481216055898D, -7.55158129739921D,
            12.7751466003399D, -21.3737135387477D, 35.6862016537667D, -59.5142819897528D,
            99.2145860007862D
        };
        
    private static readonly double[]
        GammaDerivativesRange7N = {
            -0.35, 1.38479510202651D, -1.89765322148743D,
            3.50997923688594D, -5.45315568098083D,
            8.54726981248489D, -13.2092750442655D, 20.3698140233823D, -31.3648094688375D,
            48.2704560546125D
        };
        
    private static readonly double[]
        GammaDerivativesRange6N = {
            -0.30, 1.29805533264756D, -1.58365807983323D,
            2.80542638961667D, -4.02911317569611D,
            5.89012621801116D, -8.4577012051428D, 12.1192152879838D, -17.3323878441919D,
            24.772584430011D
        };
        
    private static readonly double[]
        GammaDerivativesRange5N = {
            -0.25, 1.22541670246518D, -1.33063205864388D,
            2.27987153689211D, -3.03563250547648D,
            4.16393568279721D, -5.58283545990313D, 7.47255601476381D, -9.97739134364727D,
            13.3118755327171D
        };
        
    private static readonly double[]
        GammaDerivativesRange4N = {
            -0.20, 1.1642297137253D, -1.12349164735875D,
            1.88064759057317D, -2.32570573138395D,
            3.00999174853583D, -3.78401986966253D, 4.75293897617795D, -5.95139905960929D,
            7.44560973380385D
        };
        
    private static readonly double[]
        GammaDerivativesRange3N = {
            -0.15, 1.11248373694847D, -0.951474631732596D,
            1.57262112883441D, -1.8076083166106D,
            2.21911024240823D, -2.62497750703055D, 3.10679421302081D, -3.66254589040758D,
            4.313599542982D
        };
        
    private static readonly double[]
        GammaDerivativesRange2N = {
            -0.10, 1.06862870211932D, -0.806736606716814D,
            1.33175429379401D, -1.42237555998628D,
            1.66505223223746D, -1.85860162426187D, 2.08047822528401D, -2.31712549060625D,
            2.57813819793116D
        };
        
    private static readonly double[]
        GammaDerivativesRange1N = {
            -0.05, 1.03145331712903D, -0.683450886596578D, 
            1.14123138250983D, -1.1310941970201D,
            1.26930701966673D, -1.34008996998964D, 1.42357900301788D, -1.50248828254655D,
            1.58428059433746D
        };
        
    private static readonly double[]
        GammaDerivativesRange0 = {
            0.00, 1D, -0.577215664901533D, 0.989055995327973D,
            -0.907479076080886D, 0.9817280868344D,
            -0.981995068903145D, 0.993149114621276D, -0.996001760442431D, 0.998105693783129D
        };
        
    private static readonly double[]
        GammaDerivativesRange1 = {
            0.05, 0.973504265562776D, -0.484654222619487D,
            0.866519543079475D, -0.733402094741922D,
            0.76951527826578D, -0.730068352575782D, 0.705106859611967D, -0.673517910757183D,
            0.643103417607824D
        };
        
    private static readonly double[]
        GammaDerivativesRange2 = {
            0.10, 0.951350769866873D, -0.403139588794969D,
            0.767201321491692D, -0.596124811841251D,
            0.610752054797181D, -0.54983481640542D, 0.508648322752892D, -0.4637246400447D,
            0.422897453990549D
        };
        
    private static readonly double[]
        GammaDerivativesRange3 = {
            0.15, 0.933040931107482D, -0.330601293535208D,
            0.68630140437857D, -0.486543163883966D,
            0.490509800318163D, -0.418903275959234D, 0.372323153930381D, -0.324548950953126D,
            0.28330369492238D
        };
        
    private static readonly double[]
        GammaDerivativesRange4 = {
            0.20, 0.918168742399761D, -0.265387398357406D,
            0.620186840514084D, -0.39804931991382D,
            0.398441088563801D, -0.322440650940798D, 0.276228878711926D, -0.230558813523784D,
            0.19304025682217D
        };
        
    private static readonly double[]
        GammaDerivativesRange5 = {
            0.25, 0.906402477055477D, -0.206164446067268D,
            0.566077471605273D, -0.325778802371648D,
            0.327259338693816D, -0.250442097877781D, 0.20751580468742D, -0.166037551807169D,
            0.133605517985885D
        };
        
    private static readonly double[]
        GammaDerivativesRange6 = {
            0.30, 0.897470696306277D, -0.151843864839965D,
            0.521824909277723D, -0.266103313090753D,
            0.271758743992428D, -0.196047114034794D, 0.157734240896269D, -0.121073949387992D,
            0.0938120535339359D
        };
        
    private static readonly double[]
        GammaDerivativesRange7 = {
            0.35, 0.891151442024301D, -0.101527112789497D,
            0.485754766721896D, -0.216283339676781D,
            0.228170669693991D, -0.154479371387366D, 0.121233525593135D, -0.0893005718326651D,
            0.0667563410353867D
        };
        
    private static readonly double[]
        GammaDerivativesRange8 = {
            0.40, 0.887263817503075D, -0.0544642853642779D,
            0.456552534084234D, -0.174226534203179D,
            0.193733763640491D, -0.122366594971248D, 0.0941757813954528D, -0.0665559896195936D,
            0.0480978731915532D
        };
        
    private static readonly double[]
        GammaDerivativesRange9 = {
            0.45, 0.885661380271072D, -0.0100225184476424D,
            0.433180014487195D, -0.138317267258814D,
            0.166402288962429D, -0.0972966372492549D, 0.0739167619019626D, -0.0500770651916732D,
            0.0350599673777081D
        };
        
    private static readonly double[]
        GammaDerivativesRange10 = {
            0.50, 0.886226925452758D, 0.032338397448885D, 0.414813453688301D, -0.107294804564772D,
            0.144645359044622D, -0.0775230522998542D, 0.0586103038171763D, -0.0380019355548651D,
            0.0258376064557562D
        };
        
    private static readonly double[]
        GammaDerivativesRange11 = {
            0.55, 0.888868347803466D, 0.0730850377611226D,
            0.400797265741192D, -0.0801651656407642D,
            0.127306885634156D, -0.0617669073040915D, 0.0469515059524152D, -0.0290587521158863D,
            0.0192397787722501D
        };
        
    private static readonly double[]
        GammaDerivativesRange12 = {
            0.60, 0.89351534928769D, 0.112625333791716D,
            0.390609104633555D, -0.0561366178805398D,
            0.113506662804717D, -0.0490815229087629D, 0.0380074535123248D, -0.0223675400849471D,
            0.0144696107189241D
        };
        
    private static readonly double[]
        GammaDerivativesRange13 = {
            0.65, 0.900116816317232D, 0.151320508059681D,
            0.383833282488433D, -0.0345719557515931D,
            0.102569697134353D, -0.0387589662876536D, 0.0311040709330098D, -0.0173121313620963D,
            0.0109869666606985D
        };
        
    private static readonly double[]
        GammaDerivativesRange14 = {
            0.70, 0.90863873285329D, 0.189494676764298D,
            0.380140392898439D, -0.0149528333706117D,
            0.0939751769116959D, -0.0302646255888054D, 0.025749496445827D, -0.0134562029505838D,
            0.00842125681577184D
        };
        
    private static readonly double[]
        GammaDerivativesRange15 = {
            0.75, 0.919062526848883D, 0.227442658482271D,
            0.379271594025204D, 0.00314715778474762D,
            0.0873192566214339D, -0.0231909121301314D, 0.0215815511697291D, -0.010487492971643D,
            0.00651530589057573D
        };
        
    private static readonly double[]
        GammaDerivativesRange16 = {
            0.80, 0.931383770980243D, 0.265436395838301D,
            0.38102642509978D, 0.0200830054660084D,
            0.0822876674447171D, -0.0172241471941937D, 0.0183313112798304D, -0.0081802715094781D,
            0.00508872743379526D
        };
        
    private static readonly double[]
        GammaDerivativesRange17 = {
            0.85, 0.945611176406195D, 0.303730299975759D,
            0.385253327776652D, 0.0361540597154019D,
            0.0786353894364008D, -0.012120638567738D, 0.0157975740371346D, -0.00636979382563056D,
            0.00401372112712226D
        };
        
    private static readonly double[]
        GammaDerivativesRange18 = {
            0.90, 0.961765831907387D, 0.342565756074187D,
            0.391842257697791D, 0.0516162898063564D,
            0.0761714490274338D, -0.00768922959822524D, 0.0138287784937392D, -0.00493471626161532D,
            0.00319888753179316D
        };
        
    private static readonly double[]
        GammaDerivativesRange19 = {
            0.95, 0.979880651272581D, 0.382174974862283D,
            0.400718926787762D, 0.0666918953407411D,
            0.0747474716632996D, -0.00377845182343007D, 0.0123100828773469D, -0.00378486540133988D,
            0.0025782820740309D
        };
        
        private static readonly double[]
        GammaDerivativesRange20 =
        {
            0D, 1D, 0D, 0D, 0D, 0D, 0D, 0D, 0D, 0D
        };
        
    public static double[] GetGammaArray(double x)
    {
            switch ((int) (x / 0.05))
            {
                case -20:
                    return GammaDerivativesRange20N; // n0 = -1 :only happens 
                    		//when input (x) is non-positive Integer.
                case -19:
                    return GammaDerivativesRange19N; // n0 = -0.95 & -1.00 < x <= -0.95
                case -18:
                    return GammaDerivativesRange18N; // n0 = -0.90 & -0.95 < x <= -0.90
                case -17:
                    return GammaDerivativesRange17N; // n0 = -0.85 & -0.90 < x <= -0.85
                case -16:
                    return GammaDerivativesRange16N; // n0 =  .... & (n0 - 0.05) < x <= n0
                case -15:
                    return GammaDerivativesRange15N;
                case -14:
                    return GammaDerivativesRange14N;
                case -13:
                    return GammaDerivativesRange13N;
                case -12:
                    return GammaDerivativesRange12N;
                case -11:
                    return GammaDerivativesRange11N;
                case -10:
                    return GammaDerivativesRange10N;
                case -9:
                    return GammaDerivativesRange9N;
                case -8:
                    return GammaDerivativesRange8N;
                case -7:
                    return GammaDerivativesRange7N;
                case -6:
                    return GammaDerivativesRange6N;
                case -5:
                    return GammaDerivativesRange5N;
                case -4:
                    return GammaDerivativesRange4N;
                case -3:
                    return GammaDerivativesRange3N;
                case -2:
                    return GammaDerivativesRange2N;
                case -1:
                    return GammaDerivativesRange1N;
                case 0:
                    return GammaDerivativesRange0; // n0 = 0 & -0.05 < x < 0.05
                case 1:
                    return GammaDerivativesRange1; // n0 = 0.05 & 0.05 <= x < 0.10
                case 2:
                    return GammaDerivativesRange2; // n0 = 0.10 & 0.10 <= x < 0.15
                case 3:
                    return GammaDerivativesRange3; // n0 = 0.15 & 0.15 <= x < 0.20
                case 4:
                    return GammaDerivativesRange4; // n0 = .... & n0 <= x < (n0 + 0.5)
                case 5:
                    return GammaDerivativesRange5;
                case 6:
                    return GammaDerivativesRange6;
                case 7:
                    return GammaDerivativesRange7;
                case 8:
                    return GammaDerivativesRange8;
                case 9:
                    return GammaDerivativesRange9;
                case 10:
                    return GammaDerivativesRange10;
                case 11:
                    return GammaDerivativesRange11;
                case 12:
                    return GammaDerivativesRange12;
                case 13:
                    return GammaDerivativesRange13;
                case 14:
                    return GammaDerivativesRange14;
                case 15:
                    return GammaDerivativesRange15;
                case 16:
                    return GammaDerivativesRange16;
                case 17:
                    return GammaDerivativesRange17;
                case 18:
                    return GammaDerivativesRange18;
                case 19:
                    return GammaDerivativesRange19;
                case 20:
                    return GammaDerivativesRange20; // n0 = 1 :only happens 
                    		//when input (x) is positive Integer.
            }
            return null;
    }
}

Note: double is floating point so sometimes method will not correctly select array (example: it will choose next or previous array), but it will not affect the result more than Epsilon!

Note: If range equals to -20 or 20 means x was Integer. If you notice, we have two more arrays: GammaDerivativesRange20 and GammaDerivativesRange20N. These arrays are filled with 0 and contain exception.

GammaDerivativesRange20N[1] = Double.NaN //because factorial of non-positive Integer is not defined

GammaDerivativesRange20 [1] = 1 //the final expression would be resault*1.

The main Factorial and Gamma method:

public static double Gamma(double x)
{
    return (x - 1d).Factorial();//according to factorial rules.
}

public static double Factorial(this double x)
{
    if (x.Equals(Double.NaN)) return Double.NaN;//in case that x is NaN.
    if (x >= 170.602544778489474d)
    return Double.PositiveInfinity; //values more than this will end up to positive infinity.
    if (x <= -178.477785816716264d) return System.Math.Abs((x - (int) x))
    < Double.Epsilon? Double.NaN : 0; //values less than this will be NaN(for integer) 
      or 0(for non-Integer)
    var resault = 1d; // this will hold factorial of numbers
    if (x > 0) while (x > 1) resault *= x--;// normal factorial
    //algorithm except that resault is type of double and counter is x itself
    else while (x < -1) resault /= ++x; // factorial algorithm for negative values.
    var gammaDerivativesRange = GammaDerivative.GetGammaArray(x);//select array.
    // Gamma(k)(1+n0)/k! . constant part of series.
    x -= gammaDerivativesRange[0]; // (x-n0)^k , n0 is first index of array. look at taylor series.
    return
        ((((((((gammaDerivativesRange[9]*x + gammaDerivativesRange[8])*x + gammaDerivativesRange[7])*x +
              gammaDerivativesRange[6])*x + gammaDerivativesRange[5])*x + gammaDerivativesRange[4])*x +
           gammaDerivativesRange[3])*x + gammaDerivativesRange[2])*x +
           gammaDerivativesRange[1])*resault; // series expression using Horner's rule
}

Run the code:

static void Main()
{
    Console.WriteLine(Gamma(Double.Parse(Console.ReadLine())));
}

Compare results at www.wolframalpha.com !

Accuracy

With taylor series expanded at 41 different points up to 9 terms, this algorithm will give you 15 decimal parts all correct!

Performance

For optimization i - introduced arrays inside class and initialized them once. Performance difference: about 0.08 secs.

• Used expression when returning result except for using for loop. Performance difference is about 3.5 secs !!
• Used parenthesis inside expression: Performance difference: about 1 sec !
• Removed separate factorial algorithm for integer numbers. Performance difference: about 0.5 secs
• Used Horners rule for evaluating expression: Performance difference: about 0.05 secs

(Thanks to YvesDaoust for his suggestion.)

Take a look at this benchmark.

static void Main()
{
    for (int i = 0; i < 10; i++)
    {
        var timer = new Stopwatch();
        timer.Start();
        for (var j = -10d; j < 10d; j += 0.000001)//20,000,000 iteration.
        {
            var r = Math.Factorial(j);
        }
        Console.WriteLine(timer.Elapsed);
    }
}

Results:

00:00:00.5542497
00:00:00.5417712
00:00:00.5415342
00:00:00.5416845
00:00:00.5418725
00:00:00.5413725
00:00:00.5441604
00:00:00.5412434
00:00:00.5402805
00:00:00.5426584

History

• Returning final result expression changed based on horners scheme (thanks to YvesDaoust!)
• Added exception for double.NaN inputs
• A few changes to switch case, removed if else condition