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

use libm::ldexp;

use crate::nmath::*;
use crate::rmath::*;

/// Compute the density of the normal distribution.
pub fn dnorm4(x: f64, mu: f64, sigma: f64, give_log: bool) -> f64 {
    if x.is_nan() || mu.is_nan() || sigma.is_nan() {
        return x + mu + sigma;
    }
    if sigma < 0.0 {
        return ML_NAN; // Placeholder for warning
    }
    if !sigma.is_finite() {
        return 0.0;
    }
    if !x.is_finite() && mu == x {
        return ML_NAN;
    }
    if sigma == 0.0 {
        return if x == mu { ML_POSINF } else { 0.0 };
    }

    let x = (x - mu) / sigma;

    if !x.is_finite() {
        return 0.0;
    }

    let x = x.abs();
    if x >= 2.0 * f64::MAX.sqrt() {
        return 0.0;
    }
    if give_log {
        return -(M_LN_SQRT_2PI + 0.5 * x * x + sigma.ln());
    }

    // Following block is for the else case of MATHLIB_FAST_dnorm
    if x < 5.0 {
        M_1_SQRT_2PI * (-0.5 * x * x).exp() / sigma
    } else {
        let bound =
            (-2.0 * M_LN2 * (f64::MIN_EXP as f64 + 1.0 - f64::MANTISSA_DIGITS as f64)).sqrt();
        if x > bound {
            0.0
        } else {
            let x1 = ldexp(r_forceint(ldexp(x, 16)), -16);
            let x2 = x - x1;
            M_1_SQRT_2PI / sigma * (-0.5 * x1 * x1).exp() * ((-0.5 * x2 - x1) * x2).exp()
        }
    }
}