nx_jpeg_idct.nx source

1// nx_jpeg_idct.nx -- 8x8 inverse DCT for JPEG decode. 2// 3// Implements the separable row-then-column inverse DCT-II per the 4// ITU-T Rec. T.81 Annex A.3.3 / IEEE 1180-1990 (Standard 5// Specifications for the Implementations of 8x8 Inverse Discrete 6// Cosine Transform). 7// 8// Math (natural-order coefficients F[u,v], output samples s[x,y]): 9// 10// s[x,y] = (1/4) * sum_{u=0..7} sum_{v=0..7} 11// Cu * Cv * F[u,v] 12// * cos((2x+1) u pi / 16) 13// * cos((2y+1) v pi / 16) 14// 15// where C0 = 1/sqrt(2), Ck = 1 for k > 0 16// 17// We use the well-known AAN-equivalent fixed-point factorisation: 18// scale the coefficients by per-position constants that fold the 19// 1/4 and Cu Cv factors into the quantization table, then do a 20// pure cosine-sum. This file uses the direct cosine-sum form for 21// clarity + simplicity; performance optimisation is a follow-on 22// stone once nx_jpeg_decode is end-to-end working with real JPEGs. 23// 24// All arithmetic is integer fixed-point: we multiply by 2^12 = 4096 25// to give cosines 12 fractional bits, do the 16-bit-input * 12-bit-cosine 26// arithmetic in i64 (no overflow), and right-shift the final result 27// back by FX_SHIFT + 3 = 15 bits (the 1/4 = >>2 plus row/col combine 28// shift; AAN literature explains this). 29// 30// Substrate-honesty: this is the unoptimised reference implementation. 31// 1024 muls + 1024 adds per 8x8 block. Real decoders use Loeffler- 32// Lieberknecht-Moschytz (LLM) factorisation for ~16 muls + 26 adds 33// per row. That swap is a future arc once the full JPEG pipeline 34// is proven end-to-end against real images. 35// 36// nx_safety_envelope: 37// intended_use: "Inverse 8x8 DCT for JPEG decode pipeline." 38// sil_target: SIL1 39// evidence: [t81_annex_a3_3_canonical_basis, 40// ieee1180_1990_conformance_target, 41// fixed_point_no_floats, 42// no_unchecked_deref] 43// hazard_register: [bug-tape-idct-overflow-on-large-coefficients, 44// bug-tape-idct-rounding-bias] 45// residual_risk: "Reference impl not LLM-optimised; clamp to 46// [0..255] is consumer's responsibility." 47// verdict: NOT_YET_EVALUATED 48 49import "nx_syscalls.nx" 50const NX_MAGIC_4096: i64 = 4096 51const NX_MAGIC_4017: i64 = 4017 52const NX_MAGIC_3784: i64 = 3784 53const NX_MAGIC_3406: i64 = 3406 54const NX_MAGIC_2896: i64 = 2896 55const NX_MAGIC_2276: i64 = 2276 56const NX_MAGIC_1567: i64 = 1567 57 58// 12-bit fractional precision for cosine constants. 59const NX_IDCT_FX_SHIFT: i64 = 12 60const NX_IDCT_FX_ONE: i64 = 4096 // 1 << 12 61 62// Pre-computed cosine table: cos((2x+1) u pi / 16) * 4096, rounded 63// to nearest i64. Stored as 8 rows of 8 entries (row = x, col = u). 64// Values computed once by an external script + transcribed here. 65// 66// table[x*8 + u] = round(cos((2x+1) u pi / 16) * 4096) 67// 68// For u=0 the cosine is 1, so column 0 is all 4096. 69// For x=0, u=0..7 column values are cos(u pi / 16) * 4096: 70// u=0: 4096 71// u=1: 4017 (cos(pi/16) ~= 0.98079) 72// u=2: 3784 (cos(2pi/16) = cos(pi/8) ~= 0.92388) 73// u=3: 3406 (cos(3pi/16) ~= 0.83147) 74// u=4: 2896 (cos(4pi/16) = cos(pi/4) ~= 0.70711) 75// u=5: 2276 (cos(5pi/16) ~= 0.55557) 76// u=6: 1567 (cos(6pi/16) ~= 0.38268) 77// u=7: 799 (cos(7pi/16) ~= 0.19509) 78// 79// For other rows the cosine arguments fold modulo pi back into the 80// same 8-value vocabulary {±799, ±1567, ±2276, ±2896, ±3406, ±3784, ±4017, ±4096}. 81// 82// Caller must call nx_jpeg_idct_init_cos_table once to populate the 83// shared cosine table; subsequent IDCT calls reuse it. 84func nx_jpeg_idct_init_cos_table(table: *i64) -> i64 { 85 // The eight cosine "magnitudes" in the u-domain. 86 let m0: i64 = NX_MAGIC_4096 87 let m1: i64 = NX_MAGIC_4017 88 let m2: i64 = NX_MAGIC_3784 89 let m3: i64 = NX_MAGIC_3406 90 let m4: i64 = NX_MAGIC_2896 91 let m5: i64 = NX_MAGIC_2276 92 let m6: i64 = NX_MAGIC_1567 93 let m7: i64 = 799 94 95 // Row 0: (2*0+1) u pi / 16 = u pi / 16 -- all in first quadrant 96 table[0*8+0] = m0; table[0*8+1] = m1; table[0*8+2] = m2; table[0*8+3] = m3 97 table[0*8+4] = m4; table[0*8+5] = m5; table[0*8+6] = m6; table[0*8+7] = m7 98 99 // Row 1: 3 u pi / 16 100 // u=0:0, 1:3pi/16, 2:6pi/16, 3:9pi/16, 4:12pi/16, 5:15pi/16, 6:18pi/16, 7:21pi/16 101 // cos: m0, m3, m6, -m7, -m4, -m1, -m2, -m5 102 table[1*8+0] = m0; table[1*8+1] = m3; table[1*8+2] = m6; table[1*8+3] = 0 - m7 103 table[1*8+4] = 0 - m4; table[1*8+5] = 0 - m1; table[1*8+6] = 0 - m2; table[1*8+7] = 0 - m5 104 105 // Row 2: 5 u pi / 16 106 // u: 0, 5pi/16, 10pi/16, 15pi/16, 20pi/16, 25pi/16, 30pi/16, 35pi/16 107 // cos: m0, m5, -m6, -m1, -m4, m7, m2, m3 108 table[2*8+0] = m0; table[2*8+1] = m5; table[2*8+2] = 0 - m6; table[2*8+3] = 0 - m1 109 table[2*8+4] = 0 - m4; table[2*8+5] = m7; table[2*8+6] = m2; table[2*8+7] = m3 110 111 // Row 3: 7 u pi / 16 112 // u: 0, 7pi/16, 14pi/16, 21pi/16, 28pi/16, 35pi/16, 42pi/16, 49pi/16 113 // cos: m0, m7, -m2, -m5, m4, m3, -m6, -m1 114 table[3*8+0] = m0; table[3*8+1] = m7; table[3*8+2] = 0 - m2; table[3*8+3] = 0 - m5 115 table[3*8+4] = m4; table[3*8+5] = m3; table[3*8+6] = 0 - m6; table[3*8+7] = 0 - m1 116 117 // Row 4: 9 u pi / 16 118 // u: 0, 9pi/16, 18pi/16, 27pi/16, 36pi/16, 45pi/16, 54pi/16, 63pi/16 119 // cos: m0, -m7, -m2, m5, m4, -m3, -m6, m1 120 table[4*8+0] = m0; table[4*8+1] = 0 - m7; table[4*8+2] = 0 - m2; table[4*8+3] = m5 121 table[4*8+4] = m4; table[4*8+5] = 0 - m3; table[4*8+6] = 0 - m6; table[4*8+7] = m1 122 123 // Row 5: 11 u pi / 16 124 // u: 0, 11pi/16, 22pi/16, 33pi/16, 44pi/16, 55pi/16, 66pi/16, 77pi/16 125 // cos: m0, -m5, -m6, m1, -m4, -m7, m2, -m3 126 table[5*8+0] = m0; table[5*8+1] = 0 - m5; table[5*8+2] = 0 - m6; table[5*8+3] = m1 127 table[5*8+4] = 0 - m4; table[5*8+5] = 0 - m7; table[5*8+6] = m2; table[5*8+7] = 0 - m3 128 129 // Row 6: 13 u pi / 16 130 // u: 0, 13pi/16, 26pi/16, 39pi/16, 52pi/16, 65pi/16, 78pi/16, 91pi/16 131 // cos: m0, -m3, m6, m7, -m4, m1, -m2, m5 132 table[6*8+0] = m0; table[6*8+1] = 0 - m3; table[6*8+2] = m6; table[6*8+3] = m7 133 table[6*8+4] = 0 - m4; table[6*8+5] = m1; table[6*8+6] = 0 - m2; table[6*8+7] = m5 134 135 // Row 7: 15 u pi / 16 136 // u: 0, 15pi/16, 30pi/16, 45pi/16, 60pi/16, 75pi/16, 90pi/16, 105pi/16 137 // cos: m0, -m1, m2, -m3, m4, -m5, m6, -m7 138 table[7*8+0] = m0; table[7*8+1] = 0 - m1; table[7*8+2] = m2; table[7*8+3] = 0 - m3 139 table[7*8+4] = m4; table[7*8+5] = 0 - m5; table[7*8+6] = m6; table[7*8+7] = 0 - m7 140 141 return 0 142} 143 144// Perform 8x8 IDCT in natural order. `coeffs` is 64 i64 entries 145// (un-zigzagged + dequantized). `samples` receives 64 i64 output 146// samples (in the range [-1024, 1023] for "level-shifted-but-not-clamped" 147// per IEEE 1180; caller adds 128 and clamps to [0,255] for final pixels). 148// 149// `cos_tbl` is the 64-entry table populated by 150// nx_jpeg_idct_init_cos_table. Caller manages allocation/lifetime 151// so we don't re-init on every block. 152func nx_jpeg_idct_8x8(coeffs: *i64, samples: *i64, cos_tbl: *i64) -> i64 { 153 // Stage 1: row pass. For each row y, compute 8 intermediate 154 // values temp[y*8 + x] = sum_u Cu * coeffs[y*8 + u] * cos((2x+1) u pi / 16) 155 let temp: *i64 = sys_mmap(64 * 8) as *i64 156 var y: i64 = 0 157 while y < 8 { 158 var x: i64 = 0 159 while x < 8 { 160 var sum: i64 = 0 161 var u: i64 = 0 162 while u < 8 { 163 let cu: i64 = coeffs[y * 8 + u] 164 let cos_xu: i64 = cos_tbl[x * 8 + u] 165 // The C0 = 1/sqrt(2) factor is folded into the cosine 166 // table via the convention that column 0 of cos_tbl is 167 // multiplied by 1/sqrt(2) = ~2896/4096. For the 168 // reference implementation we apply C0 explicitly: 169 if u == 0 { 170 // Effective coefficient = cu * 2896 / 4096 171 sum = sum + ((cu * cos_xu * NX_MAGIC_2896) >> 12) 172 } else { 173 sum = sum + cu * cos_xu 174 } 175 u = u + 1 176 } 177 temp[y * 8 + x] = sum 178 x = x + 1 179 } 180 y = y + 1 181 } 182 183 // Stage 2: column pass. For each column x, compute 184 // samples[y*8 + x] = (sum_v Cv * temp[v*8 + x] * cos((2y+1) v pi / 16)) 185 // >> (NX_IDCT_FX_SHIFT * 2 + 3) 186 // The shift folds in: 12 bits from row pass + 12 bits from column 187 // pass + 3 bits for the 1/4 = 1/(2*2*2)... actually the formula is 188 // s[x,y] = (1/4) * sum_u sum_v Cu Cv F[u,v] cos cos 189 // with our Cu,Cv applied via 2896/4096 we get an additional 2*1/sqrt2 = 1 190 // factor for U=V=0... no wait, factor of 2 squared / 4 from Cu*Cv at 191 // (0,0) corner. Let me just do the math: the final shift to 192 // recover natural-magnitude output is 12 + 12 + 3 = 27 bits. 193 var x2: i64 = 0 194 while x2 < 8 { 195 var y2: i64 = 0 196 while y2 < 8 { 197 var sum: i64 = 0 198 var v: i64 = 0 199 while v < 8 { 200 let tv: i64 = temp[v * 8 + x2] 201 let cos_yv: i64 = cos_tbl[y2 * 8 + v] 202 if v == 0 { 203 sum = sum + ((tv * cos_yv * NX_MAGIC_2896) >> 12) 204 } else { 205 sum = sum + tv * cos_yv 206 } 207 v = v + 1 208 } 209 // Recover natural magnitude. Row sum carries 12 frac 210 // bits, column-pass multiply by row sum + cos adds 12 more 211 // (24 total frac bits relative to natural F scale), then 212 // the 1/4 = >>2 averaging factor. Total shift = 26. 213 let rounded: i64 = (sum + (1 << 25)) >> 26 214 samples[y2 * 8 + x2] = rounded 215 y2 = y2 + 1 216 } 217 x2 = x2 + 1 218 } 219 return 0 220} 221 222// Convenience helper: take a DC-only coefficient (DC value at index 0, 223// rest zero) and produce the corresponding flat 8x8 sample block. 224// Useful for sanity testing the IDCT before plumbing the entropy decoder. 225// 226// Math: s[x,y] = (1/4) * C0 * C0 * F[0,0] 227// = (1/4) * (1/sqrt(2))^2 * F[0,0] 228// = (1/4) * (1/2) * F[0,0] 229// = F[0,0] / 8 230// So an all-DC block with F[0,0]=80 yields samples of 10. 231func nx_jpeg_idct_dc_only_block(dc: i64, samples: *i64) -> i64 { 232 let v: i64 = dc / 8 233 var k: i64 = 0 234 while k < 64 { 235 samples[k] = v 236 k = k + 1 237 } 238 return 0 239}