nx_f64.nx source
↩ module page · 436 lines · 14500 B
1// nx_f64.nx -- IEEE 754 binary64 (double-precision float) bits-up.
2//
3// L5 of the bits-up numeric tower (docs/NISHI_BITS_UP_NUMERIC_TOWER_ROADMAP.md)
4// and ME0 keystone of the math-engine-exceed ladder (DLMF/GAMS-class special
5// functions all sit on binary64 or better). Pure i64 substrate: no libm, no
6// soft-float stubs, no compiler builtins.
7//
8// Bit layout (IEEE 754-2019 binary64):
9// bit 63: sign
10// bits 62..52: exponent (biased 1023; 0 = subnormal/zero; 2047 = inf/NaN)
11// bits 51..0: mantissa (52 bits; implicit leading 1 for normals)
12//
13// Storage: an f64 value IS the i64 holding its bit pattern (full width --
14// unlike f32-in-low-32-bits, raw f64 i64s go negative when sign=1; all field
15// extraction masks after shifting, so arithmetic-shift fill is harmless).
16//
17// The 53x53-bit significand product (106 bits) exceeds i64. We split each
18// significand 26/27 bits so every partial product and partial sum stays under
19// 2^63 (verified at runtime by _f64_probe before this file was authored).
20//
21// CORRECTNESS BAR (above the f32 v1 family): subnormal inputs are
22// pre-normalized, subnormal outputs are denormalized BEFORE rounding (no
23// double-rounding), so add/sub/mul here are correctly rounded round-to-
24// nearest-even across the whole domain. Gate: _f64_gate_authored KATs vs
25// hardware-IEEE oracle bit patterns.
26//
27// File layout note: division lives in nx_f64_div.nx and sqrt in
28// nx_f64_sqrt.nx, mirroring the f32 family split that sidesteps the nxc2
29// same-file mul+div codegen quirk (see nx_f32_div.nx header).
30//
31// genealogy_id: ieee754_2019_binary64 + goldberg_1991_csur
32// + muller_handbook_fp_2018 (referenced, not copied)
33// license_tier: ORIGINAL
34
35import "nx_syscalls.nx"
36import "nx_tier.nx"
37
38const NX_F64_CLS_ZERO: nx_int = 0
39const NX_F64_CLS_NORMAL: nx_int = 1
40const NX_F64_CLS_SUBNORMAL: nx_int = 2
41const NX_F64_CLS_INF: nx_int = 3
42const NX_F64_CLS_NAN: nx_int = 4
43
44const NX_F64_EXP_SHIFT: i64 = 52
45const NX_F64_EXP_BIAS: i64 = 1023
46const NX_F64_EXP_MAXF: i64 = 2047 // exponent field all-ones
47const NX_F64_MANT_MASK: i64 = 0x000FFFFFFFFFFFFF // low 52 bits
48const NX_F64_IMPLICIT_1: i64 = 0x0010000000000000 // 1 << 52
49const NX_F64_SIG_TOP: i64 = 0x0020000000000000 // 1 << 53
50const NX_F64_INF_RAW: i64 = 0x7FF0000000000000
51const NX_F64_NAN_RAW: i64 = 0x7FF8000000000000 // canonical quiet NaN
52const NX_F64_ABS_MASK: i64 = 0x7FFFFFFFFFFFFFFF
53
54// ===== Field extraction ===============================================
55// Raw may be negative (sign bit set); >> sign-extends, masks fix it.
56
57func nx_f64_sign(raw: i64) -> i64 {
58 return (raw >> 63) & 1
59}
60
61func nx_f64_exp_field(raw: i64) -> i64 {
62 return (raw >> NX_F64_EXP_SHIFT) & NX_F64_EXP_MAXF
63}
64
65func nx_f64_mant_field(raw: i64) -> i64 {
66 return raw & NX_F64_MANT_MASK
67}
68
69func nx_f64_classify(raw: i64) -> nx_int {
70 let e: i64 = nx_f64_exp_field(raw)
71 let m: i64 = nx_f64_mant_field(raw)
72 if e == 0 {
73 if m == 0 { return NX_F64_CLS_ZERO }
74 return NX_F64_CLS_SUBNORMAL
75 }
76 if e == NX_F64_EXP_MAXF {
77 if m == 0 { return NX_F64_CLS_INF }
78 return NX_F64_CLS_NAN
79 }
80 return NX_F64_CLS_NORMAL
81}
82
83func nx_f64_is_nan(raw: i64) -> nx_int {
84 if nx_f64_classify(raw) == NX_F64_CLS_NAN { return 1 }
85 return 0
86}
87
88func nx_f64_is_inf(raw: i64) -> nx_int {
89 if nx_f64_classify(raw) == NX_F64_CLS_INF { return 1 }
90 return 0
91}
92
93func nx_f64_is_zero(raw: i64) -> nx_int {
94 if nx_f64_classify(raw) == NX_F64_CLS_ZERO { return 1 }
95 return 0
96}
97
98func nx_f64_neg(raw: i64) -> i64 {
99 return raw ^ (1 << 63)
100}
101
102func nx_f64_abs(raw: i64) -> i64 {
103 return raw & NX_F64_ABS_MASK
104}
105
106// ===== Comparison =====================================================
107
108func nx_f64_eq(a: i64, b: i64) -> nx_int {
109 if nx_f64_is_nan(a) == 1 { return 0 }
110 if nx_f64_is_nan(b) == 1 { return 0 }
111 if a == b { return 1 }
112 // +0 == -0
113 if (a & NX_F64_ABS_MASK) == 0 {
114 if (b & NX_F64_ABS_MASK) == 0 { return 1 }
115 }
116 return 0
117}
118
119// a < b per IEEE totalOrder-free semantics (NaN compares false).
120func nx_f64_lt(a: i64, b: i64) -> nx_int {
121 if nx_f64_is_nan(a) == 1 { return 0 }
122 if nx_f64_is_nan(b) == 1 { return 0 }
123 let za: i64 = nx_f64_is_zero(a)
124 let zb: i64 = nx_f64_is_zero(b)
125 if za == 1 {
126 if zb == 1 { return 0 } // ±0 < ±0 is false
127 }
128 let sa: i64 = nx_f64_sign(a)
129 let sb: i64 = nx_f64_sign(b)
130 if sa != sb {
131 // a negative, b positive -> true, UNLESS both zero (handled above).
132 if sa == 1 {
133 if zb == 1 {
134 if za == 1 { return 0 }
135 }
136 return 1
137 }
138 return 0
139 }
140 // Same sign: compare magnitude bits (exp+mant concatenated = monotone).
141 let ma: i64 = a & NX_F64_ABS_MASK
142 let mb: i64 = b & NX_F64_ABS_MASK
143 if sa == 0 {
144 if ma < mb { return 1 }
145 return 0
146 }
147 if ma > mb { return 1 }
148 return 0
149}
150
151func nx_f64_gt(a: i64, b: i64) -> nx_int {
152 return nx_f64_lt(b, a)
153}
154
155// ===== Round + pack (shared by mul / div / sqrt / cvt / add) ==========
156//
157// Inputs: sign 0/1; mant NORMALIZED in [2^52, 2^53) (or smaller only when
158// eo == 1, the already-denormal case from add); guard / sticky in {0,1};
159// eo = biased exponent such that value = mant * 2^(eo - 1075).
160// eo may be <= 0: the value is denormalized HERE, feeding shifted-out bits
161// into guard/sticky BEFORE round-to-nearest-even -- no double rounding.
162
163func _f64_round_pack(sign: i64, eo_in: i64, mant_in: i64, guard_in: i64, sticky_in: i64) -> i64 {
164 var eo: i64 = eo_in
165 var mant: i64 = mant_in
166 var guard: i64 = guard_in
167 var sticky: i64 = sticky_in
168
169 // Denormalize while the exponent is below the minimum normal field (1).
170 if eo < 1 {
171 let k: i64 = 1 - eo
172 if k > 54 {
173 // Everything shifts out; only stickiness remains -> rounds to 0.
174 return sign << 63
175 }
176 var i: i64 = 0
177 while i < k {
178 if guard == 1 { sticky = 1 }
179 guard = mant & 1
180 mant = mant >> 1
181 i = i + 1
182 }
183 eo = 1
184 }
185
186 // Round-to-nearest-even.
187 var round_up: i64 = 0
188 if guard == 1 {
189 if sticky == 1 { round_up = 1 }
190 if round_up == 0 {
191 if (mant & 1) == 1 { round_up = 1 }
192 }
193 }
194 if round_up == 1 {
195 mant = mant + 1
196 if mant >= NX_F64_SIG_TOP {
197 mant = mant >> 1
198 eo = eo + 1
199 }
200 }
201
202 // Overflow -> inf.
203 if eo >= NX_F64_EXP_MAXF {
204 if mant >= NX_F64_IMPLICIT_1 {
205 return (sign << 63) | NX_F64_INF_RAW
206 }
207 }
208
209 // Subnormal (or zero) output: exponent field 0.
210 if mant < NX_F64_IMPLICIT_1 {
211 return (sign << 63) | mant
212 }
213
214 // Normal output.
215 return (sign << 63) | (eo << NX_F64_EXP_SHIFT) | (mant & NX_F64_MANT_MASK)
216}
217
218// ===== Multiplication =================================================
219//
220// 53x53-bit significand product = 106 bits, held as P_hi (bits 105..53)
221// and P_lo (bits 52..0) via 26/27-bit operand splits. Both operands are
222// pre-normalized to [2^52, 2^53) (subnormals shift left, exponent goes
223// negative, _f64_round_pack denormalizes back) so the product is always
224// in [2^104, 2^106) and normalization is a single top-bit test.
225
226func nx_f64_mul(a: i64, b: i64) -> i64 {
227 let cls_a: nx_int = nx_f64_classify(a)
228 let cls_b: nx_int = nx_f64_classify(b)
229 let sign_out: i64 = nx_f64_sign(a) ^ nx_f64_sign(b)
230
231 if cls_a == NX_F64_CLS_NAN { return NX_F64_NAN_RAW }
232 if cls_b == NX_F64_CLS_NAN { return NX_F64_NAN_RAW }
233
234 if cls_a == NX_F64_CLS_INF {
235 if cls_b == NX_F64_CLS_ZERO { return NX_F64_NAN_RAW }
236 return (sign_out << 63) | NX_F64_INF_RAW
237 }
238 if cls_b == NX_F64_CLS_INF {
239 if cls_a == NX_F64_CLS_ZERO { return NX_F64_NAN_RAW }
240 return (sign_out << 63) | NX_F64_INF_RAW
241 }
242
243 if cls_a == NX_F64_CLS_ZERO { return sign_out << 63 }
244 if cls_b == NX_F64_CLS_ZERO { return sign_out << 63 }
245
246 // Build significands; pre-normalize subnormals.
247 var sig_a: i64 = nx_f64_mant_field(a)
248 var exp_a: i64 = nx_f64_exp_field(a)
249 if exp_a == 0 {
250 exp_a = 1
251 while sig_a < NX_F64_IMPLICIT_1 { sig_a = sig_a << 1; exp_a = exp_a - 1 }
252 } else {
253 sig_a = sig_a | NX_F64_IMPLICIT_1
254 }
255 var sig_b: i64 = nx_f64_mant_field(b)
256 var exp_b: i64 = nx_f64_exp_field(b)
257 if exp_b == 0 {
258 exp_b = 1
259 while sig_b < NX_F64_IMPLICIT_1 { sig_b = sig_b << 1; exp_b = exp_b - 1 }
260 } else {
261 sig_b = sig_b | NX_F64_IMPLICIT_1
262 }
263
264 // 106-bit product via 26/27 split (every partial < 2^54).
265 let a0: i64 = sig_a & 0x7FFFFFF // low 27 bits
266 let a1: i64 = sig_a >> 27 // high 26 bits
267 let b0: i64 = sig_b & 0x7FFFFFF
268 let b1: i64 = sig_b >> 27
269
270 let hh: i64 = a1 * b1 // contributes at bit 54
271 let mid: i64 = a1 * b0 + a0 * b1 // contributes at bit 27; < 2^54
272 let ll: i64 = a0 * b0 // contributes at bit 0; < 2^54
273
274 // Assemble into P_hi (bits 105..53) / P_lo (bits 52..0), all sums < 2^63.
275 let m53: i64 = (1 << 53) - 1
276 let lo_lo: i64 = ll & m53
277 let lo_hi: i64 = ll >> 53 // 0 or 1
278 let mid_low: i64 = (mid & ((1 << 26) - 1)) << 27 // bits 27..52
279 let mid_high: i64 = mid >> 26 // at bit 53
280 let plo_raw: i64 = lo_lo + mid_low // < 2^54
281 let carry: i64 = plo_raw >> 53
282 let p_lo: i64 = plo_raw & m53
283 let p_hi: i64 = (hh << 1) + mid_high + lo_hi + carry // bits 105..53
284
285 var mant_out: i64 = 0
286 var guard: i64 = 0
287 var sticky: i64 = 0
288 var eo: i64 = 0
289 if p_hi >= NX_F64_IMPLICIT_1 {
290 // Product top bit at 105: mant = prod >> 53 = p_hi exactly.
291 mant_out = p_hi
292 guard = (p_lo >> 52) & 1
293 if (p_lo & ((1 << 52) - 1)) != 0 { sticky = 1 }
294 eo = exp_a + exp_b - NX_F64_EXP_BIAS + 1
295 } else {
296 // Top bit at 104: mant = prod >> 52.
297 mant_out = (p_hi << 1) | (p_lo >> 52)
298 guard = (p_lo >> 51) & 1
299 if (p_lo & ((1 << 51) - 1)) != 0 { sticky = 1 }
300 eo = exp_a + exp_b - NX_F64_EXP_BIAS
301 }
302
303 return _f64_round_pack(sign_out, eo, mant_out, guard, sticky)
304}
305
306// ===== Addition =======================================================
307//
308// Ported from the gated nx_f32_add structure, widened: 53-bit significands
309// + 3 GRS bits = 56-bit working width; carry tops at bit 56; all in i64.
310// Sticky-jam into bit 0 before the magnitude add/sub (Berkeley-SoftFloat-
311// style; with G+R+S this yields correct round-to-nearest-even).
312
313func nx_f64_add(a: i64, b: i64) -> i64 {
314 let cls_a: nx_int = nx_f64_classify(a)
315 let cls_b: nx_int = nx_f64_classify(b)
316
317 if cls_a == NX_F64_CLS_NAN { return NX_F64_NAN_RAW }
318 if cls_b == NX_F64_CLS_NAN { return NX_F64_NAN_RAW }
319
320 let sign_a: i64 = nx_f64_sign(a)
321 let sign_b: i64 = nx_f64_sign(b)
322
323 if cls_a == NX_F64_CLS_INF {
324 if cls_b == NX_F64_CLS_INF {
325 if sign_a == sign_b { return a }
326 return NX_F64_NAN_RAW
327 }
328 return a
329 }
330 if cls_b == NX_F64_CLS_INF { return b }
331
332 if cls_a == NX_F64_CLS_ZERO {
333 if cls_b == NX_F64_CLS_ZERO {
334 if sign_a == sign_b { return a }
335 return 0 // +0 + -0 = +0 (RNE)
336 }
337 return b
338 }
339 if cls_b == NX_F64_CLS_ZERO { return a }
340
341 // Significands; subnormals carry effective exponent 1, no implicit bit.
342 var sig_a_raw: i64 = nx_f64_mant_field(a)
343 var exp_a_raw: i64 = nx_f64_exp_field(a)
344 if exp_a_raw == 0 {
345 exp_a_raw = 1
346 } else {
347 sig_a_raw = sig_a_raw | NX_F64_IMPLICIT_1
348 }
349 var sig_b_raw: i64 = nx_f64_mant_field(b)
350 var exp_b_raw: i64 = nx_f64_exp_field(b)
351 if exp_b_raw == 0 {
352 exp_b_raw = 1
353 } else {
354 sig_b_raw = sig_b_raw | NX_F64_IMPLICIT_1
355 }
356
357 // Order so a holds the larger-or-equal exponent.
358 var s_a: i64 = sign_a
359 var s_b: i64 = sign_b
360 var e_a: i64 = exp_a_raw
361 var e_b: i64 = exp_b_raw
362 var m_a: i64 = sig_a_raw
363 var m_b: i64 = sig_b_raw
364 if e_b > e_a {
365 let tmp_s: i64 = s_a; s_a = s_b; s_b = tmp_s
366 let tmp_e: i64 = e_a; e_a = e_b; e_b = tmp_e
367 let tmp_m: i64 = m_a; m_a = m_b; m_b = tmp_m
368 }
369
370 var m_a_shifted: i64 = m_a << 3
371 var m_b_shifted: i64 = m_b << 3
372
373 let exp_diff: i64 = e_a - e_b
374 var sticky_b: i64 = 0
375 if exp_diff > 0 {
376 if exp_diff >= 57 {
377 if m_b_shifted != 0 { sticky_b = 1 }
378 m_b_shifted = 0
379 } else {
380 let lost_mask: i64 = (1 << exp_diff) - 1
381 if (m_b_shifted & lost_mask) != 0 { sticky_b = 1 }
382 m_b_shifted = m_b_shifted >> exp_diff
383 }
384 }
385 if sticky_b == 1 { m_b_shifted = m_b_shifted | 1 }
386
387 var result_sig: i64 = 0
388 var result_sign: i64 = s_a
389
390 if s_a == s_b {
391 result_sig = m_a_shifted + m_b_shifted
392 } else {
393 if m_a_shifted >= m_b_shifted {
394 result_sig = m_a_shifted - m_b_shifted
395 } else {
396 result_sig = m_b_shifted - m_a_shifted
397 result_sign = s_b
398 }
399 if result_sig == 0 { return 0 } // exact cancellation
400 }
401
402 var exp_out: i64 = e_a
403
404 // Carry: top moved to bit 56.
405 if result_sig >= (1 << 56) {
406 let lost: i64 = result_sig & 1
407 result_sig = result_sig >> 1
408 if lost == 1 { result_sig = result_sig | 1 }
409 exp_out = exp_out + 1
410 }
411
412 // Cancellation renormalize: top bit back to 55 (or hit the subnormal floor).
413 var renorm_done: nx_int = 0
414 while renorm_done == 0 {
415 if result_sig >= (1 << 55) { renorm_done = 1 }
416 if renorm_done == 0 {
417 if exp_out <= 1 { renorm_done = 1 }
418 }
419 if renorm_done == 0 {
420 result_sig = result_sig << 1
421 exp_out = exp_out - 1
422 }
423 }
424
425 // GRS: bit 2 = guard; bits 1..0 fold into sticky.
426 let guard: i64 = (result_sig >> 2) & 1
427 var sticky: i64 = 0
428 if (result_sig & 3) != 0 { sticky = 1 }
429 let mant_out: i64 = result_sig >> 3
430
431 return _f64_round_pack(result_sign, exp_out, mant_out, guard, sticky)
432}
433
434func nx_f64_sub(a: i64, b: i64) -> i64 {
435 return nx_f64_add(a, nx_f64_neg(b))
436}