vendor/rand-0.5.6/src/distributions/log_gamma.rs

   1 // Copyright 2016-2017 The Rust Project Developers. See the COPYRIGHT
   2 // file at the top-level directory of this distribution and at
   3 // https://rust-lang.org/COPYRIGHT.
   4 //
   5 // Licensed under the Apache License, Version 2.0 <LICENSE-APACHE or
   6 // https://www.apache.org/licenses/LICENSE-2.0> or the MIT license
   7 // <LICENSE-MIT or https://opensource.org/licenses/MIT>, at your
   8 // option. This file may not be copied, modified, or distributed
   9 // except according to those terms.
  10
  11 /// Calculates ln(gamma(x)) (natural logarithm of the gamma
  12 /// function) using the Lanczos approximation.
  13 ///
  14 /// The approximation expresses the gamma function as:
  15 /// `gamma(z+1) = sqrt(2*pi)*(z+g+0.5)^(z+0.5)*exp(-z-g-0.5)*Ag(z)`
  16 /// `g` is an arbitrary constant; we use the approximation with `g=5`.
  17 ///
  18 /// Noting that `gamma(z+1) = z*gamma(z)` and applying `ln` to both sides:
  19 /// `ln(gamma(z)) = (z+0.5)*ln(z+g+0.5)-(z+g+0.5) + ln(sqrt(2*pi)*Ag(z)/z)`
  20 ///
  21 /// `Ag(z)` is an infinite series with coefficients that can be calculated
  22 /// ahead of time - we use just the first 6 terms, which is good enough
  23 /// for most purposes.
  24 pub fn log_gamma(x: f64) -> f64 {
  25     // precalculated 6 coefficients for the first 6 terms of the series
  26     let coefficients: [f64; 6] = [
  27         76.18009172947146,
  28         -86.50532032941677,
  29         24.01409824083091,
  30         -1.231739572450155,
  31         0.1208650973866179e-2,
  32         -0.5395239384953e-5,
  33     ];
  34
  35     // (x+0.5)*ln(x+g+0.5)-(x+g+0.5)
  36     let tmp = x + 5.5;
  37     let log = (x + 0.5) * tmp.ln() - tmp;
  38
  39     // the first few terms of the series for Ag(x)
  40     let mut a = 1.000000000190015;
  41     let mut denom = x;
  42     for coeff in &coefficients {
  43         denom += 1.0;
  44         a += coeff / denom;
  45     }
  46
  47     // get everything together
  48     // a is Ag(x)
  49     // 2.5066... is sqrt(2pi)
  50     log + (2.5066282746310005 * a / x).ln()
  51 }