fastmaths 0.1.0

High-performance no_std math routines targeting glibc parity and strict f64 accuracy.
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
53
54
55
56
57
58
59
60
61
62
63

//! log10(x) implementation.
//!
//! Port of glibc/fdlibm: argument reduction via ilogb and
//! log10(x) = y*log10_2hi + (y*log10_2lo + ivln10*ln(mantissa)),
//! using high/low splits for ≤1 ULP accuracy.

use super::{TWO54, fma_internal, log::ln, log::ln_dd, two_sum};

const IVLN10: f64 = f64::from_bits(0x3fdb_cb7b_1526_e50e);
const IVLN10_HI: f64 = f64::from_bits(0x3fdb_cb7b_1800_0000);
const IVLN10_LO: f64 = f64::from_bits(0xbe26_c8d7_9000_0000);
const LOG10_2_HI: f64 = f64::from_bits(0x3fd3_4413_509f_6000);
const LOG10_2_LO: f64 = f64::from_bits(0x3d59_fef3_11f1_2b36);

#[inline]
pub fn log10(x: f64) -> f64 {
    let mut ux = x.to_bits();
    let ax = ux & 0x7fff_ffff_ffff_ffff;
    if ax == 0 {
        return f64::NEG_INFINITY;
    }
    if (ux >> 63) != 0 {
        return f64::NAN;
    }
    if ax >= 0x7ff0_0000_0000_0000 {
        return x + x;
    }

    let r = x - 1.0;
    let ar = r.abs();
    if ar <= 0.4 {
        let (lnh, lnl) = ln_dd(1.0 + r);
        let p_hi = IVLN10_HI * lnh;
        let p_lo = fma_internal(IVLN10_HI, lnh, -p_hi) + IVLN10_LO * lnh + IVLN10 * lnl;
        let (s, e) = two_sum(p_hi, p_lo);
        return s + e;
    }

    let mut k: i32 = 0;
    if ax < 0x0010_0000_0000_0000 {
        k -= 54;
        let y = x * TWO54;
        ux = y.to_bits();
    }

    let exp = ((ux >> 52) & 0x7ff) as i32;
    k += exp - 1023;
    let i = if k < 0 { 1 } else { 0 };
    let mant = (ux & 0x000f_ffff_ffff_ffff) | (((0x3ff - i) as u64) << 52);
    let y = (k + i) as f64;
    let m = f64::from_bits(mant);
    let yhi = y * LOG10_2_HI;
    let ylo = y * LOG10_2_LO;
    let (lnh, lnl) = if k == -1 || k == 0 {
        ln_dd(m)
    } else {
        (ln(m), 0.0)
    };
    let p_hi = IVLN10_HI * lnh;
    let p_lo = fma_internal(IVLN10_HI, lnh, -p_hi) + IVLN10_LO * lnh + IVLN10 * lnl;
    let (s, e) = two_sum(yhi, p_hi);
    s + (e + ylo + p_lo)
}