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
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167

//! atanh(x) implementation.
//!
//! Small |x| uses the core-math polynomial with double-double correction.
//! For the remaining range, use log1p-difference or log-ratio to keep
//! <= 1 ULP accuracy with good throughput.

use super::atanh_data::{CH, CL};
use super::{asdouble, copysign, fasttwosum, fma_internal};
use super::{log::ln, log1p};

const TINY_BITS: u64 = 0x3e4d_12ed_0af1_a27f; // 0x1.d12ed0af1a27fp-27
const SMALL_BITS: u64 = 0x3fd0_0000_0000_0000; // 0x1p-2
const ONE_BITS: u64 = 0x3ff0_0000_0000_0000;
const INF_BITS: u64 = 0x7ff0_0000_0000_0000;
const LOG1P_BOUND: f64 = 0.5;

const SMALL_C: [f64; 9] = [
    f64::from_bits(0x3fc9_9999_9999_999a),
    f64::from_bits(0x3fc2_4924_9249_2244),
    f64::from_bits(0x3fbc_71c7_1c79_715f),
    f64::from_bits(0x3fb7_45d1_6f77_7723),
    f64::from_bits(0x3fb3_b13c_a417_4634),
    f64::from_bits(0x3fb1_10c9_7249_89bd),
    f64::from_bits(0x3fae_2d17_608a_5b2e),
    f64::from_bits(0x3faa_0b56_308c_ba0b),
    f64::from_bits(0x3faf_b634_1208_ad2e),
];

const SMALL_T: f64 = f64::from_bits(0x3c75_5555_5555_5555); // 0x1.5555555555555p-56
const SMALL_BASE: f64 = f64::from_bits(0x3fd5_5555_5555_5555); // 0x1.5555555555555p-2
const EPS_X4: f64 = f64::from_bits(0x3cad_0000_0000_0000); // 0x1.dp-53
const EPS_TINY: f64 = f64::from_bits(0x3980_0000_0000_0000); // 0x1p-103

#[inline(always)]
fn muldd_acc(xh: f64, xl: f64, ch: f64, cl: f64, l: &mut f64) -> f64 {
    let ahlh = ch * xl;
    let alhh = cl * xh;
    let ahhh = ch * xh;
    let mut ahhl = fma_internal(ch, xh, -ahhh);
    ahhl += alhh + ahlh;
    let chh = ahhh + ahhl;
    *l = (ahhh - chh) + ahhl;
    chh
}

#[inline(always)]
fn muldd_acc2(xh: f64, xl: f64, ch: f64, cl: f64, l: &mut f64) -> f64 {
    let ahlh = ch * xl;
    let alhh = cl * xh;
    let ahhh = ch * xh;
    let mut ahhl = fma_internal(ch, xh, -ahhh);
    ahhl += alhh + ahlh;
    fasttwosum(ahhh, ahhl, l)
}

#[inline(always)]
fn mulddd3(xh: f64, xl: f64, ch: f64, l: &mut f64) -> f64 {
    let hh = xh * ch;
    *l = fma_internal(ch, xh, -hh) + xl * ch;
    hh
}

#[inline(always)]
fn polydd3(xh: f64, xl: f64, n: usize, c: &[[f64; 2]], l: &mut f64) -> f64 {
    let mut i = n - 1;
    let mut cl = 0.0;
    let mut ch = fasttwosum(c[i][0], *l, &mut cl);
    cl += c[i][1];
    while i > 0 {
        i -= 1;
        ch = muldd_acc2(xh, xl, ch, cl, &mut cl);
        let mut tl = 0.0;
        let th = fasttwosum(c[i][0], ch, &mut tl);
        ch = th;
        cl += tl + c[i][1];
    }
    *l = cl;
    ch
}

#[inline(always)]
fn as_atanh_zero(x: f64) -> f64 {
    let x2 = x * x;
    let x2l = fma_internal(x, x, -x2);
    let y2 = x2 * (CL[0] + x2 * (CL[1] + x2 * (CL[2] + x2 * (CL[3] + x2 * (CL[4])))));
    let mut y2v = y2;
    let mut y1 = polydd3(x2, x2l, CH.len(), CH, &mut y2v);
    y1 = mulddd3(y1, y2v, x, &mut y2v);
    y1 = muldd_acc2(x2, x2l, y1, y2v, &mut y2v);
    let mut y1e = 0.0;
    let y0 = fasttwosum(x, y1, &mut y1e);
    let mut y2e = 0.0;
    let mut y1h = fasttwosum(y1e, y2v, &mut y2e);
    let mut t = (y1h).to_bits();
    if (t & (!0u64 >> 12)) == 0 {
        let w = (y2e).to_bits();
        if ((w ^ t) >> 63) != 0 {
            t = t.wrapping_sub(1);
        } else {
            t = t.wrapping_add(1);
        }
        y1h = asdouble(t);
    }
    y0 + y1h
}

#[inline(always)]
pub fn atanh(x: f64) -> f64 {
    let ax = x.abs();
    let aix = (ax).to_bits();
    if aix == 0 {
        return x;
    }
    if aix >= ONE_BITS {
        if aix == ONE_BITS {
            return if x.is_sign_negative() {
                f64::NEG_INFINITY
            } else {
                f64::INFINITY
            };
        }
        if aix > INF_BITS {
            return x + x;
        }
        return f64::NAN;
    }

    if aix < SMALL_BITS {
        if aix < TINY_BITS {
            return fma_internal(x, f64::from_bits(0x3c80_0000_0000_0000), x);
        }
        let x2 = x * x;
        let dx2 = fma_internal(x, x, -x2);
        let x4 = x2 * x2;
        let x3 = x2 * x;
        let x8 = x4 * x4;
        let dx3 = fma_internal(x2, x, -x3) + dx2 * x;
        let p = (SMALL_C[0] + x2 * SMALL_C[1])
            + x4 * (SMALL_C[2] + x2 * SMALL_C[3])
            + x8 * ((SMALL_C[4] + x2 * SMALL_C[5])
                + x4 * (SMALL_C[6] + x2 * SMALL_C[7])
                + x8 * SMALL_C[8]);
        let t = fma_internal(x2, p, SMALL_T);
        let mut pl = 0.0;
        let mut ph = fasttwosum(SMALL_BASE, t, &mut pl);
        ph = muldd_acc(ph, pl, x3, dx3, &mut pl);
        let mut tl = 0.0;
        ph = fasttwosum(x, ph, &mut tl);
        pl += tl;
        let eps = x * (x4 * EPS_X4 + EPS_TINY);
        let lb = ph + (pl - eps);
        let ub = ph + (pl + eps);
        if lb == ub {
            return lb;
        }
        return as_atanh_zero(x);
    }

    if ax <= LOG1P_BOUND {
        let y = 0.5 * (log1p(ax) - log1p(-ax));
        return copysign(y, x);
    }

    let ratio = (1.0 + ax) / (1.0 - ax);
    let y = 0.5 * ln(ratio);
    copysign(y, x)
}