nx_av1_dct.nx source

1// nx_av1_dct.nx -- AV1/AV2 inverse DCT butterfly networks, 4 and 8 point. 2// 3// The residual path's core. Built on nx_av1_txfm.nx's 12-bit cosine table and 4// Round2, so every constant traces to one definition rather than being baked 5// in as a magic number. 6// 7// ANGLE MAPPING. The classic references express these constants as 8// cospi_N_64 = cos(N*pi/64). This table is a 256-step circle, so cos128(a) is 9// cos(a*pi/128) and the mapping is a = 2N: cospi_16_64 -> cos128(32), 10// cospi_8_64 -> cos128(16), cospi_24_64 -> cos128(48), cospi_4_64 -> 11// cos128(8), cospi_28_64 -> cos128(56). Halving that factor is the obvious 12// mistake and it produces a transform that still looks like a DCT -- the 13// basis functions are simply at the wrong frequencies, so the image decodes 14// blurred rather than broken. 15// 16// WHAT IS PROVEN HERE. The gate checks STRUCTURAL properties -- a DC-only 17// input must give a flat output, an odd basis index must give an 18// antisymmetric one, and the transform must be linear within rounding. These 19// catch a mis-wired butterfly, a wrong angle and a sign error. They do NOT 20// establish bit-exactness against libaom's reference vectors, which needs 21// those vectors; that is a separate, later step and is not claimed. 22// 23// genealogy_id: av1_spec_7_13_2_inverse_dct 24// lineage_id: nx_av1_dct_v1 25// license_tier: ORIGINAL 26 27import "nx_syscalls.nx" 28import "nx_av1_txfm.nx" 29const K_MAGIC_3784: i64 = 3784 30const K_MAGIC_1567: i64 = 1567 31 32// ===== 4-point inverse DCT ======================================== 33// 34// Operates in place on four coefficients. 35 36func nx_av1_idct4(t: *i64, io: *i64) -> i64 { 37 // The odd pair rotates by cospi_24_64 and cospi_8_64 IN THAT ORDER. 38 // Swapping them keeps the output antisymmetric -- so the symmetry test 39 // still passes -- while putting the basis at the wrong frequency, which 40 // shows up only as a monotonicity failure. That swap was a live bug here, 41 // caught by T2's ordering assertion rather than its symmetry one. 42 let c32: i64 = nx_av1_cos128(t, 32) // cospi_16_64, ~1/sqrt2 43 let c8: i64 = nx_av1_cos128(t, 16) // cospi_8_64 = K_MAGIC_3784 44 let c24: i64 = nx_av1_cos128(t, 48) // cospi_24_64 = K_MAGIC_1567 45 46 let s0: i64 = nx_av1_round2((io[0] + io[2]) * c32, NX_TX_COS_BITS) 47 let s1: i64 = nx_av1_round2((io[0] - io[2]) * c32, NX_TX_COS_BITS) 48 let s2: i64 = nx_av1_round2(io[1] * c24 - io[3] * c8, NX_TX_COS_BITS) 49 let s3: i64 = nx_av1_round2(io[1] * c8 + io[3] * c24, NX_TX_COS_BITS) 50 51 io[0] = s0 + s3 52 io[1] = s1 + s2 53 io[2] = s1 - s2 54 io[3] = s0 - s3 55 return 1 56} 57 58// ===== 8-point inverse DCT ======================================== 59// 60// Four stages. Stage 2 and stage 3 both contain the 1/sqrt2 rotation, which 61// is where a halved angle shows up first. 62 63func nx_av1_idct8(t: *i64, io: *i64) -> i64 { 64 let c32: i64 = nx_av1_cos128(t, 32) 65 let c56: i64 = nx_av1_cos128(t, 56) 66 let s56: i64 = nx_av1_sin128(t, 56) 67 let c24: i64 = nx_av1_cos128(t, 24) 68 let s24: i64 = nx_av1_sin128(t, 24) 69 let c48: i64 = nx_av1_cos128(t, 48) 70 let s48: i64 = nx_av1_sin128(t, 48) 71 72 let a: *i64 = sys_mmap(16 * 8) as *i64 73 let b: *i64 = sys_mmap(16 * 8) as *i64 74 75 // ---- stage 1: even terms pass through, odd terms rotate ---- 76 a[0] = io[0] 77 a[1] = io[4] 78 a[2] = io[2] 79 a[3] = io[6] 80 a[4] = nx_av1_round2(io[1] * c56 - io[7] * s56, NX_TX_COS_BITS) 81 a[7] = nx_av1_round2(io[1] * s56 + io[7] * c56, NX_TX_COS_BITS) 82 a[5] = nx_av1_round2(io[5] * c24 - io[3] * s24, NX_TX_COS_BITS) 83 a[6] = nx_av1_round2(io[5] * s24 + io[3] * c24, NX_TX_COS_BITS) 84 85 // ---- stage 2 ---- 86 b[0] = nx_av1_round2((a[0] + a[1]) * c32, NX_TX_COS_BITS) 87 b[1] = nx_av1_round2((a[0] - a[1]) * c32, NX_TX_COS_BITS) 88 b[2] = nx_av1_round2(a[2] * c48 - a[3] * s48, NX_TX_COS_BITS) 89 b[3] = nx_av1_round2(a[2] * s48 + a[3] * c48, NX_TX_COS_BITS) 90 b[4] = a[4] + a[5] 91 b[5] = a[4] - a[5] 92 b[6] = a[7] - a[6] 93 b[7] = a[6] + a[7] 94 95 // ---- stage 3 ---- 96 a[0] = b[0] + b[3] 97 a[1] = b[1] + b[2] 98 a[2] = b[1] - b[2] 99 a[3] = b[0] - b[3] 100 a[4] = b[4] 101 a[5] = nx_av1_round2((b[6] - b[5]) * c32, NX_TX_COS_BITS) 102 a[6] = nx_av1_round2((b[5] + b[6]) * c32, NX_TX_COS_BITS) 103 a[7] = b[7] 104 105 // ---- stage 4 ---- 106 io[0] = a[0] + a[7] 107 io[1] = a[1] + a[6] 108 io[2] = a[2] + a[5] 109 io[3] = a[3] + a[4] 110 io[4] = a[3] - a[4] 111 io[5] = a[2] - a[5] 112 io[6] = a[1] - a[6] 113 io[7] = a[0] - a[7] 114 return 1 115} 116 117// ===== the array permutation ====================================== 118// 119// The spec feeds the butterfly network a bit-reversed coefficient order. 120 121func nx_av1_dct_permute(io: *i64, n_log2: i64) -> i64 { 122 if n_log2 <= 0 { return 0 } 123 let n: i64 = 1 << n_log2 124 let copy: *i64 = sys_mmap(n * 8 + 64) as *i64 125 var i: i64 = 0 126 while i < n { copy[i] = io[i]; i = i + 1 } 127 i = 0 128 while i < n { 129 io[i] = copy[nx_av1_brev(n_log2, i)] 130 i = i + 1 131 } 132 return 1 133}