Gamma function in Rust

Question

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.

Answer 1

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

1. 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.
2. 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).
3. 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.
4. If something is a constant, make it so. Constants use const or static and are named with SCREAMING_SNAKE_CASE.
5. There is inconsistent alignment in the lookup table. The first column is aligned on the decimal point, the second and third columns are not.
6. The sm variable doesn't need an explicit type as type inference will cover it.
7. There are inplace operations like *= which are shorter than foo = foo * bar.
8. There's a f64::recip method which may be more intuitive.
9. Operators like / should have space on both sides.
10. Is there a benefit to iterating in reverse every time? I'd just reorder the table and iterate forwards. Computers generally like going forwards.
11. 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:

1. 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.
2. Long numbers could be separated with underscores: 0.111_222_333.
3. 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.

• Dan: "SCREAMING_SNAKE_CASE" Googling this led to some interesting images
• Dair: Thank you for a lot of suggestions! I appreciate it a lot. I think I can confirm that the optimizer makes the optimization. Running in rust playground the LLVM IR generated for release has the line: %.op.op = fadd double %.op, -1.000000e+00. It then unrolls the fold and uses %op.op. I was also not aware that recip existed. Interesting language feature... (One small note: Maybe use TAYLOR_COEF instead of THE_TABLE)
• Turksarama: I didn't know about recip either, and frankly I'm not sure how useful it is. 1 / x isn't unintuitive and it definitely isn't clear at first glance that recip means 'reciprocal'.