atanBinary128

fn atanBinary128(x: f128) f128 {
    const atanhi: []const f128 = &.{
        4.63647609000806116214256231461214397e-01,
        7.85398163397448309615660845819875699e-01,
        9.82793723247329067985710611014666038e-01,
        1.57079632679489661923132169163975140e+00,
    };
    const atanlo: []const f128 = &.{
        4.89509642257333492668618435220297706e-36,
        2.16795253253094525619926100651083806e-35,
        -2.31288434538183565909319952098066272e-35,
        4.33590506506189051239852201302167613e-35,
    };
    const aT: []const f128 = &.{
        3.33333333333333333333333333333333125e-01,
        -1.99999999999999999999999999999180430e-01,
        1.42857142857142857142857142125269827e-01,
        -1.11111111111111111111110834490810169e-01,
        9.09090909090909090908522355708623681e-02,
        -7.69230769230769230696553844935357021e-02,
        6.66666666666666660390096773046256096e-02,
        -5.88235294117646671706582985209643694e-02,
        5.26315789473666478515847092020327506e-02,
        -4.76190476189855517021024424991436144e-02,
        4.34782608678695085948531993458097026e-02,
        -3.99999999632663469330634215991142368e-02,
        3.70370363987423702891250829918659723e-02,
        -3.44827496515048090726669907612335954e-02,
        3.22579620681420149871973710852268528e-02,
        -3.03020767654269261041647570626778067e-02,
        2.85641979882534783223403715930946138e-02,
        -2.69824879726738568189929461383741323e-02,
        2.54194698498808542954187110873675769e-02,
        -2.35083879708189059926183138130183215e-02,
        2.04832358998165364349957325067131428e-02,
        -1.54489555488544397858507248612362957e-02,
        8.64492360989278761493037861575248038e-03,
        -2.58521121597609872727919154569765469e-03,
    };

    const hx: u128 = @bitCast(x);
    const se: u16 = @truncate(hx >> 112);
    const e = se & 0x7fff;
    const sign = se >> 15 != 0;
    // if |x| is large, atan(x)~=pi/2
    if (e >= 0x3fff + math.floatMantissaBits(f128) + 2) {
        if (math.isNan(x)) {
            return x;
        }
        return if (sign) -atanhi[3] else atanhi[3];
    }
    // Extract the exponent and the first few bits of the mantissa.
    const top: u16 = @truncate((hx >> 96) & 0x0000_ffff);
    const expman = ((@as(u32, @intCast(se)) & 0x7fff) << 8) | (@as(u32, @intCast(top)) >> 8);
    const x_, const id: ?usize = blk: {
        // |x| < 0.4375
        if (expman < ((0x3fff - 2) << 8) + 0xc0) {
            // if |x| is small, atanl(x)~=x
            if (e < 0x3fff - (math.floatMantissaBits(f128) + 2) / 2) {
                // raise underflow if subnormal
                if (e == 0) {
                    mem.doNotOptimizeAway(@as(f32, @floatCast(x)));
                }
                return x;
            }
            break :blk .{ x, null };
        } else {
            const x_ = @abs(x);
            // |x| < 1.1875
            if (expman < (0x3fff << 8) + 0x30) {
                // 7/16 <= |x| < 11/16
                if (expman < ((0x3fff - 1) << 8) + 0x60) {
                    break :blk .{ (2.0 * x_ - 1.0) / (2.0 + x_), 0 };
                }
                // 11/16 <= |x| < 19/16
                else {
                    break :blk .{ (x_ - 1.0) / (x_ + 1.0), 1 };
                }
            } else {
                // |x| < 2.4375
                if (expman < ((0x3fff + 1) << 8) + 0x38) {
                    break :blk .{ (x_ - 1.5) / (1.0 + 1.5 * x_), 2 };
                }
                // 2.4375 <= |x|
                else {
                    break :blk .{ -1.0 / x_, 3 };
                }
            }
        }
    };
    // end of argument reduction
    const z = x_ * x_;
    const w = z * z;
    // break sum aT[i]z^(i+1) into odd and even poly
    const s1 = z * (aT[0] + w * (aT[2] + w * (aT[4] + w * (aT[6] + w * (aT[8] + w * (aT[10] + w * (aT[12] + w * (aT[14] + w * (aT[16] + w * (aT[18] + w * (aT[20] + w * aT[22])))))))))));
    const s2 = w * (aT[1] + w * (aT[3] + w * (aT[5] + w * (aT[7] + w * (aT[9] + w * (aT[11] + w * (aT[13] + w * (aT[15] + w * (aT[17] + w * (aT[19] + w * (aT[21] + w * aT[23])))))))))));
    if (id) |id_| {
        const z_ = atanhi[id_] - ((x_ * (s1 + s2) - atanlo[id_]) - x_);
        return if (sign) -z_ else z_;
    } else {
        return x_ - x_ * (s1 + s2);
    }
}