Title: Calculate the gamma function in C#

The gamma function, represented by the capital Greek letter gamma Γ, was derived by Daniel Bernoulli for complex numbers with a positive real part. The interesting thing about the gamma function is that Γ(n + 1) = n! for integers n greater than zero. The function is continuous for positive integers, so it defines interpolated factorial values for non-integers. For example, 3! = Γ(3 + 1) = 6 and 4! = Γ(4 + 1) = 12. If you want to know what 3.5! would be, you can calculate Γ(3.5 + 1) ≈ 11.631728396567.

Bernoulli defined the gamma function using the following improper integral.

Calculating the Gamma Function

Bernoulli's integral doesn't have a closed-form solution, but the following function approximates its value.

// Integrate: x^(z-1)*e^(-x) dx from 0 to infinity. private double Gamma(double z, double dx, int num_slices) { double result = 0; double x = 0; for (int i = 0; i < num_slices; i++) { double new_term = Math.Pow(x, z - 1) * Math.Exp(-x); if (double.IsNaN(new_term)) break; if (double.IsInfinity(new_term)) break; result += new_term; x += dx; } return result * dx; }

This function basically adds up num_slices slices of curve that are dx thick. (This is covered in my book Essential Algorithms, Second Edition.) To get a good result, use a small value for dx and a large value for num_slices. If you use too many slices, however, rounding errors will add up so you may not improve the resulting accuracy.

The method first initializes result to zero. It then loops through the slices, calculating the height of the function for each slice. When x is big enough, a slice's height may turn out to be not a number (NaN) or infinity, at least as far as the double data type is concerned. In that case, the code breaks out of its loop and stops considering new slices because further values will not be helpful.

After is has added up all of the slices, the code multiplies the total height by dx to convert the heights into areas and returns the result.

There are other ways to approximate the gamma function. Believe it or not, this one seemed simplest. See Wikipedia for details.

Displaying Gamma Function Values

The example program uses the following code to display factorial values and their corresponding gamma function values.

private void btnCalculate_Click(object sender, EventArgs e) { lvwValues.Items.Clear(); Cursor = Cursors.WaitCursor; int minx = int.Parse(txtMinX.Text); int maxx = int.Parse(txtMaxX.Text); double dt = double.Parse(txtDt.Text); int num_slices = int.Parse(txtNumSlices.Text); for (int x = minx; x <= maxx; x++) { double x_factorial = Factorial(x); if (double.IsInfinity(x_factorial)) break; double gamma_x = Gamma(x + 1, dt, num_slices); double difference = Math.Abs(gamma_x - x_factorial); double percent_difference = difference / gamma_x * 100.0; lvwValues.AddRow( x.ToString("G4"), x_factorial.ToString("G4"), gamma_x.ToString("G4"), difference.ToString("G4"), percent_difference.ToString("G4")); } lvwValues.AutoSizeColumns(); Cursor = Cursors.Default; }

This code gets the parameters that you entered on the form and then loops through the desired X values. For each value, the code first calculates the factorial. If that value is too big to fit in a double, the code breaks out of its loop. (Because it won't be able to calculate larger factorial values, either.)

Next the code calls the gamma function for the value plus one, calculates the difference between the factorial and the gamma function, and calculates the percentage difference. It then uses the AddRow extension method (described shortly) to add the values to the program's ListView control. You can see from the picture at the top of the post that the error is less than one with the given parameters for values up to 15! = Γ(15 + 1).

ListView Extensions

This example uses two useful ListView control extensions. The following AddRow extension adds an item and subitems to create a row in the ListView control.

// Add an item and subitems to a ListView. public static void AddRow(this ListView lvw, string item, params string[] subitems) { ListViewItem list_view_item = lvw.Items.Add(item); foreach (string subitem in subitems) list_view_item.SubItems.Add(subitem); }

The AddRow extension uses its first parameter to create a new ListView item. It then loops through its other parameters and adds them as subitems.

The following code automatically sizes a ListView control's column widths.

// Autosize all columns. public static void AutoSizeColumns(this ListView lvw, params string[] values) { foreach (ColumnHeader column in lvw.Columns) column.Width = -2; }

This code simply loops through the ListView columns and sets their widths to -2. That tells each column to set its width so it is big enough to hold its values.

Conclusion

The gamma function gives you a new way to calculate factorials, but the obvious way of multiplying integers together is much faster. The real reason why the gamma function is interesting is that it lets you calculate factorials for non-integer values.

The program can calculate values up to 170!. Values larger than that are too big for the double data type.

If you look at the picture at the top of the post, you'll see that the differences between the factorial and gamma function are quite small. As the values grow larger, so do the errors. The following picture shows the results between 125! and 141!.

Notice that the first errors are large but still a small percentage of the values. For example, the error in calculating 125! is more than 1x10¹⁹⁶, but it's less than 0.0000000000068% of the value 125!.

For larger values, the percent error quickly grows much larger. At those larger values, rounding errors dominate so you can't reduce the errors too much if you use the double data type.

If you like, try implementing some of the other ways for calculating the gamma function and see which give the best results.

Download the example to experiment with it and to see additional details.