exp2

pub fn exp2(x: f64) callconv(.c) f64

File

Code

pub fn exp2(x: f64) callconv(.c) f64 {
    const tblsiz: u32 = @intCast(exp2dt.len / 2);
    const redux: f64 = 0x1.8p52 / @as(f64, @floatFromInt(tblsiz));
    const P1: f64 = 0x1.62e42fefa39efp-1;
    const P2: f64 = 0x1.ebfbdff82c575p-3;
    const P3: f64 = 0x1.c6b08d704a0a6p-5;
    const P4: f64 = 0x1.3b2ab88f70400p-7;
    const P5: f64 = 0x1.5d88003875c74p-10;

    const ux: u64 = @bitCast(x);
    const ix = @as(u32, @intCast(ux >> 32)) & 0x7FFFFFFF;

    // TODO: This should be handled beneath.
    if (math.isNan(x)) {
        return math.nan(f64);
    }

    // |x| >= 1022 or nan
    if (ix >= 0x408FF000) {
        // x >= 1024 or nan
        if (ix >= 0x40900000 and ux >> 63 == 0) {
            return if (compiler_rt.want_float_exceptions) x * 0x1p1023 else std.math.inf(f64);
        }
        // -inf or -nan
        if (ix >= 0x7FF00000) {
            return -1 / x;
        }
        // x <= -1022
        if (ux >> 63 != 0) {
            // underflow
            if (x <= -1075 or x - 0x1.0p52 + 0x1.0p52 != x) {
                if (compiler_rt.want_float_exceptions) mem.doNotOptimizeAway(@as(f32, @floatCast(-0x1.0p-149 / x)));
            }
            if (x <= -1075) {
                return 0;
            }
        }
    }
    // |x| < 0x1p-54
    else if (ix < 0x3C900000) {
        return 1.0 + x;
    }

    // NOTE: musl relies on unsafe behaviours which are replicated below
    // (addition overflow, division truncation, casting). Appears that this
    // produces the intended result but should confirm how GCC/Clang handle this
    // to ensure.

    // reduce x
    var uf: f64 = x + redux;
    // NOTE: musl performs an implicit 64-bit to 32-bit u32 truncation here
    var i_0: u32 = @truncate(@as(u64, @bitCast(uf)));
    i_0 +%= tblsiz / 2;

    const k: u32 = i_0 / tblsiz * tblsiz;
    const ik: i32 = @divTrunc(@as(i32, @bitCast(k)), tblsiz);
    i_0 %= tblsiz;
    uf -= redux;

    // r = exp2(y) = exp2t[i_0] * p(z - eps[i])
    var z: f64 = x - uf;
    const t: f64 = exp2dt[@intCast(2 * i_0)];
    z -= exp2dt[@intCast(2 * i_0 + 1)];
    const r: f64 = t + t * z * (P1 + z * (P2 + z * (P3 + z * (P4 + z * P5))));

    return math.scalbn(r, ik);
}