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
use crate::gammaln;
/// Beta function <math><mrow><mi>B</mi><mo>(</mo><mi>z</mi><mo>,</mo><mi>w</mi><mo>)</mo><mo>=</mo><mfrac><mrow><mi>Γ</mi><mo>(</mo><mi>z</mi><mo>)</mo><mi>Γ</mi><mo>(</mo><mi>w</mi><mo>)</mo></mrow><mrow><mi>Γ</mi><mo>(</mo><mi>z</mi><mo>+</mo><mi>w</mi><mo>)</mo></mrow></mfrac></mrow></math>.
///
/// This implementation evaluates the function using the natural logarithm of the Gamma function
/// to avoid intermediate overflow for large arguments:
/// <math><mrow><mi>B</mi><mo>(</mo><mi>z</mi><mo>,</mo><mi>w</mi><mo>)</mo><mo>=</mo><mi>exp</mi><mo>(</mo><mi>ln</mi><mi>Γ</mi><mo>(</mo><mi>z</mi><mo>)</mo><mo>+</mo><mi>ln</mi><mi>Γ</mi><mo>(</mo><mi>w</mi><mo>)</mo><mo>-</mo><mi>ln</mi><mi>Γ</mi><mo>(</mo><mi>z</mi><mo>+</mo><mi>w</mi><mo>)</mo><mo>)</mo></mrow></math>.
///
/// # Domain
/// - `z > 0`, `w > 0`
/// - Returns `NaN` if either argument is less than or equal to zero.
///
/// # Examples
/// ```
/// use abax::beta;
///
/// let result = beta(2.0, 3.0);
/// assert!((result - 1.0/12.0).abs() < 1e-15);
/// ```
pub fn beta(z: f64, w: f64) -> f64 {
if z <= 0.0 || w <= 0.0 {
return f64::NAN;
}
(gammaln(z) + gammaln(w) - gammaln(z + w)).exp()
}
#[cfg(test)]
mod tests {
use super::*;
use std::f64::consts::PI;
#[test]
fn test_beta_values() {
let tol = 1e-14;
// B(1, 1) = 1
assert!((beta(1.0, 1.0) - 1.0).abs() < tol);
// B(2, 1) = 1/2
assert!((beta(2.0, 1.0) - 0.5).abs() < tol);
// B(3, 2) = 1/12
assert!((beta(3.0, 2.0) - 1.0 / 12.0).abs() < tol);
// B(0.5, 0.5) = pi
assert!((beta(0.5, 0.5) - PI).abs() < tol);
}
#[test]
fn test_beta_symmetry() {
assert_eq!(beta(2.5, 4.2), beta(4.2, 2.5));
}
#[test]
fn test_beta_domain() {
assert!(beta(0.0, 1.0).is_nan());
assert!(beta(1.0, 0.0).is_nan());
assert!(beta(-0.5, 0.5).is_nan());
assert!(beta(0.5, -0.5).is_nan());
assert!(beta(f64::NAN, 1.0).is_nan());
}
}