Gamma function in Rust

The gamma function is one of a couple nice continuous extensions to the traditional factorial function. I used this Python program as a reference, which in turn, uses this Ada program. As the Ada program describes:


The implementation uses Taylor series coefficients of $$Γ(x+1)^{-1}, |x| < \infty$$ The coefficients are taken from Mathematical functions and their approximations by Yudell L. Luke.

Here is the Rust translation (Rust Playground):

fn gamma(x: f64) -> f64 {
    let _a: [f64; 29] = [ 1.00000000000000000000, 0.57721566490153286061, -0.65587807152025388108,
                         -0.04200263503409523553, 0.16653861138229148950, -0.04219773455554433675,
                         -0.00962197152787697356, 0.00721894324666309954, -0.00116516759185906511,
                         -0.00021524167411495097, 0.00012805028238811619, -0.00002013485478078824,
                         -0.00000125049348214267, 0.00000113302723198170, -0.00000020563384169776,
                          0.00000000611609510448, 0.00000000500200764447, -0.00000000118127457049,
                          0.00000000010434267117, 0.00000000000778226344, -0.00000000000369680562,
                          0.00000000000051003703, -0.00000000000002058326, -0.00000000000000534812,
                          0.00000000000000122678, -0.00000000000000011813, 0.00000000000000000119,
                          0.00000000000000000141, -0.00000000000000000023 ];
    let mut sm: f64 = 0.00000000000000000002;
    for an in _a.iter().rev() {
        sm = sm * (x - 1.0) + an;
    }
    1.0 / sm
}

fn main() {
    for i in 1..11 {
        let f = i as f64;
        println!("{}", gamma(f/3.0));
    }
}

Any suggestions on making this code better such as making the code more Rust-idiomatic would be appreciated.

There's not a lot of code here so there's not a lot to say. ^_^

• Implementations like this are generally outside of the day-to-day knowledge of most programmers. You may wish to include a reference in the code to how the algorithm was derived and how it works.

• The chosen variable names are pretty useless. I have no idea what an or sm mean. Because I don't know what they mean, I didn't change them, but you should. (Looking back at the Ada implementation, I see that sm should be sum and an is just the current value of the lookup table, but it shouldn't be required to find 2 other implementations to understand this one).

• Leading underscores in a variable name means the variable is unused (but cannot be removed for some reason). Don't use a underscore and use the variable at the same time.

• If something is a constant, make it so. Constants use const or static and are named with SCREAMING_SNAKE_CASE.

• There is inconsistent alignment in the lookup table. The first column is aligned on the decimal point, the second and third columns are not.

• The sm variable doesn't need an explicit type as type inference will cover it.

• There are inplace operations like = which are shorter than foo = foo bar.

• There's a f64::recip method which may be more intuitive.

• Operators like / should have space on both sides.

• Is there a benefit to iterating in reverse every time? I'd just reorder the table and iterate forwards. Computers generally like going forwards.

• The loop can be simplified by using Iterator::fold.

const TAYLOR_COEFFICIENTS: [f64; 29] = [
    -0.00000000000000000023,  0.00000000000000000141,  0.00000000000000000119,
    -0.00000000000000011813,  0.00000000000000122678, -0.00000000000000534812,
    -0.00000000000002058326,  0.00000000000051003703, -0.00000000000369680562,
     0.00000000000778226344,  0.00000000010434267117, -0.00000000118127457049,
     0.00000000500200764447,  0.00000000611609510448, -0.00000020563384169776,
     0.00000113302723198170, -0.00000125049348214267, -0.00002013485478078824,
     0.00012805028238811619, -0.00021524167411495097, -0.00116516759185906511,
     0.00721894324666309954, -0.00962197152787697356, -0.04219773455554433675,
     0.16653861138229148950, -0.04200263503409523553, -0.65587807152025388108,
     0.57721566490153286061,  1.00000000000000000000,
];

const INITIAL_SUM: f64 = 0.00000000000000000002;

fn gamma(x: f64) -> f64 {
    TAYLOR_COEFFICIENTS.iter().fold(INITIAL_SUM, |sum, coefficient| {
        sum * (x - 1.0) + coefficient
    }).recip()
}

fn main() {
    for i in 1..11 {
        let f = i as f64;
        println!("{}", gamma(f / 3.0));
    }
}

Notes based on the original Ada implementation:

• The lookup table was defined inside the gamma function, and the same could be done here. It's a matter of personal preference at this point in time.

• Long numbers could be separated with underscores: 0.111_222_333.

• The loop invariant 1.0 - x could be hoisted out of the loop, but I'd trust the optimizer to do that until proven otherwise.

StackExchange Code Review Q#116850, answer score: 10