asinExtended80

fn asinExtended80(x: f80) f80 {
    const pio2_hi: f80 = 1.57079632679489661926;
    const pio2_lo: f80 = -2.50827880633416601173e-20;

    const hx: u80 = @bitCast(x);
    const se: u16 = @truncate(hx >> 64);
    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(f80) + 1) / 2) {
            // return x with inexact if x!=0
            mem.doNotOptimizeAway(x + 0x1p120);
            return x;
        }
        return x + x * rationalApproxExtended80(x * x);
    }

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

    const m: u64 = @truncate(hx & 0x0000_ffff_ffff_ffff_ffff);
    if ((m >> 56) >= 0xf7) {
        const x_local = pio2_hi - (2.0 * (s + s * r) - pio2_lo);
        return if (sign) -x_local else x_local;
    }

    const hs: u80 = @bitCast(s);
    const f: f80 = @bitCast(hs & 0xffff_ffff_ffff_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;
}