sqrt

pub fn sqrt(x: f64) callconv(.c) f64

File

Code

pub fn sqrt(x: f64) callconv(.c) f64 {
    var ix: u64 = @bitCast(x);
    var top = ix >> 52;

    // special case handling.
    if (top -% 0x001 >= 0x7FF - 0x001) {
        @branchHint(.unlikely);
        // x < 0x1p-1022 or inf or nan.
        if (ix & 0x7FFF_FFFF_FFFF_FFFF == 0) return x;
        if (ix == 0x7FF0_0000_0000_0000) return x;
        if (ix > 0x7FF0_0000_0000_0000) return if (compiler_rt.want_float_exceptions) (x - x) / 0.0 else math.nan(f64);
        // x is subnormal, normalize it.
        ix = @bitCast(x * 0x1p52);
        top = (ix >> 52) -% 52;
    }

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

    // approximate r ~ 1/sqrt(m) and s ~ sqrt(m) when m in [1,4)
    //
    // initial estimate:
    // 7bit table lookup (1bit exponent and 6bit significand).
    //
    // iterative approximation:
    // using 2 goldschmidt iterations with 32bit int arithmetics
    // and a final iteration with 64bit int arithmetics.
    //
    // details:
    //
    // the relative error (e = r0 sqrt(m)-1) of a linear estimate
    // (r0 = a m + b) is |e| < 0.085955 ~ 0x1.6p-4 at best,
    // a table lookup is faster and needs one less iteration
    // 6 bit lookup table (128b) gives |e| < 0x1.f9p-8
    // 7 bit lookup table (256b) gives |e| < 0x1.fdp-9
    // for single and double prec 6bit is enough but for quad
    // prec 7bit is needed (or modified iterations). to avoid
    // one more iteration >=13bit table would be needed (16k).
    //
    // a newton-raphson iteration for r is
    //   w = r*r
    //   u = 3 - m*w
    //   r = r*u/2
    // can use a goldschmidt iteration for s at the end or
    //   s = m*r
    //
    // first goldschmidt iteration is
    //   s = m*r
    //   u = 3 - s*r
    //   r = r*u/2
    //   s = s*u/2
    // next goldschmidt iteration is
    //   u = 3 - s*r
    //   r = r*u/2
    //   s = s*u/2
    // and at the end r is not computed only s.
    //
    // they use the same amount of operations and converge at the
    // same quadratic rate, i.e. if
    //   r1 sqrt(m) - 1 = e, then
    //   r2 sqrt(m) - 1 = -3/2 e^2 - 1/2 e^3
    // the advantage of goldschmidt is that the mul for s and r
    // are independent (computed in parallel), however it is not
    // "self synchronizing": it only uses the input m in the
    // first iteration so rounding errors accumulate. at the end
    // or when switching to larger precision arithmetics rounding
    // errors dominate so the first iteration should be used.
    //
    // 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
    const three: struct { u32, u64 } = .{
        0xC000_0000,
        0xC000_0000_0000_0000,
    };
    var r: struct { u32, u64 } = undefined;
    var s: struct { u32, u64 } = undefined;
    var d: struct { u32, u64 } = undefined;
    var u: struct { u32, u64 } = undefined;
    const i: usize = @intCast((ix >> 46) & 0x7F);
    r[0] = @intCast(rsqrt_tab[i]);
    r[0] <<= 16;
    // |r sqrt(m) - 1| < 0x1.fdp-9
    s[0] = mul32(@intCast(m >> 32), r[0]);
    // |s/sqrt(m) - 1| < 0x1.fdp-9
    d[0] = mul32(s[0], r[0]);
    u[0] = three[0] - d[0];
    r[0] = mul32(r[0], u[0]) << 1;
    // |r sqrt(m) - 1| < 0x1.7bp-16
    s[0] = mul32(s[0], u[0]) << 1;
    // |s/sqrt(m) - 1| < 0x1.7bp-16
    d[0] = mul32(s[0], r[0]);
    u[0] = three[0] - d[0];
    r[0] = mul32(r[0], u[0]) << 1;
    // |r sqrt(m) - 1| < 0x1.3704p-29 (measured worst-case)
    r[1] = @intCast(r[0]);
    r[1] <<= 32;
    s[1] = mul64(m, r[1]);
    d[1] = mul64(s[1], r[1]);
    u[1] = three[1] - d[1];
    s[1] = mul64(s[1], u[1]); // repr: 3.61
    // -0x1p-57 < s - sqrt(m) < 0x1.8001p-61
    s[1] = (s[1] - 2) >> 9; // repr: 12.52
    // -0x1.09p-52 < s - sqrt(m) < -0x1.fffcp-63

    // s < sqrt(m) < s + 0x1.09p-52
    // compute nearest rounded result:
    // the nearest result to 52 bits is either s or s+0x1p-52,
    // we can decide by comparing (2^52 s + 0.5)^2 to 2^104 m.
    const d0 = (m << 42) -% s[1] *% s[1];
    const d1 = s[1] -% d0;
    const d2 = d1 +% s[1] +% 1;
    s[1] += d1 >> 63;
    s[1] &= 0x000F_FFFF_FFFF_FFFF;
    s[1] |= top << 52;
    const y: f64 = @bitCast(s[1]);

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

    return y;
}