nx_dct8.nx source
↩ module page · 195 lines · 6576 B
1// nx_dct8.nx -- 8x8 integer DCT-II in pure NishiLang. The universal
2// frequency-domain transform that underlies every modern video codec
3// (JPEG, MPEG-1/2/4, H.261/263/264/265, AV1's DCT mode). We re-derive
4// from the public math (Ahmed, Natarajan, Rao 1974; Chen, Smith,
5// Fralick 1977; Loeffler, Ligtenberg, Moschytz 1989) -- not from any
6// patented codec's reference code.
7//
8// Math: 1D DCT-II of size 8:
9// X_k = (1/2) c_k * sum_{n=0..7} x_n * cos((2n+1) k pi / 16)
10// c_0 = 1/sqrt(2); c_k = 1 (k > 0)
11//
12// 2D 8x8 DCT: apply 1D DCT row-by-row, then column-by-column on
13// transposed result. Equivalent to multiplying by the DCT matrix on
14// both sides: Y = D X D^T.
15//
16// Fixed-point: every cosine entry is precomputed as Q10
17// (scale-by-1024) and stored in a caller-supplied 64-entry matrix
18// buffer M. After each row/column pass we shift right by 10 to
19// renormalise. Error budget is bounded by a few units in i64 at full
20// scale; the smoke verifies round-trip reconstruction within 8 units.
21//
22// genealogy_id: ahmed_natarajan_rao_1974_dct + chen_smith_fralick_1977 +
23// loeffler_ligtenberg_moschytz_1989_fast_dct +
24// jpeg_iso_10918_1992
25// lineage_id: nishi_dct8_q10
26
27// nx_safety_envelope:
28// intended_use: AUTO_APPLIED -- primitive-specific tuning queued
29// sil_target: SIL1
30// evidence: [bulk_applied_2026-05-16, see-file-comment-for-detail]
31// verdict: NOT_YET_EVALUATED
32
33import "nx_syscalls.nx"
34
35// Sealed verdict for DCT operations. Most paths are infallible
36// (no I/O), so this is mainly for caller-side audit logging.
37const NX_DCT8_VERDICT_UNKNOWN: i64 = 0
38const NX_DCT8_VERDICT_OK: i64 = 1
39const NX_DCT8_VERDICT_BAD_PARAMS: i64 = 2
40const NX_DCT8_VERDICT_N: i64 = 3
41
42// Initialise the 8x8 DCT matrix M (caller pre-allocates 64 i64 slots).
43// M[k * 8 + n] holds round(0.5 * c_k * cos((2n+1) k pi / 16) * 1024)
44// where c_0 = 1/sqrt(2) and c_k = 1 otherwise. Row 0 thus scales by
45// 1/(2*sqrt(2)) = 0.3536 -> 362 in Q10.
46func nx_dct8_init(M: *i64) -> i64 {
47 // Row k=0: all entries 362
48 M[0]=362; M[1]=362; M[2]=362; M[3]=362
49 M[4]=362; M[5]=362; M[6]=362; M[7]=362
50 // Row k=1: 0.5 * cos((2n+1) pi/16)
51 M[8]=502; M[9]=426; M[10]=284; M[11]=100
52 M[12]=-100; M[13]=-284; M[14]=-426; M[15]=-502
53 // Row k=2: 0.5 * cos((2n+1) pi/8)
54 M[16]=473; M[17]=196; M[18]=-196; M[19]=-473
55 M[20]=-473; M[21]=-196; M[22]=196; M[23]=473
56 // Row k=3
57 M[24]=426; M[25]=-100; M[26]=-502; M[27]=-284
58 M[28]=284; M[29]=502; M[30]=100; M[31]=-426
59 // Row k=4: 0.5 * cos((2n+1) pi/4) = +-0.3536
60 M[32]=362; M[33]=-362; M[34]=-362; M[35]=362
61 M[36]=362; M[37]=-362; M[38]=-362; M[39]=362
62 // Row k=5
63 M[40]=284; M[41]=-502; M[42]=100; M[43]=426
64 M[44]=-426; M[45]=-100; M[46]=502; M[47]=-284
65 // Row k=6
66 M[48]=196; M[49]=-473; M[50]=473; M[51]=-196
67 M[52]=-196; M[53]=473; M[54]=-473; M[55]=196
68 // Row k=7
69 M[56]=100; M[57]=-284; M[58]=426; M[59]=-502
70 M[60]=502; M[61]=-426; M[62]=284; M[63]=-100
71 return NX_DCT8_VERDICT_OK
72}
73
74// Forward 1D DCT-II on 8 input samples x[0..7] -> y[0..7].
75// y_k = sum_n M[k*8 + n] * x_n, then >> 10 to renormalise from Q10.
76func nx_dct8_forward_1d(M: *i64, x: *i64, y: *i64) -> i64 {
77 var k: i64 = 0
78 while k < 8 {
79 var acc: i64 = 0
80 var n: i64 = 0
81 while n < 8 {
82 acc = acc + M[k * 8 + n] * x[n]
83 n = n + 1
84 }
85 // Round half-away-from-zero before shift.
86 if acc >= 0 { y[k] = (acc + 512) >> 10 }
87 else { y[k] = -(((-acc) + 512) >> 10) }
88 k = k + 1
89 }
90 return NX_DCT8_VERDICT_OK
91}
92
93// Inverse 1D DCT on 8 frequency coefficients y[0..7] -> x[0..7].
94// x_n = sum_k M[k*8 + n] * y_k, then >> 10 (M is its own inverse up
95// to scaling because it is orthonormal scaled).
96func nx_dct8_inverse_1d(M: *i64, y: *i64, x: *i64) -> i64 {
97 var n: i64 = 0
98 while n < 8 {
99 var acc: i64 = 0
100 var k: i64 = 0
101 while k < 8 {
102 acc = acc + M[k * 8 + n] * y[k]
103 k = k + 1
104 }
105 if acc >= 0 { x[n] = (acc + 512) >> 10 }
106 else { x[n] = -(((-acc) + 512) >> 10) }
107 n = n + 1
108 }
109 return NX_DCT8_VERDICT_OK
110}
111
112// Forward 2D DCT on an 8x8 block (i64, row-major). Caller supplies
113// scratch (64 i64), tmp_in (8 i64), tmp_out (8 i64) so we don't leak
114// pages on hot paths.
115func nx_dct8_forward_2d(M: *i64, block: *i64, out: *i64,
116 scratch: *i64, tmp_in: *i64, tmp_out: *i64) -> i64 {
117 // Pass 1: 1D DCT row-by-row, write into scratch row by row.
118 var r: i64 = 0
119 while r < 8 {
120 var i: i64 = 0
121 while i < 8 {
122 tmp_in[i] = block[r * 8 + i]
123 i = i + 1
124 }
125 nx_dct8_forward_1d(M, tmp_in, tmp_out)
126 var j: i64 = 0
127 while j < 8 {
128 scratch[r * 8 + j] = tmp_out[j]
129 j = j + 1
130 }
131 r = r + 1
132 }
133 // Pass 2: 1D DCT column-by-column on scratch, write into out.
134 var c: i64 = 0
135 while c < 8 {
136 var i: i64 = 0
137 while i < 8 {
138 tmp_in[i] = scratch[i * 8 + c]
139 i = i + 1
140 }
141 nx_dct8_forward_1d(M, tmp_in, tmp_out)
142 var j: i64 = 0
143 while j < 8 {
144 out[j * 8 + c] = tmp_out[j]
145 j = j + 1
146 }
147 c = c + 1
148 }
149 return NX_DCT8_VERDICT_OK
150}
151
152// Inverse 2D DCT on an 8x8 frequency block -> spatial block.
153func nx_dct8_inverse_2d(M: *i64, coeffs: *i64, out: *i64,
154 scratch: *i64, tmp_in: *i64, tmp_out: *i64) -> i64 {
155 // Pass 1: 1D IDCT column-by-column.
156 var c: i64 = 0
157 while c < 8 {
158 var i: i64 = 0
159 while i < 8 {
160 tmp_in[i] = coeffs[i * 8 + c]
161 i = i + 1
162 }
163 nx_dct8_inverse_1d(M, tmp_in, tmp_out)
164 var j: i64 = 0
165 while j < 8 {
166 scratch[j * 8 + c] = tmp_out[j]
167 j = j + 1
168 }
169 c = c + 1
170 }
171 // Pass 2: 1D IDCT row-by-row.
172 var r: i64 = 0
173 while r < 8 {
174 var i: i64 = 0
175 while i < 8 {
176 tmp_in[i] = scratch[r * 8 + i]
177 i = i + 1
178 }
179 nx_dct8_inverse_1d(M, tmp_in, tmp_out)
180 var j: i64 = 0
181 while j < 8 {
182 out[r * 8 + j] = tmp_out[j]
183 j = j + 1
184 }
185 r = r + 1
186 }
187 return NX_DCT8_VERDICT_OK
188}
189
190// Sealed-enum validity gate.
191func nx_dct8_verdict_is_valid(v: i64) -> i64 {
192 if v < 0 { return 0 }
193 if v >= NX_DCT8_VERDICT_N { return 0 }
194 return 1
195}