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

//! log1p(x) implementation.
//!
//! For tiny x uses series/polynomial to avoid cancellation; for larger x uses
//! log(1+x) with compensated argument reduction. Coefficients are derived from
//! glibc/fdlibm minimax fits.

use super::{LN2_HI, LN2_LO, TWO54, hi_word, lo_word, with_hi_lo};
const LP: [f64; 8] = [
    0.0,
    6.666_666_666_666_735_130e-01,
    3.999_999_999_940_941_908e-01,
    2.857_142_874_366_239_149e-01,
    2.222_219_843_214_978_396e-01,
    1.818_357_216_161_805_012e-01,
    1.531_383_769_920_937_332e-01,
    1.479_819_860_511_658_591e-01,
];

#[inline]
pub fn log1p(x: f64) -> f64 {
    let hx = hi_word(x) as i32;
    let ax = (hx & 0x7fff_ffff) as u32;

    let mut k = 1i32;
    let mut f = 0.0;
    let mut c = 0.0;
    let mut hu: u32 = 0;
    if hx < 0x3fda_827a {
        if ax >= 0x3ff0_0000 {
            if x == -1.0 {
                return f64::NEG_INFINITY;
            }
            return f64::NAN;
        }
        if ax < 0x3e20_0000 {
            let _ = TWO54 + x;
            if ax < 0x3c90_0000 {
                return x;
            }
            return x - x * x * 0.5;
        }
        if hx > 0 || hx <= 0xbfd2_bec3u32 as i32 {
            k = 0;
            f = x;
            hu = 1;
        }
    } else if ax >= 0x7ff0_0000 {
        return x + x;
    }

    if k != 0 {
        if hx < 0x4340_0000 {
            let mut u = 1.0 + x;
            hu = hi_word(u);
            k = ((hu >> 20) & 0x7ff) as i32 - 1023;
            c = if k > 0 { 1.0 - (u - x) } else { x - (u - 1.0) };
            c /= u;
            hu &= 0x000f_ffff;
            if hu < 0x6a09e {
                u = with_hi_lo(hu | 0x3ff0_0000, lo_word(u));
            } else {
                k += 1;
                u = with_hi_lo(hu | 0x3fe0_0000, lo_word(u));
                hu = (0x0010_0000 - hu) >> 2;
            }
            f = u - 1.0;
        } else {
            let mut u = x;
            hu = hi_word(u);
            k = ((hu >> 20) & 0x7ff) as i32 - 1023;
            c = 0.0;
            hu &= 0x000f_ffff;
            if hu < 0x6a09e {
                u = with_hi_lo(hu | 0x3ff0_0000, lo_word(u));
            } else {
                k += 1;
                u = with_hi_lo(hu | 0x3fe0_0000, lo_word(u));
                hu = (0x0010_0000 - hu) >> 2;
            }
            f = u - 1.0;
        }
    }

    let hfsq = 0.5 * f * f;
    if hu == 0 {
        if f == 0.0 {
            if k == 0 {
                return 0.0;
            }
            c += (k as f64) * LN2_LO;
            return (k as f64) * LN2_HI + c;
        }
        let r = hfsq * (1.0 - 0.666_666_666_666_666_66 * f);
        if k == 0 {
            return f - r;
        }
        return (k as f64) * LN2_HI - ((r - ((k as f64) * LN2_LO + c)) - f);
    }

    let s = f / (2.0 + f);
    let z = s * s;
    let r1 = z * LP[1];
    let z2 = z * z;
    let r2 = LP[2] + z * LP[3];
    let z4 = z2 * z2;
    let r3 = LP[4] + z * LP[5];
    let z6 = z4 * z2;
    let r4 = LP[6] + z * LP[7];
    let r = r1 + z2 * r2 + z4 * r3 + z6 * r4;
    if k == 0 {
        return f - (hfsq - s * (hfsq + r));
    }
    (k as f64) * LN2_HI - ((hfsq - (s * (hfsq + r) + ((k as f64) * LN2_LO + c))) - f)
}