deep_time/math/
log.rs

1#![allow(clippy::indexing_slicing)]
2#![allow(clippy::excessive_precision)]
3#![allow(clippy::approx_constant)]
4#![allow(clippy::eq_op)]
5
6use crate::Real;
7
8const LN2_HI: Real = 6.93147180369123816490e-01; /* 3fe62e42 fee00000 */
9const LN2_LO: Real = 1.90821492927058770002e-10; /* 3dea39ef 35793c76 */
10const LG1: Real = 6.666666666666735130e-01; /* 3FE55555 55555593 */
11const LG2: Real = 3.999999999940941908e-01; /* 3FD99999 9997FA04 */
12const LG3: Real = 2.857142874366239149e-01; /* 3FD24924 94229359 */
13const LG4: Real = 2.222219843214978396e-01; /* 3FCC71C5 1D8E78AF */
14const LG5: Real = 1.818357216161805012e-01; /* 3FC74664 96CB03DE */
15const LG6: Real = 1.531383769920937332e-01; /* 3FC39A09 D078C69F */
16const LG7: Real = 1.479819860511658591e-01; /* 3FC2F112 DF3E5244 */
17
18/// The natural logarithm of `x` (Real).
19pub const fn log(mut x: Real) -> Real {
20    let x1p54 = Real::from_bits(0x4350000000000000); // 0x1p54 === 2 ^ 54
21
22    let mut ui = x.to_bits();
23    let mut hx: u32 = (ui >> 32) as u32;
24    let mut k: i32 = 0;
25
26    if (hx < 0x00100000) || ((hx >> 31) != 0) {
27        /* x < 2**-126  */
28        if ui << 1 == 0 {
29            return -1. / (x * x); /* log(+-0)=-inf */
30        }
31        if hx >> 31 != 0 {
32            return (x - x) / 0.0; /* log(-#) = NaN */
33        }
34        /* subnormal number, scale x up */
35        k -= 54;
36        x *= x1p54;
37        ui = x.to_bits();
38        hx = (ui >> 32) as u32;
39    } else if hx >= 0x7ff00000 {
40        return x;
41    } else if hx == 0x3ff00000 && ui << 32 == 0 {
42        return 0.;
43    }
44
45    /* reduce x into [sqrt(2)/2, sqrt(2)] */
46    hx += 0x3ff00000 - 0x3fe6a09e;
47    k += ((hx >> 20) as i32) - 0x3ff;
48    hx = (hx & 0x000fffff) + 0x3fe6a09e;
49    ui = ((hx as u64) << 32) | (ui & 0xffffffff);
50    x = Real::from_bits(ui);
51
52    let f: Real = x - 1.0;
53    let hfsq: Real = 0.5 * f * f;
54    let s: Real = f / (2.0 + f);
55    let z: Real = s * s;
56    let w: Real = z * z;
57    let t1: Real = w * (LG2 + w * (LG4 + w * LG6));
58    let t2: Real = z * (LG1 + w * (LG3 + w * (LG5 + w * LG7)));
59    let r: Real = t2 + t1;
60    let dk: Real = k as Real;
61    s * (hfsq + r) + dk * LN2_LO - hfsq + f + dk * LN2_HI
62}
deep_time/math/log.rs

deep_time/math/
log.rs
