This article guides you through the process of creating a lightweight, zero-dependency Rust crate dedicated to computing decimal factorials using the Gamma Function.

With a focus on simplicity and performance, we’ll explore the inherent power of Rust and the versatility of the Gamma Function to develop a standalone solution for factorial calculations.

Explanation -

Before we delve into the implementation, let’s talk about what a factorial is.

A standard factorial $n!$ is defined as:

For instance, $4!$ expands to $24$ as follows:

A one-line implementation of a standard factorial is as follows:

1
2
3
4
5
6
7
8

let n = 4;

let standard_factorial = (1i64..=n.max(1))
        .fold(1i64, |acc, item| {
            acc * item
        });

assert_eq!(standard_factorial, 24);

This works by creating a Range from 1 to n (or 1, whichever is greater), which implements Iterator conveniently. We then use the built in fold function to multiple all elements to each other, effectively creating a basic factorial function.

You may notice that the input n is of datatype i64. This is intentional - this function does not work for non-integers! In fact, there is no such implementation that does - or is there?

The Gamma Function is a formula that allows the input of any $x \in \mathbb{R}$.

For all values $x \in \mathbb{R} \text{ and } x \ge 0$,

This holds true for all integers - and allows us to calculate decimal values of x! But hold on a minute, that’s a definite integral! Last I checked, Rust doesn’t deal with calculus. Well, thankfully, when integrated, the Gamma Function exhibits a clear pattern.

After using integration by parts, we’re left with this:

The leftmost and middle terms reduce to 0, leaving

or…

For instance, with $\Gamma(4.02 + 1)$ (or $4.02!$):

That first part won’t be too difficult to implement, but what about the second?

We’ll need to approximate that integral. We’ll use a technique called LRAM, or Left Reimann Approximation.

Setup -

Let’s start by creating our project workspace with Cargo!

1

cargo init

Implementation -

We essentially add up rectangular areas of $f(x) * step$. I implemented this as follows:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17

fn definite_integral(f: impl Fn(f64) -> f64, lower: f64, upper: f64, step: f64) -> f64 {
    let mut area: f64 = 0f64;

    let mut x: f64 = lower;
    let mut step: f64 = step;
    while x < upper {
        if (x + step) > upper {
            step = upper - x;
        }

        area += step * f(x);

        x += step;
    }

    area
}

This allows us to estimate that last part, so we can go ahead and start implementing everything.

We’ll start by modifying our standard_factorial function so that it works for floats:

1
2
3
4

let standard_factorial = (1i64..=n.max(1f64) as i64)
        .fold(1f64, |acc, item| {
            acc * (item as f64 + n % 1f64)
        });

This will calculate the $(x)(x - 1)…(x - \lfloor{x}\rfloor + 1)$ portion of our equation. Next, we need to get the factorial of the decimal with our newly implemented integral estimation.

1

let decimal_factorial = definite_integral(|x| x.powf(n - n.floor()) * ((-x).exp()), 0f64, 1000f64, 0.0001f64);

Infinity isn’t exactly an f64, so I’ll use 1000 for our upper bound with a step of 0.0001.

Finally, we need to handle two edge cases: $x < 0$ and $x = \lfloor{x}\rfloor$. The prior because the Gamma Function is inaccurate behind zero, and the latter for the sake of accuracy.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18

fn gamma_function(n: f64) -> f64 {
    if n < 0f64 {
        panic!("The Gamma Function is not accurate for negative values!");
    }

    let standard_factorial = (1i64..=n.max(1f64) as i64)
        .fold(1f64, |acc, item| {
            acc * (item as f64 + n % 1f64)
        });

    if n % 1f64 == 0f64 {
        return standard_factorial as f64;
    }

    let decimal_factorial = definite_integral(|x| x.powf(n - n.floor()) * ((-x).exp()), 0f64, 1000f64, 0.0001f64);

    standard_factorial as f64 * decimal_factorial
}

Nice! If all is done correctly, you have a fully functional implementation of decimal factorials with the Gamma Function.

It’s important that you tweak the constants to your own application for accuracy that suits your needs. I will also note that we chose Left Riemann Sum as our method of approximation. While very easy to implement, it’s not accurate! A significantly more accurate method is Simpson’s Rule - with little extra computation cost.

Result -

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40

fn definite_integral(f: impl Fn(f64) -> f64, lower: f64, upper: f64, step: f64) -> f64 {
    let mut area: f64 = 0f64;

    let mut x: f64 = lower;
    let mut step: f64 = step;
    while x < upper {
        if (x + step) > upper {
            step = upper - x;
        }

        area += step * f(x);

        x += step;
    }

    area
}
fn gamma_function(n: f64) -> f64 {
    if n < 0f64 {
        panic!("The Gamma Function is not accurate for negative values!");
    }

    let standard_factorial = (1i64..=n.max(1f64) as i64)
        .fold(1f64, |acc, item| {
            acc * (item as f64 + n % 1f64)
        });

    if n % 1f64 == 0f64 {
        return standard_factorial as f64;
    }

    let decimal_factorial = definite_integral(|x| x.powf(n - n.floor()) * ((-x).exp()), 0f64, 1000f64, 0.0001f64);

    standard_factorial as f64 * decimal_factorial
}
fn main() {
    let answer = gamma_function(4.02f64);

    println!("4.02! = {answer}")
}

Outro -

Hopefully you found this article informative or helpful. If you did, follow my GitHub!

I often write these articles after struggling to find good documentation on a concept. That means there’s often a large-scale project on there related to this article, if you’d like to see it in practice! In this case, I needed a decimal factorial for my memory model project monikaiv2 - but didn’t want the overhead of adding a crate.