nx_hevc_idct.nx source

1// nx_hevc_idct.nx -- SOVEREIGN HEVC dequant + inverse transform, RUNG 5 of the HEVC decoder. 2// Dequant (8.6.3): d = (level*16*levelScale[qP%6] << (qP/6) + round) >> bdShift, clipped. Inverse transform 3// (8.6.4): 2-stage (vertical then horizontal) integer DCT-II (4/8/16/32) or DST-VII (4x4 luma intra), with 4// the normative basis matrices. Self-test: DC-only coeff -> flat residual (hand-verified). R5a = DST4/DCT4/DCT8; 5// DCT16/32 (same mechanism, bigger matrices) = R5b. No third party. license_tier: ORIGINAL 6import "nx_syscalls.nx" 7const K_MAGIC_32768: i64 = 32768 8const K_MAGIC_32767: i64 = 32767 9const K_MAGIC_1000000: i64 = 1000000 10const K_MAGIC_3141593: i64 = 3141593 11const K_MAGIC_90509678: i64 = 90509678 12const K_MAGIC_500000000000: i64 = 500000000000 13const K_MAGIC_1000000000000: i64 = 1000000000000 14const K_MAGIC_1024: i64 = 1024 15 16func pe(s: *u8) -> i64 { var n: i64=0; while s[n]!=(0 as u8){n=n+1} return sys_write(1,s,n) } 17func pn(v: i64) -> i64 { var m: i64=v; if m<0{m=0-m} let b: *u8=sys_mmap(32); var i: i64=32; if m==0{i=i-1;b[i]=(48 as u8)} while m>0{let q: i64=m/10; i=i-1; b[i]=((48+(m-q*10)) as u8); m=q} if v<0{i=i-1;b[i]=(45 as u8)} return sys_write(1,((b as i64)+i) as *u8,32-i) } 18func parse_nums(s: *u8, out: *i64) -> i64 { var n: i64=0; var i: i64=0; var cur: i64=0; var neg: i64=0; var has: i64=0; while s[i]!=(0 as u8){ let c: i64=s[i] as i64; if c==45 {neg=1} else { if c>=48 { if c<=57 {cur=cur*10+(c-48); has=1} else { if has==1 { if neg==1 {out[n]=0-cur} else {out[n]=cur} n=n+1; cur=0; has=0; neg=0 } } } else { if has==1 { if neg==1 {out[n]=0-cur} else {out[n]=cur} n=n+1; cur=0; has=0; neg=0 } } } i=i+1 } if has==1 { if neg==1 {out[n]=0-cur} else {out[n]=cur} n=n+1 } return n } 19func clip3(lo: i64, hi: i64, v: i64) -> i64 { if v<lo {return lo} if v>hi {return hi} return v } 20func ilog2(n: i64) -> i64 { var l: i64=0; var v: i64=n; while v>1 {v=v>>1; l=l+1} return l } 21 22// dequant a block in place. levelScale={40,45,51,57,64,72}, m=16 (flat scaling list). 23func dequant(coeff: *i64, N: i64, qP: i64, bitDepth: i64, lsc: *i64) -> i64 { 24 let bdShift: i64 = bitDepth + ilog2(N) - 5 25 let ls: i64 = lsc[qP%6]; let shift: i64 = qP/6; let add: i64 = 1<<(bdShift-1) 26 var i: i64=0; while i<N*N { coeff[i] = clip3(0-K_MAGIC_32768, K_MAGIC_32767, ((coeff[i]*16*ls << shift) + add) >> bdShift); i=i+1 } 27 return 0 28} 29// 2-stage inverse transform in place. mat = NxN basis (mat[k*N+n]). residual in coeff. 30func itransform(coeff: *i64, N: i64, mat: *i64, bitDepth: i64, tmp: *i64) -> i64 { 31 // stage 1: vertical (each column), shift 7, clip 16-bit 32 var c: i64=0 33 while c<N { var r: i64=0; while r<N { var s: i64=0; var k: i64=0; while k<N { s=s+mat[k*N+r]*coeff[k*N+c]; k=k+1 } tmp[r*N+c]=clip3(0-K_MAGIC_32768,K_MAGIC_32767,(s+64)>>7); r=r+1 } c=c+1 } 34 // stage 2: horizontal (each row), shift (20-bitDepth) 35 let bd2: i64 = 20-bitDepth; let add2: i64 = 1<<(bd2-1) 36 var r2: i64=0 37 while r2<N { var c2: i64=0; while c2<N { var s: i64=0; var k: i64=0; while k<N { s=s+mat[k*N+c2]*tmp[r2*N+k]; k=k+1 } coeff[r2*N+c2]=(s+add2)>>bd2; c2=c2+1 } r2=r2+1 } 38 return 0 39} 40 41// build DCT_N from DCT_(N/2): even rows 2j = gHalf[j] mirrored (symmetric); odd rows 2j+1 = oddC[j] antisym. 42func build_dct(gN: *i64, N: i64, gHalf: *i64, oddC: *i64) -> i64 { 43 let H: i64=N/2; var j: i64=0 44 while j<H { var nn: i64=0; while nn<H { 45 gN[(2*j)*N+nn]=gHalf[j*H+nn]; gN[(2*j)*N+(N-1-nn)]=gHalf[j*H+nn] 46 gN[(2*j+1)*N+nn]=oddC[j*H+nn]; gN[(2*j+1)*N+(N-1-nn)]=0-oddC[j*H+nn] 47 nn=nn+1 } j=j+1 } 48 return 0 49} 50 51// fixed-point cos(m*pi/D) scaled to round(64*sqrt2*cos) = a HEVC DCT odd-basis integer. Range-reduce to 52// [0,pi/2] (sign-tracked), then Taylor (term-by-term, no overflow). 64*sqrt2 = 90.509678. 53func fixcos(m: i64, D: i64) -> i64 { 54 let S: i64=K_MAGIC_1000000; let PI: i64=K_MAGIC_3141593 55 var mm: i64 = m%(2*D); if mm<0 {mm=mm+2*D} 56 if mm>D {mm=2*D-mm} 57 var sign: i64=1; if 2*mm>D {mm=D-mm; sign=0-1} 58 let x: i64 = mm*PI/D; let X2: i64 = x*x/S 59 var cosf: i64=S; var term: i64=S; var k: i64=1 60 while k<=10 { term=(0-term)*X2/((2*k)*(2*k-1))/S; cosf=cosf+term; k=k+1 } 61 let mag: i64 = (K_MAGIC_90509678*cosf + K_MAGIC_500000000000)/K_MAGIC_1000000000000 62 return sign*mag 63} 64// odd-constant table for the N-point DCT: oc[j][n] = round(64sqrt2 cos((2n+1)(2j+1)pi/D)), H=N/2, D=2N. 65func gen_oc(oc: *i64, H: i64, D: i64) -> i64 { var j: i64=0; while j<H { var n: i64=0; while n<H { oc[j*H+n]=fixcos((2*n+1)*(2*j+1), D); n=n+1 } j=j+1 } return 0 } 66 67func main(argc: i64, argv: *i64) -> i64 { 68 let lsc: *i64=sys_mmap(8*8) as *i64; parse_nums("40 45 51 57 64 72" as *u8, lsc) 69 let dct4: *i64=sys_mmap(8*16) as *i64; parse_nums("64 64 64 64 83 36 -36 -83 64 -64 -64 64 36 -83 83 -36" as *u8, dct4) 70 let dst4: *i64=sys_mmap(8*16) as *i64; parse_nums("29 55 74 84 74 74 0 -74 84 -29 -74 55 55 -84 74 -29" as *u8, dst4) 71 let dct8: *i64=sys_mmap(8*64) as *i64; parse_nums("64 64 64 64 64 64 64 64 89 75 50 18 -18 -50 -75 -89 83 36 -36 -83 -83 -36 36 83 75 -18 -89 -50 50 89 18 -75 64 -64 -64 64 64 -64 -64 64 50 -89 18 75 -75 -18 89 -50 36 -83 83 -36 -36 83 -83 36 18 -50 75 -89 89 -75 50 -18" as *u8, dct8) 72 let coeff: *i64=sys_mmap(8*K_MAGIC_1024) as *i64; let tmp: *i64=sys_mmap(8*K_MAGIC_1024) as *i64 73 // test 1: DC-only 4x4, level 10, qP 24, 10-bit -> dequant d[0]=800, IDCT -> flat 25 (hand-computed) 74 var i: i64=0; while i<16 {coeff[i]=0; i=i+1} coeff[0]=10 75 dequant(coeff, 4, 24, 10, lsc) 76 pe("DC-only: dequant d[0]=" as *u8); pn(coeff[0]); pe(" (expect 800)\n" as *u8) 77 itransform(coeff, 4, dct4, 10, tmp) 78 var flat: i64=1; i=0; while i<16 {if coeff[i]!=coeff[0] {flat=0} i=i+1} 79 pe("DC-only IDCT4: residual[0]=" as *u8); pn(coeff[0]); pe(" flat=" as *u8); pn(flat); pe(" (expect 25, flat=1)\n" as *u8) 80 // test 2: DST4 DC-only flat check (level 10 qP24) -> flat (DST is not constant-basis, so DC-only is NOT flat; just run) 81 i=0; while i<16 {coeff[i]=0; i=i+1} coeff[0]=10; dequant(coeff,4,24,10,lsc); itransform(coeff,4,dst4,10,tmp) 82 pe("DST4 DC-only: residual[0]=" as *u8); pn(coeff[0]); pe(" residual[15]=" as *u8); pn(coeff[15]); pe(" (DST basis, not flat)\n" as *u8) 83 // test 3: DCT8 DC-only -> flat 84 i=0; while i<64 {coeff[i]=0; i=i+1} coeff[0]=10; dequant(coeff,8,24,10,lsc); itransform(coeff,8,dct8,10,tmp) 85 var flat8: i64=1; i=0; while i<64 {if coeff[i]!=coeff[0] {flat8=0} i=i+1} 86 pe("DCT8 DC-only: residual[0]=" as *u8); pn(coeff[0]); pe(" flat=" as *u8); pn(flat8); pe("\n" as *u8) 87 var ok: i64=1; if coeff[0]==0 {ok=0} if flat==0 {ok=0} if flat8==0 {ok=0} 88 if ok==1 { pe("R5a OK: dequant + inverse DCT4/DST4/DCT8 working (DC-only -> flat for DCT, hand-verified)\n" as *u8) } else { pe("R5a check: flat4=" as *u8); pn(flat); pe(" flat8=" as *u8); pn(flat8); pe("\n" as *u8) } 89 // ---- R5b: build DCT16/DCT32 via even/odd butterfly ---- 90 // exact spec odd-constant tables (HEVC HM, with intentional norm-preserving +-1 deviations from round(64sqrt2 cos)) 91 let oc8: *i64=sys_mmap(8*16) as *i64; parse_nums("89 75 50 18 75 -18 -89 -50 50 -89 18 75 18 -50 75 -89" as *u8, oc8) 92 let oc16: *i64=sys_mmap(8*64) as *i64; parse_nums("90 87 80 70 57 43 25 9 87 57 9 -43 -80 -90 -70 -25 80 9 -70 -87 -25 57 90 43 70 -43 -87 9 90 25 -80 -57 57 -80 -25 90 -9 -87 43 70 43 -90 57 25 -87 70 9 -80 25 -70 90 -80 43 9 -57 87 9 -25 43 -57 70 -80 87 -90" as *u8, oc16) 93 let oc32: *i64=sys_mmap(8*256) as *i64; parse_nums("90 90 88 85 82 78 73 67 61 54 46 38 31 22 13 4 90 82 67 46 22 -4 -31 -54 -73 -85 -90 -88 -78 -61 -38 -13 88 67 31 -13 -54 -82 -90 -78 -46 -4 38 73 90 85 61 22 85 46 -13 -67 -90 -73 -22 38 82 88 54 -4 -61 -90 -78 -31 82 22 -54 -90 -61 13 78 85 31 -46 -90 -67 4 73 88 38 78 -4 -82 -73 13 85 67 -22 -88 -61 31 90 54 -38 -90 -46 73 -31 -90 -22 78 67 -38 -90 -13 82 61 -46 -88 -4 85 54 67 -54 -78 38 85 -22 -90 4 90 13 -88 -31 82 46 -73 -61 61 -73 -46 82 31 -88 -13 90 -4 -90 22 85 -38 -78 54 67 54 -85 -4 88 -46 -61 82 13 -90 38 67 -78 -22 90 -31 -73 46 -90 38 54 -90 31 61 -88 22 67 -85 13 73 -82 4 78 38 -88 73 -4 -67 90 -46 -31 85 -78 13 61 -90 54 22 -82 31 -78 90 -61 4 54 -88 82 -38 -22 73 -90 67 -13 -46 85 22 -61 85 -90 73 -38 -4 46 -78 90 -82 54 -13 -31 67 -88 13 -38 61 -78 88 -90 85 -73 54 -31 4 22 -46 67 -82 90 4 -13 22 -31 38 -46 54 -61 67 -73 78 -82 85 -88 90 -90" as *u8, oc32) 94 let g8b: *i64=sys_mmap(8*64) as *i64; build_dct(g8b, 8, dct4, oc8) 95 var d8: i64=0; var z: i64=0; while z<64 {if g8b[z]!=dct8[z] {d8=d8+1} z=z+1} 96 pe("DCT8 from oc8 vs embedded dct8: diffs=" as *u8); pn(d8); pe(" (0=oc8 exact)\n" as *u8) 97 var tdiff: i64=0; var ti: i64=0 98 while ti<64 { let f: i64=fixcos((2*(ti%8)+1)*(2*(ti/8)+1),32); var dd: i64=oc16[ti]-f; if dd<0{dd=0-dd} if dd>1{tdiff=tdiff+1} ti=ti+1 } 99 ti=0; while ti<256 { let f: i64=fixcos((2*(ti%16)+1)*(2*(ti/16)+1),64); var dd: i64=oc32[ti]-f; if dd<0{dd=0-dd} if dd>1{tdiff=tdiff+1} ti=ti+1 } 100 pe("oc16/oc32 vs fixcos: entries off by >1 (typos)=" as *u8); pn(tdiff); pe(" (0=no typos; legit norm-deviations are <=1)\n" as *u8) 101 let g16: *i64=sys_mmap(8*256) as *i64; build_dct(g16, 16, g8b, oc16) 102 let g32: *i64=sys_mmap(8*K_MAGIC_1024) as *i64; build_dct(g32, 32, g16, oc32) 103 // HEVC integer DCT is NEAR-orthogonal (small nonzero cross-products by design). Verify via known rows + bounded cross. 104 pe("g16[1]= " as *u8); var q: i64=0; while q<8 {pn(g16[16+q]); pe(" " as *u8); q=q+1} pe("(expect 90 87 80 70 57 43 25 9)\n" as *u8) 105 pe("g32[1]= " as *u8); q=0; while q<8 {pn(g32[32+q]); pe(" " as *u8); q=q+1} pe("(expect 90 90 88 85 82 78 73 67)\n" as *u8) 106 var kn: i64=1 107 if g16[16]!=90 {kn=0} if g16[16+1]!=87 {kn=0} if g16[16+7]!=9 {kn=0} 108 if g32[32]!=90 {kn=0} if g32[32+2]!=88 {kn=0} if g32[32+7]!=67 {kn=0} 109 var maxc: i64=0; var k1: i64=0 110 while k1<32 { var k2: i64=k1+1; while k2<32 { var s: i64=0; var nn: i64=0; while nn<32 {s=s+g32[k1*32+nn]*g32[k2*32+nn]; nn=nn+1} if s<0 {s=0-s} if s>maxc {maxc=s} k2=k2+1 } k1=k1+1 } 111 pe("DCT32 max|cross-product|=" as *u8); pn(maxc); pe(" (vs diagonal 131072 -> near-orthogonal)\n" as *u8) 112 i=0; while i<K_MAGIC_1024 {coeff[i]=0; i=i+1} coeff[0]=10; dequant(coeff,32,24,10,lsc); itransform(coeff,32,g32,10,tmp) 113 var flat32: i64=1; i=0; while i<K_MAGIC_1024 {if coeff[i]!=coeff[0] {flat32=0} i=i+1} 114 pe("DCT32 DC-only: residual[0]=" as *u8); pn(coeff[0]); pe(" flat=" as *u8); pn(flat32); pe("\n" as *u8) 115 if d8==0 { if tdiff==0 { if kn==1 { if flat32==1 { pe("R5b OK: exact spec DCT16/32 (oc8 cross-check + fixcos typo-detect + known rows + near-orthogonal + DC-flat) -> all inverse transforms ready\n" as *u8) } } } } 116 return 0 117}