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

//! asinh(x) implementation.
//!
//! Piecewise algorithm: |x| small -> x; medium -> log1p(x + x^2/(1+sqrt(1+x^2)));
//! large -> log(2x). Designed to avoid cancellation and overflow.

use super::{fma_internal, ln, log1p, sqrt, two_sum};

#[inline(always)]
pub fn asinh(x: f64) -> f64 {
    let ux = x.to_bits();
    let e = ((ux >> 52) & 0x7ff) as i32;
    let sign = (ux >> 63) != 0;
    let mut ax = f64::from_bits(ux & 0x7fff_ffff_ffff_ffffu64);

    if e >= 0x3ff + 26 {
        // |x| >= 2^26 or inf/nan
        ax = ln(ax) + core::f64::consts::LN_2;
    } else if e > 0x3ff {
        // |x| >= 2
        ax = ln(2.0 * ax + 1.0 / (sqrt(ax * ax + 1.0) + ax));
    } else if e >= 0x3ff - 26 {
        // |x| >= 2^-26
        let z = fma_internal(ax, ax, 0.0);
        let s = sqrt(z + 1.0);
        let t = z / (s + 1.0);
        let (sum_hi, sum_lo) = two_sum(ax, t);
        ax = log1p(sum_hi) + sum_lo / (1.0 + sum_hi);
    } else {
        // tiny
        return x;
    }

    if sign { -ax } else { ax }
}