logs.rs

1use std::f64::consts::LN_2;23const MANTISSA: u64 = 0x000f_ffff_ffff_ffff;4const ONE_EXPONENT: u64 = 0x3ff0_0000_0000_0000;5const EXPONENT: u64 = 0x7ff;6const BIAS: i32 = 1023;7const SCALE: f64 = 4_503_599_627_370_496.0;8const TERMS: usize = 16;910/// Returns the natural logarithm of a number, from a series instead of libm.11pub fn ln(x: f64) -> f64 {12    if x.is_nan() || x < 0.0 {13        return f64::NAN;14    }15    if x == 0.0 {16        return f64::NEG_INFINITY;17    }18    if x.is_infinite() {19        return x;20    }21    let mut bits = x.to_bits();22    let mut e: i32 = 0;23    if (bits >> 52) & EXPONENT == 0 {24        bits = (x * SCALE).to_bits();25        e -= 52;26    }27    e += (((bits >> 52) & EXPONENT) as i32) - BIAS;28    let mut m = f64::from_bits((bits & MANTISSA) | ONE_EXPONENT);29    if m >= 1.5 {30        m *= 0.5;31        e += 1;32    }33    let s = (m - 1.0) / (m + 1.0);34    let square = s * s;35    let mut term = s;36    let mut sum = 0.0;37    for k in 0..TERMS {38        sum += term / (2 * k + 1) as f64;39        term *= square;40    }41    e as f64 * LN_2 + 2.0 * sum42}4344#[cfg(test)]45mod tests {46    use super::*;4748    #[test]49    fn one_and_two_are_exact() {50        assert_eq!(ln(1.0).to_bits(), 0.0f64.to_bits());51        assert_eq!(ln(2.0), LN_2);52    }5354    #[test]55    fn the_edges_and_the_subnormals_hold() {56        assert!(ln(f64::NAN).is_nan());57        assert!(ln(-1.0).is_nan());58        assert_eq!(ln(0.0), f64::NEG_INFINITY);59        assert_eq!(ln(-0.0), f64::NEG_INFINITY);60        assert!(ln(f64::NEG_INFINITY).is_nan());61        assert_eq!(ln(f64::INFINITY), f64::INFINITY);62        assert_eq!(ln(f64::from_bits(1)), -1074.0 * LN_2);63    }6465    #[test]66    fn the_integers_track_libm() {67        for i in 2..=100_000u32 {68            let x = i as f64;69            let want = x.ln();70            assert!(71                (ln(x) - want).abs() <= 4.0 * f64::EPSILON * want.abs(),72                "ln({x}) is {} not {want}",73                ln(x)74            );75        }76    }7778    #[test]79    fn the_powers_of_two_track_libm() {80        for k in -300..=300i32 {81            let x = f64::from_bits(((BIAS + k) as u64) << 52);82            let want = x.ln();83            assert!(84                (ln(x) - want).abs() <= 4.0 * f64::EPSILON * want.abs(),85                "ln(2^{k}) is {} not {want}",86                ln(x)87            );88        }89    }90}