asinBinary128

fn asinBinary128(x: f128) f128 {
    const pio2_hi: f128 = 1.57079632679489661923132169163975140;
    const pio2_lo: f128 = 4.33590506506189051239852201302167613e-35;

    const hx: u128 = @bitCast(x);
    const se: u16 = @truncate(hx >> 112);
    const e = se & 0x7fff;
    const sign = se >> 15 != 0;

    // |x| >= 1 or nan
    if (e >= 0x3fff) {
        // asin(+-1)=+-pi/2 with inexact
        if (x == 1.0 or x == -1.0) {
            return x * pio2_hi + 0x1p-120;
        }
        return 0.0 / (x - x);
    }

    // |x| < 0.5
    if (e < 0x3fff - 1) {
        if (e < 0x3fff - (math.floatMantissaBits(f128) + 2) / 2) {
            // return x with inexact if x!=0
            mem.doNotOptimizeAway(x + 0x1p120);
            return x;
        }
        return x + x * rationalApproxBinary128(x * x);
    }

    // 1 > |x| >= 0.5
    const z = (1.0 - @abs(x)) * 0.5;
    const s = @sqrt(z);
    const r = rationalApproxBinary128(z);

    const top: u16 = @truncate((hx >> 96) & 0x0000_ffff);
    if (top >= 0xee00) {
        const x_local = pio2_hi - (2.0 * (s + s * r) - pio2_lo);
        return if (sign) -x_local else x_local;
    }

    const hs: u128 = @bitCast(s);
    const f: f128 = @bitCast(hs & 0xffff_ffff_ffff_ffff_0000_0000_0000_0000);
    const c = (z - f * f) / (s + f);
    const x_local = 0.5 * pio2_hi - (2.0 * s * r - (pio2_lo - 2.0 * c) - (0.5 * pio2_hi - 2.0 * f));
    return if (sign) -x_local else x_local;
}