nx_fixq30_lib.nx source
↩ module page · 335 lines · 14332 B
1// nx_fixq30_lib.nx -- Q30 FIXED-POINT MATH: exact multiply and divide, sqrt, CORDIC sin/cos/atan2, exp, ln,
2// pow, cbrt. (2026-08-24, the phototwin PT1 keystone.) license_tier: ORIGINAL No hw writes (Rule 26).
3//
4// WHY Q30 AND WHY A NEW LIB. The estate's math primitives are Q10/Q14 (nx_trig, nx_exp, nx_icbrt, nx_isqrt_q10),
5// ~1e-3 to 6e-5 precision, built for attention kernels and rendering. CIEDE2000 is published to four decimals
6// (Sharma, Wu, Dalal 2005, Table I) and a hue error of 1e-4 rad on a chroma of 100 is already 0.01 dE00, so a
7// colour ruler needs ~1e-8 arithmetic. Q30 (unit 9.3e-10) on i64 leaves 33 integer bits: every CIELAB quantity
8// (|L|,|a|,|b| <= 200, C^2 <= 4e4) fits with room for the CORDIC gain. Those libs also carry a main() and so
9// cannot be imported (double-main), which is why fq_isqrt is a second copy of the digit-by-digit integer sqrt
10// -- declared here, not hidden.
11//
12// NO FLOATING POINT, NO TYPED-IN CONSTANTS. Every number below is a power of two, a formula integer (Machin's
13// 5 and 239, the atanh argument 3 for ln 2), or DERIVED at runtime by series (pi, ln 2, the CORDIC atan table,
14// the CORDIC gain). Nothing is copied from a calculator, so nothing can be mistyped; the gate KATs pi and e
15// against their published digits as the ORACLE, which is what a KAT is for.
16//
17// CONTEXT: fq_ctx() computes the derived constants ONCE and hands back an i64 slot table; every trig/exp/ln
18// call takes it. Slots: [FQ_C_PI] pi, [FQ_C_LN2] ln 2, [FQ_C_GAIN] CORDIC gain K, [FQ_C_ATAN0 + i] atan(2^-i).
19//
20// DOMAINS (each proven at the function): fq_mul is EXACT for |a|,|b| < 2^46 (65536 real); fq_div for |b| < 2^47
21// and |a/b| < 2^33 real; fq_exp saturates to FQ_SAT above x = 22.18 and underflows to 0 below x = -42.9.
22import "nx_syscalls.nx"
23import "nx_vecmath.nx"
24
25const FQ_SHIFT: i64 = 30
26const FQ_ONE: i64 = 1073741824
27const FQ_HALF: i64 = 536870912
28const FQ_LO_MASK: i64 = 1073741823
29const FQ_SPLIT: i64 = 15 // fq_div refines the remainder in two 15-bit steps (2*15 = FQ_SHIFT)
30const FQ_DIV_MAX_DEN: i64 = 140737488355328 // 2^47: the remainder step (r << 15, r < denominator) must stay below 2^63
31const FQ_ISQRT_TARGET_BITS: i64 = 62 // fq_sqrt pre-scales its argument to this many bits before the integer sqrt
32const FQ_ISQRT_HALF_SHIFT: i64 = 15 // sqrt(x * 2^30) = sqrt(x) * 2^15
33const FQ_SERIES_MAX_TERMS: i64 = 64 // a series that has not converged in 64 terms is being misused
34const FQ_MACHIN_A: i64 = 5 // Machin (1706): pi/4 = 4 atan(1/5) - atan(1/239)
35const FQ_MACHIN_B: i64 = 239
36const FQ_MACHIN_KA: i64 = 4
37const FQ_ATANH_LN2_ARG: i64 = 3 // ln 2 = 2 atanh(1/3)
38const FQ_CORDIC_ITERS: i64 = 30 // = FQ_SHIFT: atan(2^-i) is below one Q30 unit past i = 30
39const FQ_EXP_MAX_POW2: i64 = 32 // exp saturates once its 2^n factor would exceed 2^31 (x > 22.18)
40const FQ_SAT: i64 = 4611686018427387904 // 2^62: returned instead of overflowing
41const FQ_NEG_SAT: i64 = 0 - 4611686018427387904
42const FQ_DEG_HALF_TURN: i64 = 180
43const FQ_MICRO: i64 = 1000000
44const FQ_DECIMAL: i64 = 10
45const FQ_MICRO_DIGITS: i64 = 6
46const FQ_ASCII_ZERO: i64 = 48
47const FQ_ASCII_NINE: i64 = 57
48const FQ_ASCII_MINUS: i64 = 45
49const FQ_ASCII_DOT: i64 = 46
50const FQ_TWO: i64 = 2
51const FQ_THREE: i64 = 3
52const FQ_CBRT_NEWTON_STEPS: i64 = 2 // after the exp(ln/3) seed two Newton steps reach the Q30 floor
53const FQ_C_PI: i64 = 0
54const FQ_C_LN2: i64 = 1
55const FQ_C_GAIN: i64 = 2
56const FQ_C_ATAN0: i64 = 8
57const FQ_CTX_SLOTS: i64 = 40 // FQ_C_ATAN0 + FQ_CORDIC_ITERS, rounded up
58const FQ_SLOT_BYTES: i64 = 8
59
60func fq_abs(v: i64) -> i64 { if v < 0 { return 0 - v } return v }
61func fq_from_int(n: i64) -> i64 { return n * FQ_ONE }
62func fq_to_int(x: i64) -> i64 { if x >= 0 { return x >> FQ_SHIFT } return 0 - ((0 - x) >> FQ_SHIFT) }
63// right shift that truncates toward zero for either sign (CORDIC needs a sign-safe halving, not an arithmetic shift)
64func fq_shr(v: i64, n: i64) -> i64 {
65 if n <= 0 { return v }
66 if n >= FQ_ISQRT_TARGET_BITS { return 0 }
67 if v >= 0 { return v >> n }
68 return 0 - ((0 - v) >> n)
69}
70func fq_bitlen(v: i64) -> i64 { var n: i64 = 0; var x: i64 = v; while x > 0 { x = x >> 1; n = n + 1 } return n }
71// micro-units (1e-6) <-> Q30. from_micro: m << 30 fits while |m| < 2^33 (8.6e3 real); to_micro: x * 1e6 fits while |x| < 2^43.
72func fq_from_micro(m: i64) -> i64 {
73 if m >= 0 { return (m << FQ_SHIFT) / FQ_MICRO }
74 return 0 - (((0 - m) << FQ_SHIFT) / FQ_MICRO)
75}
76func fq_to_micro(x: i64) -> i64 {
77 if x >= 0 { return (x * FQ_MICRO) >> FQ_SHIFT }
78 return 0 - (((0 - x) * FQ_MICRO) >> FQ_SHIFT)
79}
80
81// EXACT Q30 multiply without i128 for |a|,|b| < 2^46: split each operand at the Q30 boundary.
82// xh*yh < 2^32 (shifted back up: < 2^62), xh*yl and xl*yh < 2^46, xl*yl < 2^60 -- every partial fits i64.
83func fq_mul(a: i64, b: i64) -> i64 {
84 var neg: i64 = 0
85 var x: i64 = a
86 var y: i64 = b
87 if x < 0 { x = 0 - x; neg = 1 - neg }
88 if y < 0 { y = 0 - y; neg = 1 - neg }
89 let xh: i64 = x >> FQ_SHIFT
90 let xl: i64 = x & FQ_LO_MASK
91 let yh: i64 = y >> FQ_SHIFT
92 let yl: i64 = y & FQ_LO_MASK
93 let r: i64 = ((xh * yh) << FQ_SHIFT) + (xh * yl) + (xl * yh) + ((xl * yl) >> FQ_SHIFT)
94 if neg == 1 { return 0 - r }
95 return r
96}
97// Q30 divide: integer quotient, then the remainder refined in two 15-bit steps (r < |b| < 2^48 keeps r << 15 in i64).
98func fq_div(a: i64, b: i64) -> i64 {
99 if b == 0 { if a >= 0 { return FQ_SAT } return FQ_NEG_SAT }
100 var neg: i64 = 0
101 var x: i64 = a
102 var y: i64 = b
103 if x < 0 { x = 0 - x; neg = 1 - neg }
104 if y < 0 { y = 0 - y; neg = 1 - neg }
105 // DOMAIN GUARD (measured 2026-08-24): a 2^54 denominator overflowed the remainder step silently and a colour LUT came
106 // out garbage while every KAT stayed green. Out of domain, both operands are halved together until the denominator
107 // fits -- the ratio is preserved to the bits that survive, and nothing wraps.
108 while y >= FQ_DIV_MAX_DEN { x = x >> 1; y = y >> 1 }
109 let q: i64 = x / y
110 var r: i64 = x - q * y
111 let t1: i64 = (r << FQ_SPLIT) / y
112 r = (r << FQ_SPLIT) - t1 * y
113 let t2: i64 = (r << FQ_SPLIT) / y
114 let res: i64 = (q << FQ_SHIFT) + (t1 << FQ_SPLIT) + t2
115 if neg == 1 { return 0 - res }
116 return res
117}
118// floor(sqrt(n)) for n >= 0, digit-by-digit in base 4 (the nx_isqrt algorithm; that file carries a main and
119// cannot be imported, so this copy is declared rather than hidden).
120func fq_isqrt(n: i64) -> i64 { return vm_isqrt(n) }
121// sqrt in Q30: scale x up to ~62 bits so the integer sqrt keeps ~31 significant bits, then rescale.
122// sqrt(x_real * 2^(30+2k)) = sqrt(x_real) * 2^(15+k); the Q30 answer is that times 2^(15-k).
123func fq_sqrt(x: i64) -> i64 {
124 if x <= 0 { return 0 }
125 let n: i64 = fq_bitlen(x)
126 var k: i64 = 0
127 if n < FQ_ISQRT_TARGET_BITS { k = (FQ_ISQRT_TARGET_BITS - n) / FQ_TWO }
128 let y: i64 = x << (k * FQ_TWO)
129 let r: i64 = fq_isqrt(y)
130 if k <= FQ_ISQRT_HALF_SHIFT { return r << (FQ_ISQRT_HALF_SHIFT - k) }
131 return r >> (k - FQ_ISQRT_HALF_SHIFT)
132}
133func fq_hypot(a: i64, b: i64) -> i64 { return fq_sqrt(fq_mul(a, a) + fq_mul(b, b)) }
134
135// atan(x) = x - x^3/3 + x^5/5 - ... for 0 <= x <= 1/2 (terms shrink by at least 1/4 per step)
136func fq_atan_small(x: i64) -> i64 {
137 let x2: i64 = fq_mul(x, x)
138 var term: i64 = x
139 var sum: i64 = x
140 var k: i64 = 1
141 var sign: i64 = 0 - 1
142 var i: i64 = 0
143 while i < FQ_SERIES_MAX_TERMS {
144 term = fq_mul(term, x2)
145 k = k + FQ_TWO
146 let t: i64 = term / k
147 if t == 0 { i = FQ_SERIES_MAX_TERMS } else { sum = sum + sign * t; sign = 0 - sign; i = i + 1 }
148 }
149 return sum
150}
151// atanh(t) = t + t^3/3 + t^5/5 + ... for 0 <= t <= 1/3 (terms shrink by at least 1/9 per step)
152func fq_atanh_small(t: i64) -> i64 {
153 let t2: i64 = fq_mul(t, t)
154 var term: i64 = t
155 var sum: i64 = t
156 var k: i64 = 1
157 var i: i64 = 0
158 while i < FQ_SERIES_MAX_TERMS {
159 term = fq_mul(term, t2)
160 k = k + FQ_TWO
161 let u: i64 = term / k
162 if u == 0 { i = FQ_SERIES_MAX_TERMS } else { sum = sum + u; i = i + 1 }
163 }
164 return sum
165}
166
167// THE CONTEXT: pi (Machin), ln 2 (atanh), the CORDIC table atan(2^-i), the CORDIC gain K = 1/sqrt(prod(1+2^-2i)).
168func fq_ctx() -> *i64 {
169 let c: *i64 = sys_mmap(FQ_CTX_SLOTS * FQ_SLOT_BYTES) as *i64
170 let inv_a: i64 = FQ_ONE / FQ_MACHIN_A
171 let inv_b: i64 = FQ_ONE / FQ_MACHIN_B
172 let quarter: i64 = FQ_MACHIN_KA * fq_atan_small(inv_a) - fq_atan_small(inv_b)
173 c[FQ_C_PI] = quarter * FQ_MACHIN_KA
174 c[FQ_C_LN2] = FQ_TWO * fq_atanh_small(FQ_ONE / FQ_ATANH_LN2_ARG)
175 var i: i64 = 0
176 var p: i64 = FQ_ONE
177 while i < FQ_CORDIC_ITERS {
178 if i == 0 { c[FQ_C_ATAN0] = quarter } else { c[FQ_C_ATAN0 + i] = fq_atan_small(FQ_ONE >> i) }
179 p = fq_mul(p, FQ_ONE + (FQ_ONE >> (i * FQ_TWO)))
180 i = i + 1
181 }
182 c[FQ_C_GAIN] = fq_div(FQ_ONE, fq_sqrt(p))
183 return c
184}
185
186// sin and cos of an angle in Q30 radians (any magnitude): out2[0] = sin, out2[1] = cos. CORDIC rotation mode.
187func fq_sincos(ctx: *i64, ang: i64, out2: *i64) -> i64 {
188 let pi: i64 = ctx[FQ_C_PI]
189 let two_pi: i64 = pi * FQ_TWO
190 let half_pi: i64 = pi / FQ_TWO
191 var z: i64 = ang
192 if z > pi { let n1: i64 = (z + pi) / two_pi; z = z - n1 * two_pi }
193 if z < (0 - pi) { let n2: i64 = (pi - z) / two_pi; z = z + n2 * two_pi }
194 var flip: i64 = 0
195 if z > half_pi { z = z - pi; flip = 1 }
196 if z < (0 - half_pi) { z = z + pi; flip = 1 }
197 var x: i64 = ctx[FQ_C_GAIN]
198 var y: i64 = 0
199 var i: i64 = 0
200 while i < FQ_CORDIC_ITERS {
201 let at: i64 = ctx[FQ_C_ATAN0 + i]
202 if z >= 0 {
203 let xn: i64 = x - fq_shr(y, i)
204 y = y + fq_shr(x, i)
205 x = xn
206 z = z - at
207 } else {
208 let xm: i64 = x + fq_shr(y, i)
209 y = y - fq_shr(x, i)
210 x = xm
211 z = z + at
212 }
213 i = i + 1
214 }
215 if flip == 1 { x = 0 - x; y = 0 - y }
216 out2[0] = y
217 out2[1] = x
218 return 0
219}
220func fq_sin(ctx: *i64, ang: i64) -> i64 { let o: *i64 = sys_mmap(FQ_SLOT_BYTES * FQ_TWO) as *i64; fq_sincos(ctx, ang, o); return o[0] }
221func fq_cos(ctx: *i64, ang: i64) -> i64 { let o: *i64 = sys_mmap(FQ_SLOT_BYTES * FQ_TWO) as *i64; fq_sincos(ctx, ang, o); return o[1] }
222
223// atan2(y, x) in Q30 radians, range (-pi, pi]. CORDIC vectoring mode; x < 0 handled by point reflection.
224func fq_atan2(ctx: *i64, y: i64, x: i64) -> i64 {
225 let pi: i64 = ctx[FQ_C_PI]
226 if x == 0 {
227 if y > 0 { return pi / FQ_TWO }
228 if y < 0 { return 0 - (pi / FQ_TWO) }
229 return 0
230 }
231 var xx: i64 = x
232 var yy: i64 = y
233 var addz: i64 = 0
234 if x < 0 {
235 xx = 0 - x
236 yy = 0 - y
237 if y >= 0 { addz = pi } else { addz = 0 - pi }
238 }
239 var z: i64 = 0
240 var i: i64 = 0
241 while i < FQ_CORDIC_ITERS {
242 let at: i64 = ctx[FQ_C_ATAN0 + i]
243 if yy >= 0 {
244 let xn: i64 = xx + fq_shr(yy, i)
245 yy = yy - fq_shr(xx, i)
246 xx = xn
247 z = z + at
248 } else {
249 let xm: i64 = xx - fq_shr(yy, i)
250 yy = yy + fq_shr(xx, i)
251 xx = xm
252 z = z - at
253 }
254 i = i + 1
255 }
256 return z + addz
257}
258func fq_deg2rad(ctx: *i64, deg: i64) -> i64 { return fq_div(fq_mul(deg, ctx[FQ_C_PI]), fq_from_int(FQ_DEG_HALF_TURN)) }
259func fq_rad2deg(ctx: *i64, rad: i64) -> i64 { return fq_div(fq_mul(rad, fq_from_int(FQ_DEG_HALF_TURN)), ctx[FQ_C_PI]) }
260
261// exp(x): x = n ln2 + r with r in [0, ln2); exp(r) by Taylor to the Q30 floor; then scale by 2^n.
262func fq_exp(ctx: *i64, x: i64) -> i64 {
263 let ln2: i64 = ctx[FQ_C_LN2]
264 var n: i64 = x / ln2
265 if x < 0 { if n * ln2 != x { n = n - 1 } }
266 let r: i64 = x - n * ln2
267 var term: i64 = FQ_ONE
268 var sum: i64 = FQ_ONE
269 var k: i64 = 1
270 while k < FQ_SERIES_MAX_TERMS {
271 term = fq_mul(term, r) / k
272 if term == 0 { k = FQ_SERIES_MAX_TERMS } else { sum = sum + term; k = k + 1 }
273 }
274 if n >= 0 {
275 if n >= FQ_EXP_MAX_POW2 { return FQ_SAT }
276 return sum << n
277 }
278 let m: i64 = 0 - n
279 if m >= FQ_ISQRT_TARGET_BITS { return 0 }
280 return sum >> m
281}
282// ln(x) for x > 0: x = m 2^e with m in [1, 2); ln x = e ln2 + 2 atanh((m-1)/(m+1)), the argument <= 1/3.
283func fq_ln(ctx: *i64, x: i64) -> i64 {
284 if x <= 0 { return FQ_NEG_SAT }
285 var m: i64 = x
286 var e: i64 = 0
287 while m >= FQ_ONE * FQ_TWO { m = m >> 1; e = e + 1 }
288 while m < FQ_ONE { m = m << 1; e = e - 1 }
289 let t: i64 = fq_div(m - FQ_ONE, m + FQ_ONE)
290 return e * ctx[FQ_C_LN2] + FQ_TWO * fq_atanh_small(t)
291}
292// x^p for x > 0 (x <= 0 returns 0: the caller decides what a non-positive base means)
293func fq_pow(ctx: *i64, x: i64, p: i64) -> i64 {
294 if x <= 0 { return 0 }
295 return fq_exp(ctx, fq_mul(p, fq_ln(ctx, x)))
296}
297// cube root, sign-symmetric: seed exp(ln(x)/3), then Newton y = (2y + x/y^2)/3
298func fq_cbrt(ctx: *i64, x: i64) -> i64 {
299 if x == 0 { return 0 }
300 var neg: i64 = 0
301 var v: i64 = x
302 if v < 0 { v = 0 - v; neg = 1 }
303 var y: i64 = fq_exp(ctx, fq_ln(ctx, v) / FQ_THREE)
304 var i: i64 = 0
305 while i < FQ_CBRT_NEWTON_STEPS {
306 let y2: i64 = fq_mul(y, y)
307 if y2 == 0 { i = FQ_CBRT_NEWTON_STEPS } else { y = (FQ_TWO * y + fq_div(v, y2)) / FQ_THREE; i = i + 1 }
308 }
309 if neg == 1 { return 0 - y }
310 return y
311}
312
313// parse a signed decimal (up to 6 fractional digits) from buf[off..off+len) into MICRO units; returns micro
314func fq_parse_micro(buf: *u8, off: i64, len: i64) -> i64 {
315 var i: i64 = 0
316 var neg: i64 = 0
317 var ip: i64 = 0
318 var fp: i64 = 0
319 var fd: i64 = 0
320 var infrac: i64 = 0
321 while i < len {
322 let c: i64 = buf[off + i] as i64
323 if c == FQ_ASCII_MINUS { neg = 1 }
324 if c == FQ_ASCII_DOT { infrac = 1 }
325 if c >= FQ_ASCII_ZERO { if c <= FQ_ASCII_NINE {
326 if infrac == 0 { ip = ip * FQ_DECIMAL + (c - FQ_ASCII_ZERO) }
327 else { if fd < FQ_MICRO_DIGITS { fp = fp * FQ_DECIMAL + (c - FQ_ASCII_ZERO); fd = fd + 1 } }
328 } }
329 i = i + 1
330 }
331 while fd < FQ_MICRO_DIGITS { fp = fp * FQ_DECIMAL; fd = fd + 1 }
332 let v: i64 = ip * FQ_MICRO + fp
333 if neg == 1 { return 0 - v }
334 return v
335}