`nx_f64_div.nx` source

1// nx_f64_div.nx -- IEEE 754 binary64 division (standalone file).
2//
3// Separate from nx_f64.nx mirroring the f32 family split (same-file mul+div
4// codegen quirk, see nx_f32_div.nx header).
5//
6// CORRECTNESS BAR: unlike f32 v1 (which zeroed subnormal operands and lost
7// the guard bit in the sig_a < sig_b case), this is correctly rounded RNE
8// over the full domain: subnormal operands pre-normalize, the restoring
9// division produces 56 quotient bits + exact remainder, guard is explicit
10// in both alignment cases, sticky = remainder nonzero.
11//
12// Restoring division: q = floor(sig_a * 2^55 / sig_b) bit-serial.
13// Working values: r < 2^54 before doubling, < 2^55 after; q < 2^56.
14// Everything fits i64 with headroom (probe-verified arithmetic).
15//
16// license_tier: ORIGINAL
17
18import "nx_syscalls.nx"
19import "nx_tier.nx"
20import "nx_f64.nx"
21
22func nx_f64_div(a: i64, b: i64) -> i64 {
23    let cls_a: nx_int = nx_f64_classify(a)
24    let cls_b: nx_int = nx_f64_classify(b)
25    let sign_out: i64 = nx_f64_sign(a) ^ nx_f64_sign(b)
26
27    if cls_a == NX_F64_CLS_NAN { return NX_F64_NAN_RAW }
28    if cls_b == NX_F64_CLS_NAN { return NX_F64_NAN_RAW }
29
30    if cls_b == NX_F64_CLS_ZERO {
31        if cls_a == NX_F64_CLS_ZERO { return NX_F64_NAN_RAW }
32        return (sign_out << 63) | NX_F64_INF_RAW      // x/0 = inf (incl. inf/0)
33    }
34    if cls_b == NX_F64_CLS_INF {
35        if cls_a == NX_F64_CLS_INF { return NX_F64_NAN_RAW }
36        return sign_out << 63                          // finite/inf = 0
37    }
38    if cls_a == NX_F64_CLS_ZERO { return sign_out << 63 }
39    if cls_a == NX_F64_CLS_INF {
40        return (sign_out << 63) | NX_F64_INF_RAW
41    }
42
43    // Significands pre-normalized to [2^52, 2^53); subnormal exponents go
44    // negative and _f64_round_pack denormalizes the result back.
45    var sig_a: i64 = nx_f64_mant_field(a)
46    var exp_a: i64 = nx_f64_exp_field(a)
47    if exp_a == 0 {
48        exp_a = 1
49        while sig_a < NX_F64_IMPLICIT_1 { sig_a = sig_a << 1; exp_a = exp_a - 1 }
50    } else {
51        sig_a = sig_a | NX_F64_IMPLICIT_1
52    }
53    var sig_b: i64 = nx_f64_mant_field(b)
54    var exp_b: i64 = nx_f64_exp_field(b)
55    if exp_b == 0 {
56        exp_b = 1
57        while sig_b < NX_F64_IMPLICIT_1 { sig_b = sig_b << 1; exp_b = exp_b - 1 }
58    } else {
59        sig_b = sig_b | NX_F64_IMPLICIT_1
60    }
61
62    // q = floor(sig_a * 2^55 / sig_b), r = exact remainder.
63    var q: i64 = 0
64    var r: i64 = sig_a
65    if r >= sig_b { q = 1; r = r - sig_b }
66    var i: i64 = 0
67    while i < 55 {
68        r = r << 1
69        q = q << 1
70        if r >= sig_b { q = q + 1; r = r - sig_b }
71        i = i + 1
72    }
73
74    // sig_a/sig_b in (1/2, 2) so q in (2^54, 2^56).
75    var mant_out: i64 = 0
76    var guard: i64 = 0
77    var sticky: i64 = 0
78    var eo: i64 = 0
79    if q >= (1 << 55) {
80        // Ratio in [1,2): top 53 bits = q >> 3.
81        mant_out = q >> 3
82        guard = (q >> 2) & 1
83        if (q & 3) != 0 { sticky = 1 }
84        if r != 0 { sticky = 1 }
85        eo = exp_a - exp_b + NX_F64_EXP_BIAS
86    } else {
87        // Ratio in (1/2, 1): top 53 bits = q >> 2, exponent one lower.
88        mant_out = q >> 2
89        guard = (q >> 1) & 1
90        if (q & 1) != 0 { sticky = 1 }
91        if r != 0 { sticky = 1 }
92        eo = exp_a - exp_b + NX_F64_EXP_BIAS - 1
93    }
94
95    return _f64_round_pack(sign_out, eo, mant_out, guard, sticky)
96}
nx_f64_div.nx source

`nx_f64_div.nx` source
