`nx_f64_sqrt.nx` source

1// nx_f64_sqrt.nx -- IEEE 754 binary64 square root (standalone file).
2//
3// Digit-pair (2-bits-per-iteration) restoring square root on the radicand
4// sig * 2^54, fed MSB-first as 54 bit-pairs (27 pairs from the 54-bit
5// adjusted significand, then 27 zero pairs).  Produces a 54-bit root =
6// 53-bit mantissa + explicit guard; sticky = remainder nonzero.  Correctly
7// rounded RNE for all finite positive inputs including subnormals.
8//
9// Bounds: root < 2^54, trial = (root<<2)|1 < 2^56, rem <= 2*root + pair
10// so rem << 2 < 2^57 -- all comfortably inside i64.
11//
12// IEEE specials: sqrt(+0) = +0, sqrt(-0) = -0, sqrt(+inf) = +inf,
13// sqrt(negative) = NaN, sqrt(NaN) = NaN.
14//
15// license_tier: ORIGINAL
16
17import "nx_syscalls.nx"
18import "nx_tier.nx"
19import "nx_f64.nx"
20const K_MAGIC_1049: i64 = 1049
21
22func nx_f64_sqrt(a: i64) -> i64 {
23    let cls_a: nx_int = nx_f64_classify(a)
24
25    if cls_a == NX_F64_CLS_NAN { return NX_F64_NAN_RAW }
26    if cls_a == NX_F64_CLS_ZERO { return a }           // preserves -0
27    if nx_f64_sign(a) == 1 { return NX_F64_NAN_RAW }   // negative -> NaN
28    if cls_a == NX_F64_CLS_INF { return a }            // +inf
29
30    // Pre-normalize significand to [2^52, 2^53).
31    var sig: i64 = nx_f64_mant_field(a)
32    var ef: i64 = nx_f64_exp_field(a)
33    if ef == 0 {
34        ef = 1
35        while sig < NX_F64_IMPLICIT_1 { sig = sig << 1; ef = ef - 1 }
36    } else {
37        sig = sig | NX_F64_IMPLICIT_1
38    }
39
40    // value = sig * 2^(E - 52), E unbiased.  Make E even so the exponent
41    // halves cleanly; the odd case doubles sig (now < 2^54).
42    var e_unb: i64 = ef - NX_F64_EXP_BIAS
43    // Parity via bit 0 (two's complement keeps bit 0 = parity for negatives;
44    // E - 2*(E/2) would be wrong there under truncating division).
45    let parity: i64 = e_unb & 1
46    if parity != 0 {
47        sig = sig << 1
48        e_unb = e_unb - 1
49    }
50    // Now E even; sqrt(value) = sqrt(sig) * 2^((E-52)/2)... computed as
51    // root = floor(sqrt(sig * 2^54)) (in [2^53, 2^54)), giving
52    // sqrt(value) = (root + frac) * 2^(F - 27) with F = (E - 52) / 2.
53    // E - 52 is even and may be negative; arithmetic >> 1 halves correctly.
54    let f_half: i64 = (e_unb - 52) >> 1
55
56    // Digit-pair restoring sqrt: 27 pairs from sig (54 bits), 27 zero pairs.
57    var rem: i64 = 0
58    var root: i64 = 0
59    var idx: i64 = 26
60    while idx >= 0 {
61        let pair: i64 = (sig >> (idx * 2)) & 3
62        rem = (rem << 2) | pair
63        let trial: i64 = (root << 2) | 1
64        if rem >= trial {
65            rem = rem - trial
66            root = (root << 1) | 1
67        } else {
68            root = root << 1
69        }
70        idx = idx - 1
71    }
72    var j: i64 = 0
73    while j < 27 {
74        rem = rem << 2
75        let trial2: i64 = (root << 2) | 1
76        if rem >= trial2 {
77            rem = rem - trial2
78            root = (root << 1) | 1
79        } else {
80            root = root << 1
81        }
82        j = j + 1
83    }
84
85    // root has 54 bits: mantissa 53 + guard 1.  value = root * 2^(f_half-27).
86    let mant_out: i64 = root >> 1
87    let guard: i64 = root & 1
88    var sticky: i64 = 0
89    if rem != 0 { sticky = 1 }
90    // mant * 2^(eo-1075) = root/2 * 2^(eo-1075)  =>  eo = f_half + 1049.
91    let eo: i64 = f_half + K_MAGIC_1049
92
93    return _f64_round_pack(0, eo, mant_out, guard, sticky)
94}
nx_f64_sqrt.nx source

`nx_f64_sqrt.nx` source
