sqrtq

pub fn sqrtq(x: f128) callconv(.c) f128

File

Code

pub fn sqrtq(x: f128) callconv(.c) f128 {
    var ix: u128 = @bitCast(x);
    var top = ix >> 112;

    // special case handling.
    if (top -% 0x0001 >= 0x7FFF - 0x0001) {
        @branchHint(.unlikely);
        // x < 0x1p-16382 or inf or nan.
        if (ix & 0x7FFF_FFFF_FFFF_FFFF_FFFF_FFFF_FFFF_FFFF == 0) return x;
        if (ix == 0x7FFF_0000_0000_0000_0000_0000_0000_0000) return x;
        if (ix > 0x7FFF_0000_0000_0000_0000_0000_0000_0000) return math.nan(f128);
        // x is subnormal, normalize it.
        ix = @bitCast(x * 0x1p112);
        top = (ix >> 112) -% 112;
    }

    // argument reduction:
    // x = 4^e m; with integer e, and m in [1, 4)
    // m: fixed point representation [2.126]
    // 2^e is the exponent part of the result.
    const even = (top & 1) != 0;
    const m = if (even) (ix << 14) & 0x7FFF_FFFF_FFFF_FFFF_FFFF_FFFF_FFFF_FFFF else (ix << 15) | 0x8000_0000_0000_0000_0000_0000_0000_0000;
    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 128 bit
    //   m: 2.126 r: 0.128, s: 2.126, d: 2.126, u: 2.126, three: 2.126
    const three: struct { u32, u64, u128 } = .{
        0xC000_0000,
        0xC000_0000_0000_0000,
        0xC000_0000_0000_0000_0000_0000_0000_0000,
    };
    var r: struct { u32, u64, u128 } = undefined;
    var s: struct { u32, u64, u128 } = undefined;
    var d: struct { u32, u64, u128 } = undefined;
    var u: struct { u32, u64, u128 } = undefined;
    const i: usize = @intCast((ix >> 106) & 0x7F);
    r[0] = @intCast(rsqrt_tab[i]);
    r[0] <<= 16;
    // |r sqrt(m) - 1| < 0x1p-8
    s[0] = mul32(@intCast(m >> 96), 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 >> 64), 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 128bit
    r[2] = @intCast(r[1]);
    r[2] <<= 64;
    s[2] = mul128(m, r[2]);
    d[2] = mul128(s[2], r[2]);
    u[2] = three[2] - d[2];
    s[2] = mul128(u[2], s[2]); // repr: 3.125
    // -0x1p-116 < s - sqrt(m) < 0x3.8001p-125
    s[2] = (s[2] - 4) >> 13; // repr: 16.122
    // s < sqrt(m) < s + 1 ULP + tiny

    // compute nearest rounded result:
    // the nearest result to 122 bits is either s or s+0x1p-122,
    // we can decide by comparing (2^122 s + 0.5)^2 to 2^244 m
    const d0 = (m << 98) -% s[2] *% s[2];
    const d1 = s[2] -% d0;
    const d2 = d1 +% s[2] +% 1;
    s[2] += d1 >> 127;
    s[2] &= 0x0000_FFFF_FFFF_FFFF_FFFF_FFFF_FFFF_FFFF;
    s[2] |= top << 112;
    const y: f128 = @bitCast(s[2]);

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

    return y;
}