atanBinary32

fn atanBinary32(x: f32) f32

File

Code

fn atanBinary32(x: f32) f32 {
    const atanhi: []const f32 = &.{
        4.6364760399e-01, // atan(0.5)hi 0x3eed6338
        7.8539812565e-01, // atan(1.0)hi 0x3f490fda
        9.8279368877e-01, // atan(1.5)hi 0x3f7b985e
        1.5707962513e+00, // atan(inf)hi 0x3fc90fda
    };
    const atanlo: []const f32 = &.{
        5.0121582440e-09, // atan(0.5)lo 0x31ac3769
        3.7748947079e-08, // atan(1.0)lo 0x33222168
        3.4473217170e-08, // atan(1.5)lo 0x33140fb4
        7.5497894159e-08, // atan(inf)lo 0x33a22168
    };
    const aT: []const f32 = &.{
        3.3333328366e-01,
        -1.9999158382e-01,
        1.4253635705e-01,
        -1.0648017377e-01,
        6.1687607318e-02,
    };

    const hx: u32 = @bitCast(x);
    const ix = hx & 0x7fff_ffff;
    const sign = (hx >> 31) != 0;
    // if |x| >= 2^26
    if (ix >= 0x4c80_0000) {
        if (math.isNan(x)) {
            return x;
        }
        const z = atanhi[3] + 0x1p-120;
        return if (sign) -z else z;
    }
    const x_, const id: ?usize = blk: {
        // |x| < 0.4375
        if (ix < 0x3ee00000) {
            // |x| < 2^(-12)
            if (ix < 0x39800000) {
                if (ix < 0x00800000) {
                    // raise underflow for subnormal x
                    mem.doNotOptimizeAway(x * x);
                }
                return x;
            }
            break :blk .{ x, null };
        } else {
            const x_ = @abs(x);
            // |x| < 1.1875
            if (ix < 0x3f98_0000) {
                // 7/16 <= |x| < 11/16
                if (ix < 0x3f30_0000) {
                    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 (ix < 0x401c_0000) {
                    break :blk .{ (x_ - 1.5) / (1.0 + 1.5 * x_), 2 };
                }
                // 2.4375 <= |x| < 2^26
                else {
                    break :blk .{ -1.0 / x_, 3 };
                }
            }
        }
    };
    // end of argument reduction
    const z = x_ * x_;
    const w = z * z;
    // break sum from i=0 to 10 aT[i]z^(i+1) into odd and even poly
    const s1 = z * (aT[0] + w * (aT[2] + w * aT[4]));
    const s2 = w * (aT[1] + w * aT[3]);
    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);
    }
}