__sqrtx

pub fn __sqrtx(x: f80) callconv(.c) f80

File

Code

pub fn __sqrtx(x: f80) callconv(.c) f80 {
    var ix: u80 = @bitCast(x);
    var top = ix >> 64;

    // special case handling.
    if (top -% 0x0001 >= 0x7FFF - 0x0001) {
        @branchHint(.unlikely);
        // x < 0x1p-16382 or inf or nan.
        if (ix & 0x7FFF_FFFF_FFFF_FFFF_FFFF == 0) return x;
        if (ix == 0x7FFF_8000_0000_0000_0000) return x;
        if (ix > 0x7FFF_8000_0000_0000_0000) return if (compiler_rt.want_float_exceptions) (x - x) / 0.0 else math.nan(f80);
        // x is subnormal, normalize it.
        ix = @bitCast(x * 0x1p63);
        top = (ix >> 64) -% 63;
    }

    // argument reduction:
    // x = 4^e m; with integer e, and m in [1, 4)
    // m: fixed point representation [2.78]
    // 2^e is the exponent part of the result.
    const even = (top & 1) != 0;
    const m = if (even) (ix << 15) & 0x7FFF_FFFF_FFFF_FFFF_FFFF else ix << 16;
    top = (top +% 0x3FFF) >> 1;

    // approximate r ~ 1/sqrt(m) and s ~ sqrt(m) when m in [1,4)
    // the fixed point representations are
    //   m: 2.30 r: 0.32, s: 2.30, d: 2.30, u: 2.30, three: 2.30
    // and after switching to 64 bit
    //   m: 2.62 r: 0.64, s: 2.62, d: 2.62, u: 2.62, three: 2.62
    // and after switching to 80 bit
    //   m: 2.78 r: 0.80, s: 2.78, d: 2.78, u: 2.78, three: 2.78
    const three: struct { u32, u64, u80 } = .{
        0xC000_0000,
        0xC000_0000_0000_0000,
        0xC000_0000_0000_0000_0000,
    };
    var r: struct { u32, u64, u80 } = undefined;
    var s: struct { u32, u64, u80 } = undefined;
    var d: struct { u32, u64, u80 } = undefined;
    var u: struct { u32, u64, u80 } = undefined;
    var i: usize = @intCast((ix >> 57) & 0x3F);
    if (even) i += 64;
    r[0] = @intCast(rsqrt_tab[i]);
    r[0] <<= 16;
    // |r sqrt(m) - 1| < 0x1p-8
    s[0] = mul32(@intCast(m >> 48), r[0]);
    d[0] = mul32(s[0], r[0]);
    u[0] = three[0] - d[0];
    r[0] = mul32(u[0], r[0]) << 1;
    // |r sqrt(m) - 1| < 0x1.7bp-16, switch to 64bit
    r[1] = @intCast(r[0]);
    r[1] <<= 32;
    s[1] = mul64(@intCast(m >> 16), r[1]);
    d[1] = mul64(s[1], r[1]);
    u[1] = three[1] - d[1];
    r[1] = mul64(u[1], r[1]) << 1;
    // |r sqrt(m) - 1| < 0x1.a5p-31
    s[1] = mul64(u[1], s[1]) << 1;
    d[1] = mul64(s[1], r[1]);
    u[1] = three[1] - d[1];
    r[1] = mul64(u[1], r[1]) << 1;
    // |r sqrt(m) - 1| < 0x1.c001p-59, switch to 80bit
    r[2] = @intCast(r[1]);
    r[2] <<= 16;
    s[2] = mul80(m, r[2]);
    d[2] = mul80(s[2], r[2]);
    u[2] = three[2] - d[2];
    s[2] = mul80(u[2], s[2]); // repr: 3.77
    s[2] = (s[2] - 4) >> 14; // repr: 17.63
    // s < sqrt(m) < s + 1 ULP + tiny

    // compute nearest rounded result:
    // the nearest result to 63 bits is either s or s+0x1p-63,
    // we can decide by comparing (2^63 s + 0.5)^2 to 2^126 m
    const d0 = (m << 48) -% mul80_tail(s[2], s[2]);
    const d1 = s[2] -% d0;
    const d2 = d1 +% s[2] +% 1;
    s[2] += d1 >> 79;
    s[2] &= 0x0000_7FFF_FFFF_FFFF_FFFF;
    s[2] |= 0x0000_8000_0000_0000_0000;
    s[2] |= top << 64;
    const y: f80 = @bitCast(s[2]);

    // handle rounding modes and inexact exception:
    // only (s+1)^2 == 2^48 m case is exact otherwise
    // add a tiny value to cause the fenv effects.
    if (d2 != 0) {
        @branchHint(.likely);
        var tiny: u80 = 0x0001_8000_0000_0000_0000;
        tiny |= (d1 ^ d2) & 0x8000_0000_0000_0000_0000;
        const t: f80 = @bitCast(tiny);
        return y + t;
    }

    return y;
}