crabc-libc 0.7.6

// Translated from musl: exp.c, exp2.c, expf.c, exp_data.c, exp2f_data.c

// ============================================================
// Data tables and constants for exp (double precision)
// ============================================================

const EXP_TABLE_BITS: u32 = 7;
const EXP_N: u64 = 1 << EXP_TABLE_BITS;

const EXP_INVLN2N: f64 = asdouble(0x40671547652B82FE);
const EXP_NEGLN2HI_N: f64 = asdouble(0xBF762E42FEFA0000);
const EXP_NEGLN2LO_N: f64 = asdouble(0xBD0CF79ABC9E3B3A);
const EXP_SHIFT: f64 = asdouble(0x4338000000000000);

const EXP_C2: f64 = asdouble(0x3FDFFFFFFFFFFDBD);
const EXP_C3: f64 = asdouble(0x3FC555555555543C);
const EXP_C4: f64 = asdouble(0x3FA55555CF172B91);
const EXP_C5: f64 = asdouble(0x3F81111167A4D017);

const EXP_1P1009: f64 = asdouble(0x7F00000000000000);
const EXP_1P_NEG1022: f64 = asdouble(0x0010000000000000);

// 2^(k/N) ~= H[k]*(1 + T[k]) for k in [0, N)
// tab[2*k] = asuint64(T[k]), tab[2*k+1] = asuint64(H[k]) - (k << 52)/N
const EXP_TAB: [u64; 256] = [
    0x0, 0x3ff0000000000000,
    0x3c9b3b4f1a88bf6e, 0x3feff63da9fb3335,
    0xbc7160139cd8dc5d, 0x3fefec9a3e778061,
    0xbc905e7a108766d1, 0x3fefe315e86e7f85,
    0x3c8cd2523567f613, 0x3fefd9b0d3158574,
    0xbc8bce8023f98efa, 0x3fefd06b29ddf6de,
    0x3c60f74e61e6c861, 0x3fefc74518759bc8,
    0x3c90a3e45b33d399, 0x3fefbe3ecac6f383,
    0x3c979aa65d837b6d, 0x3fefb5586cf9890f,
    0x3c8eb51a92fdeffc, 0x3fefac922b7247f7,
    0x3c3ebe3d702f9cd1, 0x3fefa3ec32d3d1a2,
    0xbc6a033489906e0b, 0x3fef9b66affed31b,
    0xbc9556522a2fbd0e, 0x3fef9301d0125b51,
    0xbc5080ef8c4eea55, 0x3fef8abdc06c31cc,
    0xbc91c923b9d5f416, 0x3fef829aaea92de0,
    0x3c80d3e3e95c55af, 0x3fef7a98c8a58e51,
    0xbc801b15eaa59348, 0x3fef72b83c7d517b,
    0xbc8f1ff055de323d, 0x3fef6af9388c8dea,
    0x3c8b898c3f1353bf, 0x3fef635beb6fcb75,
    0xbc96d99c7611eb26, 0x3fef5be084045cd4,
    0x3c9aecf73e3a2f60, 0x3fef54873168b9aa,
    0xbc8fe782cb86389d, 0x3fef4d5022fcd91d,
    0x3c8a6f4144a6c38d, 0x3fef463b88628cd6,
    0x3c807a05b0e4047d, 0x3fef3f49917ddc96,
    0x3c968efde3a8a894, 0x3fef387a6e756238,
    0x3c875e18f274487d, 0x3fef31ce4fb2a63f,
    0x3c80472b981fe7f2, 0x3fef2b4565e27cdd,
    0xbc96b87b3f71085e, 0x3fef24dfe1f56381,
    0x3c82f7e16d09ab31, 0x3fef1e9df51fdee1,
    0xbc3d219b1a6fbffa, 0x3fef187fd0dad990,
    0x3c8b3782720c0ab4, 0x3fef1285a6e4030b,
    0x3c6e149289cecb8f, 0x3fef0cafa93e2f56,
    0x3c834d754db0abb6, 0x3fef06fe0a31b715,
    0x3c864201e2ac744c, 0x3fef0170fc4cd831,
    0x3c8fdd395dd3f84a, 0x3feefc08b26416ff,
    0xbc86a3803b8e5b04, 0x3feef6c55f929ff1,
    0xbc924aedcc4b5068, 0x3feef1a7373aa9cb,
    0xbc9907f81b512d8e, 0x3feeecae6d05d866,
    0xbc71d1e83e9436d2, 0x3feee7db34e59ff7,
    0xbc991919b3ce1b15, 0x3feee32dc313a8e5,
    0x3c859f48a72a4c6d, 0x3feedea64c123422,
    0xbc9312607a28698a, 0x3feeda4504ac801c,
    0xbc58a78f4817895b, 0x3feed60a21f72e2a,
    0xbc7c2c9b67499a1b, 0x3feed1f5d950a897,
    0x3c4363ed60c2ac11, 0x3feece086061892d,
    0x3c9666093b0664ef, 0x3feeca41ed1d0057,
    0x3c6ecce1daa10379, 0x3feec6a2b5c13cd0,
    0x3c93ff8e3f0f1230, 0x3feec32af0d7d3de,
    0x3c7690cebb7aafb0, 0x3feebfdad5362a27,
    0x3c931dbdeb54e077, 0x3feebcb299fddd0d,
    0xbc8f94340071a38e, 0x3feeb9b2769d2ca7,
    0xbc87deccdc93a349, 0x3feeb6daa2cf6642,
    0xbc78dec6bd0f385f, 0x3feeb42b569d4f82,
    0xbc861246ec7b5cf6, 0x3feeb1a4ca5d920f,
    0x3c93350518fdd78e, 0x3feeaf4736b527da,
    0x3c7b98b72f8a9b05, 0x3feead12d497c7fd,
    0x3c9063e1e21c5409, 0x3feeab07dd485429,
    0x3c34c7855019c6ea, 0x3feea9268a5946b7,
    0x3c9432e62b64c035, 0x3feea76f15ad2148,
    0xbc8ce44a6199769f, 0x3feea5e1b976dc09,
    0xbc8c33c53bef4da8, 0x3feea47eb03a5585,
    0xbc845378892be9ae, 0x3feea34634ccc320,
    0xbc93cedd78565858, 0x3feea23882552225,
    0x3c5710aa807e1964, 0x3feea155d44ca973,
    0xbc93b3efbf5e2228, 0x3feea09e667f3bcd,
    0xbc6a12ad8734b982, 0x3feea012750bdabf,
    0xbc6367efb86da9ee, 0x3fee9fb23c651a2f,
    0xbc80dc3d54e08851, 0x3fee9f7df9519484,
    0xbc781f647e5a3ecf, 0x3fee9f75e8ec5f74,
    0xbc86ee4ac08b7db0, 0x3fee9f9a48a58174,
    0xbc8619321e55e68a, 0x3fee9feb564267c9,
    0x3c909ccb5e09d4d3, 0x3feea0694fde5d3f,
    0xbc7b32dcb94da51d, 0x3feea11473eb0187,
    0x3c94ecfd5467c06b, 0x3feea1ed0130c132,
    0x3c65ebe1abd66c55, 0x3feea2f336cf4e62,
    0xbc88a1c52fb3cf42, 0x3feea427543e1a12,
    0xbc9369b6f13b3734, 0x3feea589994cce13,
    0xbc805e843a19ff1e, 0x3feea71a4623c7ad,
    0xbc94d450d872576e, 0x3feea8d99b4492ed,
    0x3c90ad675b0e8a00, 0x3feeaac7d98a6699,
    0x3c8db72fc1f0eab4, 0x3feeace5422aa0db,
    0xbc65b6609cc5e7ff, 0x3feeaf3216b5448c,
    0x3c7bf68359f35f44, 0x3feeb1ae99157736,
    0xbc93091fa71e3d83, 0x3feeb45b0b91ffc6,
    0xbc5da9b88b6c1e29, 0x3feeb737b0cdc5e5,
    0xbc6c23f97c90b959, 0x3feeba44cbc8520f,
    0xbc92434322f4f9aa, 0x3feebd829fde4e50,
    0xbc85ca6cd7668e4b, 0x3feec0f170ca07ba,
    0x3c71affc2b91ce27, 0x3feec49182a3f090,
    0x3c6dd235e10a73bb, 0x3feec86319e32323,
    0xbc87c50422622263, 0x3feecc667b5de565,
    0x3c8b1c86e3e231d5, 0x3feed09bec4a2d33,
    0xbc91bbd1d3bcbb15, 0x3feed503b23e255d,
    0x3c90cc319cee31d2, 0x3feed99e1330b358,
    0x3c8469846e735ab3, 0x3feede6b5579fdbf,
    0xbc82dfcd978e9db4, 0x3feee36bbfd3f37a,
    0x3c8c1a7792cb3387, 0x3feee89f995ad3ad,
    0xbc907b8f4ad1d9fa, 0x3feeee07298db666,
    0xbc55c3d956dcaeba, 0x3feef3a2b84f15fb,
    0xbc90a40e3da6f640, 0x3feef9728de5593a,
    0xbc68d6f438ad9334, 0x3feeff76f2fb5e47,
    0xbc91eee26b588a35, 0x3fef05b030a1064a,
    0x3c74ffd70a5fddcd, 0x3fef0c1e904bc1d2,
    0xbc91bdfbfa9298ac, 0x3fef12c25bd71e09,
    0x3c736eae30af0cb3, 0x3fef199bdd85529c,
    0x3c8ee3325c9ffd94, 0x3fef20ab5fffd07a,
    0x3c84e08fd10959ac, 0x3fef27f12e57d14b,
    0x3c63cdaf384e1a67, 0x3fef2f6d9406e7b5,
    0x3c676b2c6c921968, 0x3fef3720dcef9069,
    0xbc808a1883ccb5d2, 0x3fef3f0b555dc3fa,
    0xbc8fad5d3ffffa6f, 0x3fef472d4a07897c,
    0xbc900dae3875a949, 0x3fef4f87080d89f2,
    0x3c74a385a63d07a7, 0x3fef5818dcfba487,
    0xbc82919e2040220f, 0x3fef60e316c98398,
    0x3c8e5a50d5c192ac, 0x3fef69e603db3285,
    0x3c843a59ac016b4b, 0x3fef7321f301b460,
    0xbc82d52107b43e1f, 0x3fef7c97337b9b5f,
    0xbc892ab93b470dc9, 0x3fef864614f5a129,
    0x3c74b604603a88d3, 0x3fef902ee78b3ff6,
    0x3c83c5ec519d7271, 0x3fef9a51fbc74c83,
    0xbc8ff7128fd391f0, 0x3fefa4afa2a490da,
    0xbc8dae98e223747d, 0x3fefaf482d8e67f1,
    0x3c8ec3bc41aa2008, 0x3fefba1bee615a27,
    0x3c842b94c3a9eb32, 0x3fefc52b376bba97,
    0x3c8a64a931d185ee, 0x3fefd0765b6e4540,
    0xbc8e37bae43be3ed, 0x3fefdbfdad9cbe14,
    0x3c77893b4d91cd9d, 0x3fefe7c1819e90d8,
    0x3c5305c14160cc89, 0x3feff3c22b8f71f1,
];

// ============================================================
// Data tables and constants for expf (single precision)
// ============================================================

const EXP2F_TABLE_BITS: u32 = 5;
const EXP2F_N: u64 = 1 << EXP2F_TABLE_BITS;

const EXP2F_SHIFT_SCALED: f64 = asdouble(0x42E8000000000000);
const EXP2F_POLY: [f64; 3] = [
    asdouble(0x3FAC6AF84B912394),
    asdouble(0x3FCEBFCE50FAC4F3),
    asdouble(0x3FE62E42FF0C52D6),
];
const EXP2F_SHIFT: f64 = asdouble(0x4338000000000000);
const EXP2F_INVLN2_SCALED: f64 = asdouble(0x40471547652B82FE);
const EXP2F_POLY_SCALED: [f64; 3] = [
    asdouble(0x3EBC6AF84B912394),
    asdouble(0x3F2EBFCE50FAC4F3),
    asdouble(0x3F962E42FF0C52D6),
];

// tab[i] = uint(2^(i/N)) - (i << 52-BITS)
const EXP2F_TAB: [u64; 32] = [
    0x3ff0000000000000, 0x3fefd9b0d3158574, 0x3fefb5586cf9890f, 0x3fef9301d0125b51,
    0x3fef72b83c7d517b, 0x3fef54873168b9aa, 0x3fef387a6e756238, 0x3fef1e9df51fdee1,
    0x3fef06fe0a31b715, 0x3feef1a7373aa9cb, 0x3feedea64c123422, 0x3feece086061892d,
    0x3feebfdad5362a27, 0x3feeb42b569d4f82, 0x3feeab07dd485429, 0x3feea47eb03a5585,
    0x3feea09e667f3bcd, 0x3fee9f75e8ec5f74, 0x3feea11473eb0187, 0x3feea589994cce13,
    0x3feeace5422aa0db, 0x3feeb737b0cdc5e5, 0x3feec49182a3f090, 0x3feed503b23e255d,
    0x3feee89f995ad3ad, 0x3feeff76f2fb5e47, 0x3fef199bdd85529c, 0x3fef3720dcef9069,
    0x3fef5818dcfba487, 0x3fef7c97337b9b5f, 0x3fefa4afa2a490da, 0x3fefd0765b6e4540,
];

const EXP2F_OVFL: f32 = asfloat(0x42B17217);
const EXP2F_UVFL: f32 = asfloat(0xC2CFF1B4);

// ============================================================
// exp (double precision)
// ============================================================

#[inline]
fn exp_top12(x: f64) -> u32 {
    (asuint64(x) >> 52) as u32
}

#[inline]
fn exp_specialcase(tmp: f64, sbits: u64, ki: u64) -> f64 {
    if (ki & 0x80000000) == 0 {
        // k > 0, the exponent of scale might have overflowed by <= 460.
        let sbits = sbits.wrapping_sub(1009u64 << 52);
        let scale = asdouble(sbits);
        let y = EXP_1P1009 * (scale + scale * tmp);
        return eval_as_double(y);
    }
    // k < 0, need special care in the subnormal range.
    let sbits = sbits + (1022u64 << 52);
    let scale = asdouble(sbits);
    let mut y = scale + scale * tmp;
    if y < 1.0 {
        let lo = scale - y + scale * tmp;
        let hi = 1.0 + y;
        let lo = 1.0 - hi + y + lo;
        y = eval_as_double(hi + lo) - 1.0;
        if y == 0.0 {
            y = 0.0; // avoid -0.0 with downward rounding
        }
        // signal underflow exception explicitly
        fp_force_eval(fp_barrier(EXP_1P_NEG1022) * EXP_1P_NEG1022);
    }
    y = EXP_1P_NEG1022 * y;
    eval_as_double(y)
}

#[no_mangle]
pub extern "C" fn exp(x: f64) -> f64 {
    let mut abstop = exp_top12(x) & 0x7ff;

    if predict_false(abstop.wrapping_sub(0x3C9) >= 0x408 - 0x3C9) {
        if abstop.wrapping_sub(0x3C9) >= 0x80000000 {
            // tiny x (including 0): avoid spurious underflow
            return 1.0 + x;
        }
        if abstop >= 0x409 {
            if x == f64::NEG_INFINITY {
                return 0.0;
            }
            if abstop >= 0x7FF {
                return 1.0 + x; // inf or nan
            }
            if asuint64(x) >> 63 != 0 {
                return __math_uflow(0);
            } else {
                return __math_oflow(0);
            }
        }
        // large but finite: force specialcase below
        abstop = 0;
    }

    let z = EXP_INVLN2N * x;
    let kd = eval_as_double(z + EXP_SHIFT);
    let ki = asuint64(kd);
    let kd = kd - EXP_SHIFT;
    let r = x + kd * EXP_NEGLN2HI_N + kd * EXP_NEGLN2LO_N;
    let idx = (2 * (ki % EXP_N)) as usize;
    let top = ki << (52 - EXP_TABLE_BITS);
    let tail = asdouble(EXP_TAB[idx]);
    let sbits = EXP_TAB[idx + 1].wrapping_add(top);
    let r2 = r * r;
    let tmp = tail + r + r2 * (EXP_C2 + r * EXP_C3) + r2 * r2 * (EXP_C4 + r * EXP_C5);

    if predict_false(abstop == 0) {
        return exp_specialcase(tmp, sbits, ki);
    }
    let scale = asdouble(sbits);
    eval_as_double(scale + scale * tmp)
}

// ============================================================
// exp2 (double precision)
// ============================================================

const EXP2_SHIFT: f64 = asdouble(0x42C8000000000000);
const EXP2_1P_NEG54: f64 = asdouble(0x3C90000000000000);
const EXP2_C1: f64 = asdouble(0x3FE62E42FEFA39EF);
const EXP2_C2: f64 = asdouble(0x3FCEBFBDFF82C424);
const EXP2_C3: f64 = asdouble(0x3FAC6B08D70CF4B5);
const EXP2_C4: f64 = asdouble(0x3F83B2ABD24650CC);
const EXP2_C5: f64 = asdouble(0x3F55D7E09B4E3A84);

#[inline]
fn exp2_specialcase(tmp: f64, sbits: u64, ki: u64) -> f64 {
    if (ki & 0x80000000) == 0 {
        // k > 0, the exponent of scale might have overflowed by 1.
        let sbits = sbits.wrapping_sub(1u64 << 52);
        let scale = asdouble(sbits);
        let y = 2.0 * (scale + scale * tmp);
        return eval_as_double(y);
    }
    // k < 0, need special care in the subnormal range.
    let sbits = sbits.wrapping_add(1022u64 << 52);
    let scale = asdouble(sbits);
    let mut y = scale + scale * tmp;
    if y < 1.0 {
        let lo = scale - y + scale * tmp;
        let hi = 1.0 + y;
        let lo = 1.0 - hi + y + lo;
        y = eval_as_double(hi + lo) - 1.0;
        // Avoid -0.0 with downward rounding.
        if y == 0.0 {
            y = 0.0;
        }
        // The underflow exception needs to be signaled explicitly.
        fp_force_eval(fp_barrier(EXP_1P_NEG1022) * EXP_1P_NEG1022);
    }
    y = EXP_1P_NEG1022 * y;
    eval_as_double(y)
}

pub(crate) fn exp2_impl(x: f64) -> f64 {
    let top12 = |v: f64| (asuint64(v) >> 52) as u32;
    let mut abstop = top12(x) & 0x7ff;
    if predict_false(abstop.wrapping_sub(top12(EXP2_1P_NEG54)) >= top12(512.0) - top12(EXP2_1P_NEG54)) {
        if abstop.wrapping_sub(top12(EXP2_1P_NEG54)) >= 0x80000000 {
            // Avoid spurious underflow for tiny x. Note: 0 is common input.
            return 1.0 + x;
        }
        if abstop >= top12(1024.0) {
            if asuint64(x) == asuint64(f64::NEG_INFINITY) {
                return 0.0;
            }
            if abstop >= top12(f64::INFINITY) {
                return 1.0 + x;
            }
            if asuint64(x) >> 63 == 0 {
                return __math_oflow(0);
            } else if asuint64(x) >= asuint64(-1075.0) {
                return __math_uflow(0);
            }
        }
        if 2u64.wrapping_mul(asuint64(x)) > 2u64.wrapping_mul(asuint64(928.0)) {
            // Large x is special cased below.
            abstop = 0;
        }
    }

    // exp2(x) = 2^(k/N) * 2^r, with 2^r in [2^(-1/2N),2^(1/2N)].
    let kd = eval_as_double(x + EXP2_SHIFT);
    let ki = asuint64(kd);
    let kd = kd - EXP2_SHIFT;
    let r = x - kd;
    let idx = (2 * (ki % EXP_N)) as usize;
    let top = ki << (52 - EXP_TABLE_BITS);
    let tail = asdouble(EXP_TAB[idx]);
    let sbits = EXP_TAB[idx + 1].wrapping_add(top);
    let r2 = r * r;
    let tmp = tail
        + r * EXP2_C1
        + r2 * (EXP2_C2 + r * EXP2_C3)
        + r2 * r2 * (EXP2_C4 + r * EXP2_C5);
    if predict_false(abstop == 0) {
        return exp2_specialcase(tmp, sbits, ki);
    }
    let scale = asdouble(sbits);
    eval_as_double(scale + scale * tmp)
}

// ============================================================
// expf (single precision)
// ============================================================

#[inline]
fn expf_top12(x: f32) -> u32 {
    asuint(x) >> 20
}

#[no_mangle]
pub extern "C" fn expf(x: f32) -> f32 {
    let xd = x as f64;
    let abstop = expf_top12(x) & 0x7ff;

    if predict_false(abstop >= 0x42B) {
        if x == f32::NEG_INFINITY {
            return 0.0;
        }
        if abstop >= 0x7F8 {
            return x + x; // inf or nan
        }
        if x > EXP2F_OVFL {
            return __math_oflowf(0);
        }
        if x < EXP2F_UVFL {
            return __math_uflowf(0);
        }
    }

    let z = EXP2F_INVLN2_SCALED * xd;
    let kd = eval_as_double(z + EXP2F_SHIFT);
    let ki = asuint64(kd);
    let kd = kd - EXP2F_SHIFT;
    let r = z - kd;

    let t = EXP2F_TAB[(ki % EXP2F_N) as usize].wrapping_add(ki << (52 - EXP2F_TABLE_BITS));
    let s = asdouble(t);
    let z = EXP2F_POLY_SCALED[0] * r + EXP2F_POLY_SCALED[1];
    let r2 = r * r;
    let mut y = EXP2F_POLY_SCALED[2] * r + 1.0;
    y = z * r2 + y;
    y = y * s;
    y as f32
}