tract_linalg/generic/
ln.rs

1#![allow(clippy::excessive_precision)]
2
3/// Coefficients of the minimax fit of `ln(1 + f)`, less its first two Taylor terms, over
4/// `[SPLIT / 2 - 1, SPLIT - 1]`, highest power first.
5///
6/// Shared by every ln kernel whatever its width: the fit is one decision, and a kernel
7/// keeping its own copy of it could come to disagree with the rest silently.
8pub const POLY: [f32; 9] = [
9    7.0376836292e-2,
10    -1.1514610310e-1,
11    1.1676998740e-1,
12    -1.2420140846e-1,
13    1.4249322787e-1,
14    -1.6668057665e-1,
15    2.0000714765e-1,
16    -2.4999993993e-1,
17    3.3333331174e-1,
18];
19
20/// `ln 2`, split so that `e * LN2_HI` is exact for every exponent an f32 carries and
21/// `LN2_LO` holds what it drops. One constant instead costs the mantissa its last bits
22/// wherever `|e|` is large.
23pub const LN2_HI: f32 = 0.693_359_375;
24pub const LN2_LO: f32 = -2.121_944_4e-4;
25
26/// Where the mantissa splits: above it the fit's argument would leave the interval
27/// [`POLY`] was fitted over, so a power of two moves into the exponent instead.
28pub const SPLIT: f32 = std::f32::consts::SQRT_2;
29
30/// `2^SUBNORMAL_SHIFT`, what a subnormal input is scaled by before its exponent field is
31/// read, the exponent then paying the shift back. A subnormal carries no implicit leading
32/// one, so its field says nothing about its magnitude until it is normal.
33pub const SUBNORMAL_SCALE: f32 = 16_777_216.0;
34pub const SUBNORMAL_SHIFT: i32 = 24;
35
36/// f32 natural logarithm: the exponent, plus the [`POLY`] fit of the mantissa reduced to
37/// the interval around one where the fit holds.
38///
39/// Within one ulp of a correctly rounded `ln` over the whole f32 range, subnormals
40/// included. The specials are IEEE's: `±0` gives `-inf`, a negative gives NaN, `+inf`
41/// gives `+inf`, and a NaN input gives a NaN of the kernel's own rather than its payload
42/// back.
43pub fn sln(x: f32) -> f32 {
44    let subnormal = x < f32::MIN_POSITIVE;
45    let scaled = if subnormal { x * SUBNORMAL_SCALE } else { x };
46    let bits = scaled.to_bits();
47    let mut e = ((bits >> 23) & 0xff) as i32 - 127;
48    if subnormal {
49        e -= SUBNORMAL_SHIFT;
50    }
51    let mut m = f32::from_bits((bits & 0x007fffff) | 0x3f800000);
52    if m > SPLIT {
53        m *= 0.5;
54        e += 1;
55    }
56    let f = m - 1.0;
57    let f2 = f * f;
58    let mut p = POLY[0];
59    for c in &POLY[1..] {
60        p = p.mul_add(f, *c);
61    }
62    let e = e as f32;
63    let y = p * f2 * f;
64    let y = e.mul_add(LN2_LO, y);
65    let y = (-0.5f32).mul_add(f2, y) + f;
66    let y = e.mul_add(LN2_HI, y);
67    if x <= 0.0 || x.is_nan() {
68        if x == 0.0 { f32::NEG_INFINITY } else { f32::NAN }
69    } else if x.is_infinite() {
70        f32::INFINITY
71    } else {
72        y
73    }
74}
75
76routine_ew_rust!(generic;
77    f32,
78    generic_ln_f32_4n,
79    4,
80    4,
81    fn run(x: &mut [f32], _: ()) {
82        debug_assert!(x.len() % Self::nr() == 0);
83        debug_assert!(x.as_ptr() as usize % Self::alignment_bytes() == 0);
84        x.iter_mut().for_each(|px| *px = sln(*px))
85    },
86    func(Ln)
87);
tract_linalg/generic/ln.rs

tract_linalg/generic/
ln.rs
