atanExtended80

fn atanExtended80(x: f80) f80

File

Code

fn atanExtended80(x: f80) f80 {
    const atanhi: []const f80 = &.{
        4.63647609000806116202e-01,
        7.85398163397448309628e-01,
        9.82793723247329067960e-01,
        1.57079632679489661926e+00,
    };
    const atanlo: []const f80 = &.{
        1.18469937025062860669e-20,
        -1.25413940316708300586e-20,
        2.55232234165405176172e-20,
        -2.50827880633416601173e-20,
    };
    const aT: []const f80 = &.{
        3.33333333333333333017e-01,
        -1.99999999999999632011e-01,
        1.42857142857046531280e-01,
        -1.11111111100562372733e-01,
        9.09090902935647302252e-02,
        -7.69230552476207730353e-02,
        6.66661718042406260546e-02,
        -5.88158892835030888692e-02,
        5.25499891539726639379e-02,
        -4.70119845393155721494e-02,
        4.03539201366454414072e-02,
        -2.91303858419364158725e-02,
        1.24822046299269234080e-02,
    };

    const hx: u80 = @bitCast(x);
    const se: u16 = @truncate(hx >> 64);
    const e = se & 0x7fff;
    const sign = se >> 15 != 0;
    // if |x| is large, atan(x)~=pi/2
    if (e >= 0x3fff + math.floatMantissaBits(f80) + 1) {
        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 m: u64 = @truncate(hx & 0x0000_ffff_ffff_ffff_ffff);
    const expman = ((@as(u32, @intCast(se)) & 0x7fff) << 8) | (@as(u32, @truncate(m >> 55)) & 0xff);
    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(f80) + 1) / 2) {
                // raise underflow if subnormal
                if (e == 0) {
                    std.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]))))));
    const s2 = w * (aT[1] + w * (aT[3] + w * (aT[5] + w * (aT[7] + w * (aT[9] + w * aT[11])))));
    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);
    }
}