fx.nx source
↩ module page · 358 lines · 13910 B
1// fx.nx -- deterministic Q16.16 fixed-point math.
2//
3// Why fixed-point and not IEEE-754 float:
4// Rollback netcode (GGPO 2006, Skullgirls, Killer Instinct) requires
5// bit-identical simulation across machines. IEEE-754 + Math.sin
6// cannot deliver this -- different JS engines (V8, SpiderMonkey,
7// JavaScriptCore) ship different sin/cos implementations, and even
8// addition order can diverge under JIT inlining. Fixed-point
9// integer math is the only way to guarantee "same input → same bits
10// on every machine."
11//
12// Representation: Q16.16
13// - 64-bit signed integer: high 16 bits integer part, low 16 bits
14// fractional part. Negative values use two's complement.
15// - Range: ±32767.99998... blocks. Precision: 1/65536 ≈ 15 μm at
16// 1 m = 1 block. Adequate for voxel games, physics, gameplay
17// sim up to ~32 km arenas.
18//
19// Invariants (enforced, not hoped):
20// FX1 No IEEE-754 anywhere. NishiLang has no float type today, so
21// this is trivial -- but if floats ever land, fx.nx does not
22// use them.
23// FX2 fx_mul and fx_div round toward zero (arithmetic shift, not
24// round-to-nearest). Deterministic and simple; callers needing
25// round-to-nearest do `(a + (b >> 1)) / b` explicitly.
26// FX3 sin/cos via CORDIC algorithm -- 16 iterations, ~Q16.16
27// precision, no Math.sin dependency. Converges through a
28// precomputed atan table. Every implementation (NishiVM,
29// native C VM, future Nishi silicon) produces bit-identical
30// results by construction.
31// FX4 Overflow in fx_mul is possible if both operands exceed 16-bit
32// integer range. Callers are responsible for clamping; we do
33// NOT implement saturation arithmetic by default (silent
34// saturation masks bugs; explicit clamp is better).
35//
36// References:
37// - Volder 1959, "The CORDIC Trigonometric Computing Technique"
38// - Kota-Kuroda-Shimamura 1989, "A High-Speed Fixed-Point Multiplier
39// Using Carry-Propagation-Free Adders"
40// - GGPO 2006 whitepaper (rollback determinism rationale)
41
42import "syscalls.nx"
43
44// ---- scale knobs ---------------------------------------------------
45
46const FX_SHIFT: i64 = 16
47const FX_ONE: i64 = 65536 // 2^16
48const FX_HALF: i64 = 32768 // 0.5 in Q16.16
49const FX_FRAC_MASK: i64 = 0xFFFF // low 16 bits
50const FX_TWO_PI: i64 = 411775 // round(2π * 65536); 1 rev in Q16.16
51const FX_PI: i64 = 205887 // round(π * 65536)
52const FX_HALF_PI: i64 = 102944 // round(π/2 * 65536)
53
54// CORDIC scaling factor K_n for n=16 iterations: ≈ 0.6072529350088814
55// Stored in Q16.16 = round(0.6072529 * 65536) = 39797.
56const FX_CORDIC_K: i64 = 39797
57
58// ---- basic arithmetic ---------------------------------------------
59
60// Addition / subtraction are plain integer ops -- the same-format
61// operands, same-format result invariant means no scaling needed.
62func fx_add(a: i64, b: i64) -> i64 { return a + b }
63func fx_sub(a: i64, b: i64) -> i64 { return a - b }
64func fx_neg(a: i64) -> i64 { return 0 - a }
65
66// Q16.16 * Q16.16 = Q32.32; we shift right by 16 to get Q16.16 back.
67// Careful: intermediate product can overflow i64 if a and b both
68// exceed ~ Q16.0 range. FX4 responsibility.
69func fx_mul(a: i64, b: i64) -> i64 {
70 return (a * b) >> FX_SHIFT
71}
72
73// Q16.16 / Q16.16 = Q16.16. Shift LHS up by 16 before dividing so
74// the result carries the correct fractional bits. Does integer
75// truncation toward zero (FX2).
76func fx_div(a: i64, b: i64) -> i64 {
77 return (a << FX_SHIFT) / b
78}
79
80// ---- conversions ---------------------------------------------------
81
82func fx_from_int(n: i64) -> i64 { return n << FX_SHIFT }
83
84// Truncate Q16.16 to integer (toward zero).
85func fx_to_int(a: i64) -> i64 { return a >> FX_SHIFT }
86
87// Construct Q16.16 from (integer, numerator, denominator): useful for
88// exact rationals like "1/3 of a block".
89func fx_from_frac(num: i64, den: i64) -> i64 {
90 return (num << FX_SHIFT) / den
91}
92
93// ---- binary logarithm (Q16.16) -----------------------------------
94//
95// fx_log2(x): base-2 logarithm of a Q16.16 value, result in Q16.16.
96//
97// Why this exists: rank-discounted retrieval gain (DCG/nDCG, Jarvelin
98// & Kekalainen 2002) needs 1/log2(rank+1). IEEE-754 log2 cannot be
99// bit-reproducible across machines (FX1); this fixed-point form is.
100//
101// Algorithm (classic, no tables, no float):
102// 1. Characteristic: shift x into [1.0, 2.0); the shift count is the
103// integer part of log2 (negative if x < 1.0).
104// 2. Mantissa: FX_SHIFT iterations of "square the value; if it
105// crossed 2.0, emit a fractional bit and halve it" -- each bit is
106// worth half the previous (1/2, 1/4, ...). Deterministic by
107// construction; same bits on every target.
108//
109// Domain: x must be > 0. log2 of <= 0 is undefined; we return 0 and
110// leave it to the caller's guard (defensive at the boundary, Rule #12).
111// Exact on powers of two: fx_log2(2^k) == k * FX_ONE.
112func fx_log2(x: i64) -> i64 {
113 if x <= 0 { return 0 }
114 var v: i64 = x
115 var int_log: i64 = 0
116 while v >= (FX_ONE << 1) { v = v >> 1; int_log = int_log + 1 }
117 while v < FX_ONE { v = v << 1; int_log = int_log - 1 }
118 // v now in [1.0, 2.0). Seed result with the integer part.
119 var result: i64 = int_log << FX_SHIFT
120 var b: i64 = FX_HALF // weight of the next fractional bit
121 var i: i64 = 0
122 while i < FX_SHIFT {
123 v = fx_mul(v, v)
124 if v >= (FX_ONE << 1) {
125 v = v >> 1
126 result = result + b
127 }
128 b = b >> 1
129 i = i + 1
130 }
131 return result
132}
133
134// Convenience: log2 of a plain integer n (n >= 1), result in Q16.16.
135// fx_log2_int(2) == FX_ONE; the workhorse for DCG rank discounts.
136func fx_log2_int(n: i64) -> i64 {
137 return fx_log2(fx_from_int(n))
138}
139
140// ---- CORDIC sin/cos -----------------------------------------------
141//
142// Rotation-mode CORDIC. Input: angle in Q16.16 radians (call
143// fx_normalize_angle first if you can't guarantee |angle| <= π/2).
144// Output: writes sin to *sin_out, cos to *cos_out, both Q16.16.
145//
146// Precomputed atan(2^-i) in Q16.16 for i = 0..15. Each entry is
147// round(atan(2^-i) * 65536). These are constant across all targets
148// by spec (not by local IEEE-754 evaluation).
149func fx_cordic_atan(i: i64) -> i64 {
150 if i == 0 { return 51472 } // atan(1) * 2^16
151 if i == 1 { return 30386 } // atan(1/2)
152 if i == 2 { return 16055 } // atan(1/4)
153 if i == 3 { return 8150 } // atan(1/8)
154 if i == 4 { return 4091 } // atan(1/16)
155 if i == 5 { return 2047 } // atan(1/32)
156 if i == 6 { return 1024 } // atan(1/64)
157 if i == 7 { return 512 } // atan(1/128)
158 if i == 8 { return 256 } // atan(1/256)
159 if i == 9 { return 128 }
160 if i == 10 { return 64 }
161 if i == 11 { return 32 }
162 if i == 12 { return 16 }
163 if i == 13 { return 8 }
164 if i == 14 { return 4 }
165 return 2 // atan(2^-15)
166}
167
168// Normalise angle into [-π, π] (Q16.16). Strips off full-rotation
169// multiples of 2π so CORDIC convergence stays stable.
170func fx_normalize_angle(a: i64) -> i64 {
171 var r: i64 = a
172 while r > FX_PI { r = r - FX_TWO_PI }
173 while r < (0 - FX_PI) { r = r + FX_TWO_PI }
174 return r
175}
176
177// Compute sin(theta) and cos(theta) for theta in Q16.16. Writes both
178// outputs to caller-supplied slots so callers can take both without
179// a second pass.
180//
181// CORDIC range is [-π/2, π/2] in the basic form. To cover the full
182// [-π, π] we reduce via quadrant rules:
183// θ ∈ [π/2, π] -> θ' = π - θ; sin = sin(θ'); cos = -cos(θ')
184// θ ∈ [-π, -π/2] -> θ' = -π - θ; sin = -sin(θ'); cos = -cos(θ')
185func fx_sin_cos(theta_raw: i64, sin_out: *i64, cos_out: *i64) -> i64 {
186 let theta: i64 = fx_normalize_angle(theta_raw)
187 var t: i64 = theta
188 var flip_cos: i64 = 0
189 var flip_sin: i64 = 0
190
191 if t > FX_HALF_PI {
192 t = FX_PI - t
193 flip_cos = 1
194 }
195 if t < (0 - FX_HALF_PI) {
196 // theta in [-pi, -pi/2]: reflect through -pi to theta' = -pi - theta in [-pi/2, 0].
197 // sin(theta) = sin(-pi-theta') = sin(theta') -> NO sin flip (the original flip_sin=1 was a latent
198 // third-quadrant sign bug, never caught because haversine uses sin^2; found via the f32 FFT vs full-
199 // range DFT gate 2026-06-14). cos(theta) = cos(-pi-theta') = -cos(theta') -> flip cos.
200 t = (0 - FX_PI) - t
201 flip_cos = 1
202 }
203
204 // Core CORDIC loop: start at (K, 0), iteratively converge on
205 // angle t.
206 var x: i64 = FX_CORDIC_K
207 var y: i64 = 0
208 var z: i64 = t
209 var i: i64 = 0
210 while i < 16 {
211 let atan_i: i64 = fx_cordic_atan(i)
212 let dx: i64 = y >> i
213 let dy: i64 = x >> i
214 if z >= 0 {
215 x = x - dx
216 y = y + dy
217 z = z - atan_i
218 } else {
219 x = x + dx
220 y = y - dy
221 z = z + atan_i
222 }
223 i = i + 1
224 }
225
226 if flip_sin == 1 { y = 0 - y }
227 if flip_cos == 1 { x = 0 - x }
228
229 *sin_out = y
230 *cos_out = x
231 return 0
232}
233
234// Convenience wrappers for callers that only want one component.
235func fx_sin(theta: i64) -> i64 {
236 let s_raw: *u8 = sys_mmap(16)
237 let c_raw: *u8 = sys_mmap(16)
238 let s: *i64 = s_raw as *i64
239 let c: *i64 = c_raw as *i64
240 fx_sin_cos(theta, s, c)
241 return *s
242}
243
244func fx_cos(theta: i64) -> i64 {
245 let s_raw: *u8 = sys_mmap(16)
246 let c_raw: *u8 = sys_mmap(16)
247 let s: *i64 = s_raw as *i64
248 let c: *i64 = c_raw as *i64
249 fx_sin_cos(theta, s, c)
250 return *c
251}
252
253// ---- game-sim convenience: u8 yaw angle --------------------------
254//
255// Many game protocols (including ours) transmit yaw as a single byte:
256// 256 discrete angles equally spaced around the circle. This helper
257// converts a yaw byte to Q16.16 radians using the exact formula
258// angle = yaw * 2π / 256 (Q16.16)
259// implemented as (yaw * FX_TWO_PI) / 256, where FX_TWO_PI is the
260// same integer on every target -- no target-dependent rounding.
261func fx_yaw_u8_to_angle(yaw: i64) -> i64 {
262 let y: i64 = yaw & 0xFF
263 return (y * FX_TWO_PI) / 256
264}
265
266// ---- base-2 exponential (Q16.16) ---------------------------------
267//
268// fx_exp2(y) = 2^y, the exact inverse of fx_log2. Dual bit-by-bit
269// algorithm: the integer part is a shift; the fractional part f in [0,1)
270// is built as a PRODUCT of precomputed 2^(2^-k) factors, one per set
271// fractional bit (the mirror of fx_log2's square-and-emit-a-bit loop).
272// Exact on integers: fx_exp2(k*FX_ONE) == 2^k * FX_ONE. Round-trips:
273// fx_exp2(fx_log2(x)) ~= x. Unblocks MA-predicted log-domain gains and
274// the neural postfilter (both need a fixed-point antilog).
275//
276// 2^(2^-k) in Q16.16, k=1..16. Constant on every target by spec
277// (computed offline, NOT by local IEEE-754) -- matches fx_cordic_atan.
278func fx_exp2_factor(k: i64) -> i64 {
279 if k==1 { return 92682 }
280 if k==2 { return 77936 }
281 if k==3 { return 71468 }
282 if k==4 { return 68438 }
283 if k==5 { return 66971 }
284 if k==6 { return 66250 }
285 if k==7 { return 65892 }
286 if k==8 { return 65714 }
287 if k==9 { return 65625 }
288 if k==10 { return 65580 }
289 if k==11 { return 65558 }
290 if k==12 { return 65547 }
291 if k==13 { return 65542 }
292 if k==14 { return 65539 }
293 if k==15 { return 65537 }
294 if k==16 { return 65537 }
295 return FX_ONE
296}
297
298// 2^f for f in Q16.16 within [0, 1). Result in [FX_ONE, 2*FX_ONE).
299func fx_exp2_frac(f: i64) -> i64 {
300 var r: i64 = FX_ONE
301 var k: i64 = 1
302 while k <= FX_SHIFT {
303 if (f & (1 << (FX_SHIFT - k))) != 0 { r = fx_mul(r, fx_exp2_factor(k)) }
304 k = k + 1
305 }
306 return r
307}
308
309// 2^y for Q16.16 y. Saturates above ~2^46 (sentinel) and to 0 below ~2^-31.
310func fx_exp2(y: i64) -> i64 {
311 let ip: i64 = y >> FX_SHIFT // floor via arithmetic shift
312 let fp: i64 = y & FX_FRAC_MASK // fractional part in [0, FX_ONE)
313 let m: i64 = fx_exp2_frac(fp) // 2^fp in [1.0, 2.0)
314 if ip >= 0 {
315 if ip > 46 { return 1 << 62 }
316 return m << ip
317 }
318 let s: i64 = 0 - ip
319 if s > 31 { return 0 }
320 return m >> s
321}
322
323// Change-of-base helpers: 10^x = 2^(x*log2 10); e^x = 2^(x*log2 e);
324// b^e = 2^(e*log2 b). Constants in Q16.16.
325func fx_exp10(x: i64) -> i64 { return fx_exp2(fx_mul(x, 217706)) }
326func fx_expe(x: i64) -> i64 { return fx_exp2(fx_mul(x, 94548)) }
327func fx_pow(base: i64, ex: i64) -> i64 {
328 if base <= 0 { return 0 }
329 return fx_exp2(fx_mul(ex, fx_log2(base)))
330}
331
332// ---- compile-only smoke ------------------------------------------
333//
334// Exercises each primitive. Real validation (bit-exact comparison
335// with reference values computed offline) comes in fx_test.nx once
336// an execution harness is available on this host.
337func main() -> i64 {
338 // Trig smoke: sin(0) should be 0, cos(0) should be FX_ONE.
339 let zero: i64 = 0
340 let s0: i64 = fx_sin(zero)
341 let c0: i64 = fx_cos(zero)
342
343 // sin(π/2) ≈ 1.0 → close to FX_ONE.
344 let s1: i64 = fx_sin(FX_HALF_PI)
345
346 // Multiply: 2.5 * 4.0 = 10.0 in Q16.16.
347 let a: i64 = FX_ONE * 2 + FX_HALF // 2.5
348 let b: i64 = FX_ONE * 4 // 4.0
349 let prod: i64 = fx_mul(a, b) // expect 10 << 16
350
351 // Divide: 1.0 / 3.0 ≈ 21845 in Q16.16.
352 let third: i64 = fx_div(FX_ONE, FX_ONE * 3)
353
354 // Yaw conversion: yaw=64 (quarter-turn) → angle ≈ π/2.
355 let q: i64 = fx_yaw_u8_to_angle(64)
356
357 return s0 + c0 + s1 + prod + third + q
358}