rmathlib 1.1.0

A Rust Port of R's C Library of Special Functions
Documentation 
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
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59

use libm::log;
use libm::log1p;

use crate::gammafn;
use crate::lgammacor;
use crate::lgammafn;
use crate::nmath::ML_NAN;
use crate::nmath::ML_NEGINF;
use crate::nmath::ML_POSINF;
use crate::nmath::M_LN_SQRT_2PI;

/// This function returns the value of the log beta function
///
/// log B(a,b) = log G(a) + log G(b) - log G(a+b)
///
/// ## NOTES
///
/// This routine is a translation into C of a Fortran subroutine
/// by W. Fullerton of Los Alamos Scientific Laboratory.
pub fn lbeta(a: f64, b: f64) -> f64 {
    let (mut p, mut q) = (a, a);
    if b < p {
        p = b;
    }
    if b > q {
        q = b;
    }

    // both arguments must be >= 0
    if p < 0.0 {
        return ML_NAN;
    } else if p == 0.0 {
        return ML_POSINF;
    } else if !q.is_finite() {
        return ML_NEGINF;
    }

    if p >= 10.0 {
        // p and q are big.
        let corr = lgammacor(p) + lgammacor(q) - lgammacor(p + q);
        log(q) * -0.5
            + M_LN_SQRT_2PI
            + corr
            + (p - 0.5) * log(p / (p + q))
            + q * log1p(-p / (p + q))
    } else if q >= 10.0 {
        // p is small, but q is big.
        let corr = lgammacor(q) - lgammacor(p + q);
        lgammafn(p) + corr + p - p * log(p + q) + (q - 0.5) * log1p(-p / (p + q))
    } else {
        // p and q are small: p <= q < 10.
        // R change for very small args
        if p < 1e-306 {
            lgammafn(p) + (lgammafn(q) - lgammafn(p + q))
        } else {
            log(gammafn(p) * (gammafn(q) / gammafn(p + q)))
        }
    }
}