__sqrth

pub fn __sqrth(x: f16) callconv(.c) f16

File

Code

pub fn __sqrth(x: f16) callconv(.c) f16 {
    var ix: u16 = @bitCast(x);
    var top = ix >> 10;

    // special case handling.
    if (top -% 0x01 >= 0x1F - 0x01) {
        @branchHint(.unlikely);
        // x < 0x1p-14 or inf or nan.
        if (ix & 0x7FFF == 0) return x;
        if (ix == 0x7C00) return x;
        if (ix > 0x7C00) return math.nan(f16);
        // x is subnormal, normalize it.
        ix = @bitCast(x * 0x1p10);
        top = (ix >> 10) -% 10;
    }

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

    // approximate r ~ 1/sqrt(m) and s ~ sqrt(m) when m in [1,4)
    // the fixed point representations are
    //   m: 2.14 r: 0.16, s: 2.14, d: 2.14, u: 2.14, three: 2.14
    const three: u16 = 0xC000;
    const i: usize = @intCast((ix >> 4) & 0x7F);
    const r = rsqrt_tab[i];
    // |r*sqrt(m) - 1| < 0x1p-8
    var s = mul16(m, r);
    // |s/sqrt(m) - 1| < 0x1p-8
    const d = mul16(s, r);
    const u = three - d;
    s = mul16(s, u); // repr: 3.13
    // -0x1.20p-13 < s/sqrt(m) - 1 < 0x7Dp-16
    s = (s - 1) >> 3; // repr: 6.10
    // s < sqrt(m) < s + 0x1.24p-10

    // compute nearest rounded result:
    // the nearest result to 10 bits is either s or s+0x1p-10,
    // we can decide by comparing (2^10 s + 0.5)^2 to 2^20 m.
    const d0 = (m << 6) -% s *% s;
    const d1 = s -% d0;
    const d2 = d1 +% s +% 1;
    s += d1 >> 15;
    s &= 0x03FF;
    s |= top << 10;
    const y: f16 = @bitCast(s);

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

    return y;
}