oxideav-av1 0.1.21

//! Forward 2D transform dispatcher — the encoder counterpart of the
//! §7.13.3 2D inverse-transform dispatcher implemented in
//! [`crate::transform::inverse_transform_2d`].
//!
//! ## Why "forward dispatcher"
//!
//! Rounds 219, 222, 225, 226 built the per-axis forward 1D / 2D
//! primitives in [`super::forward_transform`] (DCT, sizes 4..64),
//! [`super::forward_adst`] (ADST + FLIPADST, sizes 4 / 8 / 16),
//! [`super::forward_identity`] (IDTX, sizes 4 / 8 / 16 / 32), and
//! [`super::forward_wht`] (WHT, size 4 only — the lossless arm).
//!
//! Those primitives are the forward kernels for a **fixed**
//! row / column kernel choice. The §7.13.3 inverse dispatcher
//! composes a per-(`PlaneTxType`, `tx_size`) row kernel with a
//! column kernel and wraps the composition in the §7.13.3 shift
//! envelope (`Round2(_, rowShift)` per row after the row pass,
//! `Round2(_, colShift)` per column after the column pass) plus a
//! between-stage `Clip3` clamp at `colClampRange` bits. The encoder
//! mirror — this module — composes the **forward** kernels in the
//! transpose order (column-then-row) and leaves the shift envelope
//! to be absorbed by the decoder's `Round2(_, shift)` post-scales
//! against the kernel's intrinsic `N`-times gain (see "Shift
//! envelope" below for the derivation).
//!
//! ## Square scope plus the four `|log2W - log2H| == 1` pairs and
//! two `|log2W - log2H| == 2` pairs (rounds 250 + 251)
//!
//! The five **square** transform sizes are landed:
//! `TX_4X4`, `TX_8X8`, `TX_16X16`, `TX_32X32`, `TX_64X64`. Round 235
//! added the **`|log2W - log2H| == 1`** rectangular pair at the
//! shorter side — specifically `TX_4X8` and `TX_8X4`. Round 238
//! extended the `|log2W - log2H| == 1` arc by one more pair —
//! `TX_8X16` and `TX_16X8` — covering the `min(log2W, log2H) == 3`
//! short-side-8 family. Round 241 extended the arc by the
//! short-side-16 pair `TX_16X32` / `TX_32X16` (`min_log2 == 4`).
//! Round 244 closed the `|log2W - log2H| == 1` chain with the
//! largest-kernel pair — **`TX_32X64`** and **`TX_64X32`**, at
//! `min(log2W, log2H) == 5`. Round 250 opened the second arc — the
//! **`|log2W - log2H| == 2`** family — at the pair with the
//! smallest spec-aligned delta from the r244 chain: **`TX_16X64`**
//! and **`TX_64X16`**, at `min(log2W, log2H) == 4`, `max(log2W,
//! log2H) == 6`. Round 251 extends the `|log2W - log2H| == 2` arc
//! at the opposite end — the smallest-kernel pair at `min(log2W,
//! log2H) == 2`, `max(log2W, log2H) == 4`: **`TX_4X16`** and
//! **`TX_16X4`**. Round 252 closes the `|log2W - log2H| == 2`
//! family with its remaining mid-pair at `min(log2W, log2H) == 3`,
//! `max(log2W, log2H) == 5`: **`TX_8X32`** and **`TX_32X8`**.
//!
//! The three `|log2W - log2H| == 2` shapes landed in rounds
//! 250 / 251 / 252 split into two reachability regimes by §5.11.40
//! `compute_tx_type()`. `TX_16X64` / `TX_64X16` have
//! `txSzSqrUp == TX_64X64`, so they are **forced to `DCT_DCT`** by
//! the `txSzSqrUp > TX_32X32 ⇒ return DCT_DCT` short-circuit,
//! leaving only one reachable kernel pair. `TX_4X16` / `TX_16X4`
//! have `txSzSqrUp == TX_16X16` (≤ `TX_32X32`) so the §5.11.40
//! short-circuit does NOT fire — and the §5.11.48 `get_tx_set()`
//! routes through `TX_SET_INTRA_1` / `TX_SET_INTER_1` (their
//! `txSzSqr == TX_4X4`), giving the full §6.10.19 16-ordinal
//! tx_type matrix. `TX_8X32` / `TX_32X8` have
//! `txSzSqrUp == TX_32X32`, so §5.11.40's `txSzSqrUp > TX_32X32`
//! short-circuit does NOT fire (strict inequality), but §5.11.48's
//! intra arm `txSzSqrUp == TX_32X32 ⇒ TX_SET_DCTONLY` does — so
//! intra-coded blocks of this shape see only `DCT_DCT`. The inter
//! arm of §5.11.48 routes `txSzSqrUp == TX_32X32` (or
//! `reduced_tx_set`) to `TX_SET_INTER_3`, whose §6.10.19 inversion
//! table `Tx_Type_Inter_Inv_Set3` = `{IDTX, DCT_DCT}` admits exactly
//! two ordinals. So the reachable per-shape tx_type set on
//! `TX_8X32` / `TX_32X8` is `{DCT_DCT, IDTX}` — two ordinals, not
//! the 16-ordinal full matrix the smaller pair allows. Each axis
//! runs its own log-size forward kernel:
//!
//! | tx_size    | row kernel        | col kernel        |
//! | ---------- | ----------------- | ----------------- |
//! | `TX_4X8`   | `forward_dct_4`   | `forward_dct_8`   |
//! | `TX_8X4`   | `forward_dct_8`   | `forward_dct_4`   |
//! | `TX_8X16`  | `forward_dct_8`   | `forward_dct_16`  |
//! | `TX_16X8`  | `forward_dct_16`  | `forward_dct_8`   |
//! | `TX_16X32` | `forward_dct_16`  | `forward_dct_32`  |
//! | `TX_32X16` | `forward_dct_32`  | `forward_dct_16`  |
//! | `TX_32X64` | `forward_dct_32`  | `forward_dct_64`  |
//! | `TX_64X32` | `forward_dct_64`  | `forward_dct_32`  |
//! | `TX_16X64` | `forward_dct_16`  | `forward_dct_64`  |
//! | `TX_64X16` | `forward_dct_64`  | `forward_dct_16`  |
//! | `TX_4X16`  | `forward_dct_4`   | `forward_dct_16`  |
//! | `TX_16X4`  | `forward_dct_16`  | `forward_dct_4`   |
//! | `TX_8X32`  | `forward_dct_8`   | `forward_dct_32`  |
//! | `TX_32X8`  | `forward_dct_32`  | `forward_dct_8`   |
//!
//! (kernel column shown for the `DCT_DCT` arm; the ADST / IDTX
//! arms route through the per-axis ADST / IDTX dispatcher at the
//! same `log2_w` / `log2_h` slots — see `forward_row_kernel` /
//! `forward_col_kernel` below).
//!
//! For the `|log2W - log2H| == 1` family the per-row §7.13.3
//! rectangular scale `Round2(T[j] * 2896, 12)` is applied after the
//! row kernel pass — the encoder mirror of the decoder's pre-row-
//! kernel rectangular scale at the same constant. Both encoder and
//! decoder thus contribute one factor of `2896 / 4096`, giving a
//! combined rectangular gain of `(2896 / 4096)^2 ≈ 1/2` per
//! rectangular axis pair. This holds uniformly across the
//! `|log2W - log2H| == 1` family — the per-tx-size differentiator
//! is the row-shift envelope (`Transform_Row_Shift[TX_8X16] =
//! Transform_Row_Shift[TX_16X8] = Transform_Row_Shift[TX_16X32] =
//! Transform_Row_Shift[TX_32X16] = Transform_Row_Shift[TX_32X64] =
//! Transform_Row_Shift[TX_64X32] = 1` vs the `0` for `TX_4X8` /
//! `TX_8X4`), not the rectangular scale constant itself. The
//! `TX_32X64` / `TX_64X32` pair also pulls the §7.12.3
//! `dqDenom = 4` branch (any 64-axis transform);
//! [`super::forward_quantize`] already routes that path via
//! [`crate::cdf::dequant_denom`].
//!
//! For the `|log2W - log2H| == 2` family (the `TX_16X64` /
//! `TX_64X16` pair landed in round 250, the `TX_4X16` / `TX_16X4`
//! pair landed in round 251, and the `TX_8X32` / `TX_32X8` pair
//! landed in round 252) the `× 2896` rectangular row-pass step
//! does NOT fire — per av1-spec p.305 §7.13.3 line "if (Abs(log2W -
//! log2H) == 1) T[j] = Round2(T[j] * 2896, 12)", the rectangular
//! post-scale gates on **exactly** `|log2W - log2H| == 1`. This
//! family runs the bare per-axis kernel composition without the
//! rectangular `2896` factor on either side; the per-cell round-trip
//! scale is the analytic `(N_w * N_h) / 2^(row_shift + col_shift)`
//! ratio, undamped. For the `TX_16X64` / `TX_64X16` pair `N_w * N_h
//! = 16 * 64 = 1024` and `Transform_Row_Shift = 2`, so the analytic
//! scale is `1024 / 2^(2 + 4) = 16`. For the `TX_4X16` / `TX_16X4`
//! pair (round 251) `N_w * N_h = 4 * 16 = 64` and
//! `Transform_Row_Shift = 1`, so the analytic scale is
//! `64 / 2^(1 + 4) = 2`. For the `TX_8X32` / `TX_32X8` pair (round
//! 252) `N_w * N_h = 8 * 32 = 256` and `Transform_Row_Shift = 2`
//! (av1-spec p.408), so the analytic scale is
//! `256 / 2^(2 + 4) = 4`. Same `dqDenom = 4` branch as r244 for the
//! r250 pair (length-64 axis); the r251 pair sits at `dqDenom = 1`
//! per §7.12.3 (the dqDenom rule cites TX_32X32 / TX_16X32 /
//! TX_32X16 / TX_16X64 / TX_64X16 → 2 and the 64-axis sizes → 4 —
//! `TX_4X16` / `TX_16X4` fall to the default `1` arm). The r252
//! pair (`TX_8X32` / `TX_32X8`) is also outside both §7.12.3
//! enumerated sets and thus sits at `dqDenom = 1` — neither shape
//! is in the `{TX_32X32, TX_16X32, TX_32X16, TX_16X64, TX_64X16}`
//! list, nor in the `{TX_64X64, TX_32X64, TX_64X32}` list.
//!
//! With round 252, all 12 §5.11 rectangular `TX_SIZES_ALL` shapes
//! are dispatcher-supported.
//!
//! Within the square sizes, the per-kernel coverage is:
//!
//! | Kernel | TX_4X4 | TX_8X8 | TX_16X16 | TX_32X32 | TX_64X64 |
//! | ------ | ------ | ------ | -------- | -------- | -------- |
//! | DCT    | ✓      | ✓      | ✓        | ✓        | ✓        |
//! | ADST   | ✓      | ✓      | ✓        | —        | —        |
//! | IDTX   | ✓      | ✓      | ✓        | ✓        | —        |
//! | WHT    | ✓      | —      | —        | —        | —        |
//!
//! ADST stops at size 16 because the §7.13.2.9 inverse-ADST
//! dispatcher itself only routes `n in 2..=4` (sizes 4 / 8 / 16);
//! IDTX stops at size 32 because §7.13.2.15 inverse-identity routes
//! `n in 2..=5` (sizes 4 / 8 / 16 / 32). The spec's §6.10.19
//! `tx_type` derivation forces `DCT_DCT` for the unreachable
//! combinations (e.g. `TX_32X32` with non-DCT row kernel is
//! disallowed by the §5.11.48 tx-type-set restrictions).
//!
//! ## FLIPADST handling
//!
//! Per §7.13.3 + §7.12.3 step 3, the FLIPADST family runs the same
//! butterfly schedule as ADST; the **flip** is purely a destination-
//! coordinate transform applied during the decoder's frame-buffer
//! write (`xx = flipLR ? w - j - 1 : j`, `yy = flipUD ? h - i - 1 :
//! i`). The encoder mirror is to flip the spatial residual buffer
//! **before** the forward transform runs and then run the plain
//! ADST kernel on the flipped data.
//!
//! Per-tx-type flip axis:
//!
//! | tx_type             | flip rows (vertical) | flip cols (horizontal) |
//! | ------------------- | -------------------- | ---------------------- |
//! | `FLIPADST_DCT`      | yes                  | no                     |
//! | `DCT_FLIPADST`      | no                   | yes                    |
//! | `FLIPADST_FLIPADST` | yes                  | yes                    |
//! | `ADST_FLIPADST`     | no                   | yes                    |
//! | `FLIPADST_ADST`     | yes                  | no                     |
//! | `V_FLIPADST`        | yes                  | no                     |
//! | `H_FLIPADST`        | no                   | yes                    |
//!
//! "Row kernel = FLIPADST" ⇒ the **horizontal** axis runs ADST on
//! a horizontally-flipped row, i.e. flip columns. "Column kernel =
//! FLIPADST" ⇒ the **vertical** axis runs ADST on a vertically-
//! flipped column, i.e. flip rows. The table above is the union of
//! the two axis flips per tx_type.
//!
//! ## Shift envelope (lossy arm)
//!
//! §7.13.3 lossy decoder pipeline (per row of the row pass):
//!
//! ```text
//!   T[j]            = Dequant[i][j]                       (input)
//!   T[j]            = kernel_row(T[j], log2W, r=BD+8)     (1D pass)
//!   Residual[i][j]  = Round2(T[j], rowShift)              (per-row right-shift)
//! ```
//!
//! and per column:
//!
//! ```text
//!   T[i]            = Residual[i][j]
//!   T[i]            = kernel_col(T[i], log2H, r=Max(BD+6, 16))
//!   Residual[i][j]  = Round2(T[i], colShift=4)            (per-column right-shift)
//! ```
//!
//! Encoder mirror (column-then-row, the transpose of the decoder's
//! row-then-column composition — the decoder's column pass runs
//! against `Residual` last, so the encoder's column pass runs
//! against the spatial residual **first**):
//!
//! ```text
//!   T[i]                = input[i*w + j]                       (residual cell)
//!   T[i]                = fwd_kernel_col(T[i], log2H)          (1D pass)
//!   intermediate[i*w+j] = T[i]
//!   T[j]                = intermediate[i*w + j]
//!   T[j]                = fwd_kernel_row(T[j], log2W)
//!   coeff[i*w + j]      = T[j]
//! ```
//!
//! Note the encoder does **not** apply the `<< rowShift` / `<<
//! colShift` pre-scales that would, in theory, cancel the decoder's
//! `Round2(_, shift)` post-scales bit-exactly. Two reasons:
//!
//!   1. The decoder applies an explicit between-stage `Clip3` at
//!      `colClampRange = Max(BitDepth + 6, 16)` bits. Pre-shifting
//!      the encoder's coefficients by `2^4` or `2^6` for larger
//!      block sizes pushes the inverse pipeline's intermediate
//!      values past the decoder clamp, breaking the round-trip
//!      catastrophically (the saturation truncates the kernel
//!      output rather than preserving it).
//!
//!   2. The spec is set up so that the per-axis kernel matrix has
//!      `M^T · M ≈ N · I` (squared L2 norm per column ≈ `N` for the
//!      sqrt(N)-normalised kernels). The decoder's
//!      `2^(rowShift+colShift)` divisor is chosen to balance this
//!      `N^2 ≈` kernel × kernel gain back toward unity — for the
//!      five square sizes the inverse-only effective gain per cell
//!      is `N^2 / 2^(rowShift + colShift)` = `{1, 2, 4, 16, 64}`
//!      for `TX_{4..64}X{4..64}` respectively. A round-trip with
//!      no encoder pre-shift then has per-cell gain
//!      `N^2 / 2^(rowShift + colShift)` (same number).
//!
//! The §7.13.3 inverse pipeline's combined per-cell gain on the
//! lossy arm is therefore the round-trip scale this dispatcher
//! produces. A real encoder driver that pairs this dispatcher with
//! [`super::forward_quantize`] doesn't see the gain — the quantizer
//! divides by the per-stage gain to recover bit-correct coefficients
//! against the spec's per-tx-size quantizer step.
//!
//! ## Lossless arm
//!
//! §7.13.3 routes the `Lossless = 1` path through a `TX_4X4`-only
//! WHT pipeline:
//!
//! ```text
//!   row pass: inverse_wht4(T, shift = 2)
//!   col pass: inverse_wht4(T, shift = 0)
//! ```
//!
//! The forward [`super::forward_wht::forward_wht_4x4`] already
//! implements the **bit-exact** inverse of this pipeline (column
//! pass with `shift = 0` first, then row pass with `shift = 2` —
//! pre-multiplied by `<< shift` to cancel the inverse's `>> shift`).
//! This dispatcher delegates `lossless == true` to that primitive.
//!
//! ## Round-trip behaviour
//!
//! For the **lossless** arm (`TX_4X4` only, any input residual):
//! `inverse_transform_2d(forward_transform_2d(x), TX_4X4, _, true)
//! == x` exactly. The WHT integer butterflies + pre-shift envelope
//! preserve every input bit.
//!
//! For the **lossy** arms (all other tx_size × tx_type combinations
//! with a sufficiently small residual magnitude that the inverse-
//! side between-stage `Clip3` doesn't saturate): `inverse_transform_2d
//! (forward_transform_2d(x), tx_size, tx_type, false) ≈ scale * x`
//! with a small per-cell rounding error. The `scale` factor per
//! tx-size for the DCT family is `N^2 / 2^(rowShift + colShift)` —
//! evaluated:
//!
//! | tx_size    | N  | rowShift | colShift | scale |
//! | ---------- | -- | -------- | -------- | ----- |
//! | `TX_4X4`   | 4  | 0        | 4        | 1     |
//! | `TX_8X8`   | 8  | 1        | 4        | 2     |
//! | `TX_16X16` | 16 | 2        | 4        | 4     |
//! | `TX_32X32` | 32 | 2        | 4        | 16    |
//! | `TX_64X64` | 64 | 2        | 4        | 64    |
//!
//! For ADST × ADST the per-axis matrix norm is the same shape so the
//! round-trip scale is identical to the DCT scale at each size. For
//! IDTX × IDTX the per-axis scale is the inverse-identity multiplier
//! from §7.13.2.11..§7.13.2.14 — n = 3 / n = 5 are exact integer
//! multiplies (`× 2`, `× 4`), while n = 2 / n = 4 have small Round2-
//! per-cell floor error. The roundtrip tests in this module use a
//! conservative per-cell `max_abs_error` bound; the typical worst
//! case is a few LSBs per cell from stacked `Round2(_, 12)`
//! operations.
//!
//! For TX_32X32 / TX_64X64 the inverse-side between-stage clamp
//! (`colClampRange = Max(BitDepth + 6, 16) = 16` bits at `BitDepth
//! = 8`, i.e. `±32768`) bounds the intermediate magnitude after the
//! inverse row pass to `±32768`. For arbitrary `±128` residual
//! inputs this clamp can trip on the larger sizes; the round-trip
//! tests therefore restrict TX_32X32 / TX_64X64 input residuals to
//! a magnitude small enough that the inverse pipeline doesn't
//! saturate (`±4` for TX_64X64, `±16` for TX_32X32). The smaller
//! sizes can use the full `±128` range.
//!
//! ## Rectangular scaling (rounds 235 + 238 + 241 + 244 + 250)
//!
//! For the `|log2W - log2H| == 1` sizes the §7.13.3 decoder inserts
//! a per-row `Round2(T[j] * 2896, 12)` step BEFORE the row kernel
//! runs. The encoder mirror — time-reversed of the decoder pipeline
//! — places the SAME `Round2(T[j] * 2896, 12)` step per row AFTER
//! the row kernel runs (the encoder's last pass). The constant
//! `2896` is exactly `Cos128_Lookup[32]` (`round(4096 * cos(pi/4)) =
//! round(4096 / sqrt(2))`), so the per-pair contribution is
//! `(2896 / 4096)^2 ≈ 1/2`. Combined with the per-axis kernel norm
//! and the §7.13.3 row + col shift envelope, the empirical
//! round-trip per-cell scale on a constant-DC input evaluates to
//!
//! | tx_size    | empirical per-cell round-trip scale |
//! | ---------- | ----------------------------------- |
//! | `TX_4X8`   | `1/4`                               |
//! | `TX_8X4`   | `1/4`                               |
//! | `TX_8X16`  | `1/2`                               |
//! | `TX_16X8`  | `1/2`                               |
//! | `TX_16X32` | `2`                                 |
//! | `TX_32X16` | `2`                                 |
//! | `TX_32X64` | `8`                                 |
//! | `TX_64X32` | `8`                                 |
//! | `TX_16X64` | `4`                                 |
//! | `TX_64X16` | `4`                                 |
//! | `TX_4X16`  | `1/2`                               |
//! | `TX_16X4`  | `1/2`                               |
//! | `TX_8X32`  | `1`                                 |
//! | `TX_32X8`  | `1`                                 |
//!
//! TX_4X8 / TX_8X4 inherit the TX_4X4 `1/4` scale (short-side-4
//! kernel norm dominates the row/col composition); TX_8X16 /
//! TX_16X8 land at `1/2` — one step up — because the larger
//! `N_w * N_h = 8 * 16 = 128` kernel norm gains `4×` over the
//! short-side-4 pair's `N_w * N_h = 4 * 8 = 32`, while the larger
//! row-shift envelope (`Transform_Row_Shift[TX_8X16] = 1` vs
//! `0` for TX_4X8) eats back `2×`, netting a `2×` larger
//! round-trip per-cell scale. TX_16X32 / TX_32X16 land at `2` —
//! two steps up from TX_8X16 — because `N_w * N_h = 16 * 32 = 512`
//! gains another `4×` over the short-side-8 pair while the
//! row-shift envelope stays at `Transform_Row_Shift = 1` (same as
//! TX_8X16 / TX_16X8), so the full `4×` lands in the per-cell
//! round-trip scale. TX_32X64 / TX_64X32 land at `8` — two more
//! steps up — because `N_w * N_h = 32 * 64 = 2048` gains another
//! `4×` over the short-side-16 pair while the row-shift envelope
//! stays at `Transform_Row_Shift = 1`.
//!
//! Walking the analytic derivation for TX_32X64: round-trip
//! per-cell gain = `(N_w * N_h) / 2^(row_shift + col_shift) *
//! (2896 / 4096)^2` = `2048 / 2^(1 + 4) * 1/2` = `2048 / 32 * 1/2` =
//! `32`. The factor of `4` between the analytic `32` and empirical
//! `8` matches the `4×` ratio already documented across the
//! `|log2W - log2H| == 1` family — i.e. the constant-DC probe's
//! input-mass factor on the inverse-side after the row-pass is
//! invariant from TX_8X16 forward.
//!
//! Walking the analytic derivation for TX_16X64 (round 250, the
//! first `|log2W - log2H| == 2` shape): round-trip per-cell gain =
//! `(N_w * N_h) / 2^(row_shift + col_shift)` with NO rectangular
//! `(2896 / 4096)^2` factor (the `× 2896` step gates on
//! `|log2W - log2H| == 1` only). `(16 * 64) / 2^(2 + 4)` =
//! `1024 / 64` = `16`. The factor of `4` between the analytic `16`
//! and empirical `4` matches the same `4×` constant-DC input-mass
//! factor documented across the rectangular family. TX_64X16 is
//! the transpose and produces the same per-cell scale.
//!
//! The quantizer in [`super::forward_quantize`] absorbs this fixed
//! per-tx-size gain through the `dqDenom` factor (the same path
//! used for the square sizes): `dequant_denom(TX_16X32) =
//! dequant_denom(TX_32X16) = 2` per §7.12.3 (any 32-axis transform
//! promotes `dqDenom` to `2`), and `dequant_denom(TX_32X64) =
//! dequant_denom(TX_64X32) = dequant_denom(TX_16X64) =
//! dequant_denom(TX_64X16) = 4` (any 64-axis transform promotes
//! `dqDenom` to `4`). TX_8X16 / TX_16X8 still sit at `dqDenom = 1`.
//! A downstream encoder pipeline that pairs this dispatcher with
//! the existing quantizer sees a bit-correct coefficient stream
//! against the spec's per-tx-size quantizer step.

use crate::cdf::{
    ADST_ADST, ADST_DCT, ADST_FLIPADST, DCT_ADST, DCT_DCT, DCT_FLIPADST, FLIPADST_ADST,
    FLIPADST_DCT, FLIPADST_FLIPADST, H_ADST, H_DCT, H_FLIPADST, IDTX, TX_HEIGHT, TX_SIZES_ALL,
    TX_WIDTH, TX_WIDTH_LOG2, V_ADST, V_DCT, V_FLIPADST,
};

use super::forward_adst::{forward_adst_16, forward_adst_4, forward_adst_8};
use super::forward_identity::{forward_idtx_16, forward_idtx_32, forward_idtx_4, forward_idtx_8};
use super::forward_transform::{
    forward_dct_16, forward_dct_32, forward_dct_4, forward_dct_64, forward_dct_8,
};
use super::forward_wht::forward_wht_4x4;

/// §7.13.3 row-pass forward kernel selector — the encoder mirror of
/// [`crate::transform::apply_row_kernel`]. DCT for
/// `{ DCT_DCT, ADST_DCT, FLIPADST_DCT, H_DCT }`; ADST for
/// `{ DCT_ADST, ADST_ADST, DCT_FLIPADST, FLIPADST_FLIPADST,
/// ADST_FLIPADST, FLIPADST_ADST, H_ADST, H_FLIPADST }`; identity
/// for `{ IDTX, V_DCT, V_ADST, V_FLIPADST }`.
fn forward_row_kernel(t: &mut [i64], tx_type: usize, log2_w: u32) {
    let r = 32; // signature parity only; forward kernels ignore r.
    if matches!(tx_type, x if x == DCT_DCT || x == ADST_DCT || x == FLIPADST_DCT || x == H_DCT) {
        forward_dct_dispatch(t, log2_w, r);
    } else if matches!(
        tx_type,
        x if x == DCT_ADST
            || x == ADST_ADST
            || x == DCT_FLIPADST
            || x == FLIPADST_FLIPADST
            || x == ADST_FLIPADST
            || x == FLIPADST_ADST
            || x == H_ADST
            || x == H_FLIPADST
    ) {
        forward_adst_dispatch(t, log2_w, r);
    } else {
        debug_assert!(
            matches!(tx_type, x if x == IDTX || x == V_DCT || x == V_ADST || x == V_FLIPADST)
        );
        forward_idtx_dispatch(t, log2_w);
    }
}

/// §7.13.3 column-pass forward kernel selector — the encoder mirror
/// of [`crate::transform::apply_col_kernel`]. DCT for
/// `{ DCT_DCT, DCT_ADST, DCT_FLIPADST, V_DCT }`; ADST for
/// `{ ADST_DCT, ADST_ADST, FLIPADST_DCT, FLIPADST_FLIPADST,
/// ADST_FLIPADST, FLIPADST_ADST, V_ADST, V_FLIPADST }`; identity
/// for `{ IDTX, H_DCT, H_ADST, H_FLIPADST }`.
fn forward_col_kernel(t: &mut [i64], tx_type: usize, log2_h: u32) {
    let r = 32;
    if matches!(tx_type, x if x == DCT_DCT || x == DCT_ADST || x == DCT_FLIPADST || x == V_DCT) {
        forward_dct_dispatch(t, log2_h, r);
    } else if matches!(
        tx_type,
        x if x == ADST_DCT
            || x == ADST_ADST
            || x == FLIPADST_DCT
            || x == FLIPADST_FLIPADST
            || x == ADST_FLIPADST
            || x == FLIPADST_ADST
            || x == V_ADST
            || x == V_FLIPADST
    ) {
        forward_adst_dispatch(t, log2_h, r);
    } else {
        debug_assert!(
            matches!(tx_type, x if x == IDTX || x == H_DCT || x == H_ADST || x == H_FLIPADST)
        );
        forward_idtx_dispatch(t, log2_h);
    }
}

fn forward_dct_dispatch(t: &mut [i64], n: u32, r: u32) {
    match n {
        2 => forward_dct_4(t, r),
        3 => forward_dct_8(t, r),
        4 => forward_dct_16(t, r),
        5 => forward_dct_32(t, r),
        6 => forward_dct_64(t, r),
        _ => panic!("oxideav-av1 forward_dct_dispatch: n must be 2..=6, got {n}"),
    }
}

fn forward_adst_dispatch(t: &mut [i64], n: u32, r: u32) {
    match n {
        2 => forward_adst_4(t, r),
        3 => forward_adst_8(t, r),
        4 => forward_adst_16(t, r),
        _ => panic!(
            "oxideav-av1 forward_adst_dispatch: ADST is only defined for n in 2..=4 \
             (the §7.13.2.9 inverse-ADST dispatcher's range), got {n} — \
             §6.10.19 tx_type derivation forces DCT_DCT outside this range",
        ),
    }
}

fn forward_idtx_dispatch(t: &mut [i64], n: u32) {
    match n {
        2 => forward_idtx_4(t),
        3 => forward_idtx_8(t),
        4 => forward_idtx_16(t),
        5 => forward_idtx_32(t),
        _ => panic!(
            "oxideav-av1 forward_idtx_dispatch: IDTX is only defined for n in 2..=5 \
             (the §7.13.2.15 inverse-identity dispatcher's range), got {n} — \
             §6.10.19 tx_type derivation forces DCT_DCT outside this range",
        ),
    }
}

/// Per-tx-type flip-axis decoder. Returns `(flip_rows, flip_cols)`
/// for the §7.12.3 step-3 frame-buffer flip the decoder applies
/// post-inverse — the encoder must apply the same flip on the
/// spatial residual before the forward transform runs.
fn flip_axes(tx_type: usize) -> (bool, bool) {
    // Vertical-axis flip (= flip rows in the row-major layout): the
    // *column* kernel is FLIPADST.
    let flip_rows = matches!(
        tx_type,
        x if x == FLIPADST_DCT
            || x == FLIPADST_FLIPADST
            || x == FLIPADST_ADST
            || x == V_FLIPADST
    );
    // Horizontal-axis flip (= flip cols in the row-major layout):
    // the *row* kernel is FLIPADST.
    let flip_cols = matches!(
        tx_type,
        x if x == DCT_FLIPADST
            || x == FLIPADST_FLIPADST
            || x == ADST_FLIPADST
            || x == H_FLIPADST
    );
    (flip_rows, flip_cols)
}

fn apply_flip(input: &[i64], w: usize, h: usize, flip_rows: bool, flip_cols: bool) -> Vec<i64> {
    if !flip_rows && !flip_cols {
        return input.to_vec();
    }
    let mut out = vec![0i64; w * h];
    for i in 0..h {
        let src_i = if flip_rows { h - 1 - i } else { i };
        for j in 0..w {
            let src_j = if flip_cols { w - 1 - j } else { j };
            out[i * w + j] = input[src_i * w + src_j];
        }
    }
    out
}

/// §7.13.3-equivalent forward 2D transform dispatcher — the encoder
/// counterpart of [`crate::transform::inverse_transform_2d`].
///
/// Consumes `input` (a row-major spatial residual buffer of length
/// `w * h` where `w = Tx_Width[tx_size]`, `h = Tx_Height[tx_size]`)
/// and returns the row-major coefficient buffer of the same length.
/// The forward composition is **column pass first, then row pass**
/// (the transpose of the decoder's row-then-column composition).
///
/// `tx_size` must be one of the five square sizes
/// ([`crate::cdf::TX_4X4`] / [`crate::cdf::TX_8X8`] /
/// [`crate::cdf::TX_16X16`] / [`crate::cdf::TX_32X32`] /
/// [`crate::cdf::TX_64X64`]) **or** one of the eight
/// `|log2W - log2H| == 1` rectangular pairs at the short-side-4,
/// short-side-8, short-side-16 and short-side-32 sizes —
/// [`crate::cdf::TX_4X8`] / [`crate::cdf::TX_8X4`] (round 235),
/// [`crate::cdf::TX_8X16`] / [`crate::cdf::TX_16X8`] (round 238),
/// [`crate::cdf::TX_16X32`] / [`crate::cdf::TX_32X16`] (round 241),
/// and [`crate::cdf::TX_32X64`] / [`crate::cdf::TX_64X32`] (round
/// 244) **or** all three `|log2W - log2H| == 2` pairs —
/// [`crate::cdf::TX_16X64`] / [`crate::cdf::TX_64X16`] (round 250 —
/// the long-side-64 endpoint of the family),
/// [`crate::cdf::TX_4X16`] / [`crate::cdf::TX_16X4`] (round 251 —
/// the short-side-4 endpoint), and
/// [`crate::cdf::TX_8X32`] / [`crate::cdf::TX_32X8`] (round 252 —
/// the mid pair at `min_log2 == 3` / `max_log2 == 5`). All twelve
/// `TX_SIZES_ALL` rectangular shapes are now in arc.
///
/// `plane_tx_type` must be one of the 16 §6.10.19 ordinals; the per-
/// tx-type / per-tx-size kernel coverage is the intersection of the
/// per-axis [`forward_dct_dispatch`] / [`forward_adst_dispatch`] /
/// [`forward_idtx_dispatch`] ranges (see the module docs' coverage
/// table). FLIPADST flips the spatial residual on the appropriate
/// axis before the plain ADST kernel runs (per §7.13.3 + §7.12.3
/// step 3).
///
/// `lossless` is the per-block §6.8.11 `Lossless` flag. When `true`,
/// `tx_size` must be `TX_4X4` and the dispatcher routes through the
/// bit-exact WHT path in [`super::forward_wht::forward_wht_4x4`].
///
/// # Panics
///
/// * `tx_size >= TX_SIZES_ALL`.
/// * `tx_size` is a rectangular family currently out of arc scope
///   (any rectangular size beyond `TX_4X8` / `TX_8X4` / `TX_8X16` /
///   `TX_16X8` / `TX_16X32` / `TX_32X16` / `TX_32X64` / `TX_64X32` /
///   `TX_16X64` / `TX_64X16` / `TX_4X16` / `TX_16X4` /
///   `TX_8X32` / `TX_32X8`).
/// * `input.len() != w * h`.
/// * `lossless == true` with `tx_size != TX_4X4`.
/// * `(tx_size, plane_tx_type)` selects an out-of-range kernel size
///   (e.g. ADST at `TX_32X32`).
pub fn forward_transform_2d(
    input: &[i64],
    tx_size: usize,
    plane_tx_type: usize,
    lossless: bool,
) -> Vec<i64> {
    assert!(
        tx_size < TX_SIZES_ALL,
        "oxideav-av1 forward_transform_2d: tx_size {tx_size} out of range (TX_SIZES_ALL = {TX_SIZES_ALL})"
    );
    let log2_w = TX_WIDTH_LOG2[tx_size] as u32;
    let log2_h = (TX_HEIGHT[tx_size] as u32).trailing_zeros();
    let w = TX_WIDTH[tx_size];
    let h = TX_HEIGHT[tx_size];
    // Square (|log2W - log2H| == 0) is universally supported. The
    // |log2W - log2H| == 1 short-side-4 pair (TX_4X8 / TX_8X4) was
    // landed in round 235; round 238 extended the rectangular arc
    // by the next |log2W - log2H| == 1 pair at short-side-8 —
    // TX_8X16 and TX_16X8. Round 241 extended the arc to the
    // short-side-16 pair (TX_16X32 / TX_32X16). Round 244 closed the
    // |log2W - log2H| == 1 chain with the largest-kernel pair —
    // short-side-32 (TX_32X64 / TX_64X32). Round 250 opened the
    // |log2W - log2H| == 2 family at the long-side-64 endpoint —
    // TX_16X64 / TX_64X16 (min_log2 == 4 / max_log2 == 6). Round
    // 251 extended the |log2W - log2H| == 2 family at the opposite
    // endpoint — TX_4X16 / TX_16X4 at min_log2 == 2 / max_log2 == 4
    // (short-side-4 / long-side-16). Round 252 closes the
    // |log2W - log2H| == 2 family with the mid pair —
    // TX_8X32 / TX_32X8 at min_log2 == 3 / max_log2 == 5
    // (short-side-8 / long-side-32). All twelve rectangular
    // §5.11 TX_SIZES_ALL shapes are now in arc.
    let is_square = log2_w == log2_h;
    let min_log2 = core::cmp::min(log2_w, log2_h);
    let max_log2 = core::cmp::max(log2_w, log2_h);
    let diff = log2_w.abs_diff(log2_h);
    let is_supported_rect_diff1 = diff == 1
        && (min_log2 == 2 || min_log2 == 3 || min_log2 == 4 || min_log2 == 5)
        && max_log2 == min_log2 + 1;
    // The |log2W - log2H| == 2 arm: round 250 lands min_log2 == 4 /
    // max_log2 == 6 (TX_16X64 / TX_64X16); round 251 adds min_log2
    // == 2 / max_log2 == 4 (TX_4X16 / TX_16X4); round 252 adds the
    // mid pair min_log2 == 3 / max_log2 == 5 (TX_8X32 / TX_32X8).
    let is_supported_rect_diff2 = diff == 2
        && ((min_log2 == 2 && max_log2 == 4)
            || (min_log2 == 3 && max_log2 == 5)
            || (min_log2 == 4 && max_log2 == 6));
    let is_supported_rect = is_supported_rect_diff1 || is_supported_rect_diff2;
    assert!(
        is_square || is_supported_rect,
        "oxideav-av1 forward_transform_2d: rectangular tx_size {tx_size} (w={w}, h={h}) \
         not supported in this arc — supported shapes are the five square sizes \
         (TX_4X4 / TX_8X8 / TX_16X16 / TX_32X32 / TX_64X64), the eight \
         |log2W - log2H| == 1 pairs (TX_4X8 / TX_8X4 / TX_8X16 / TX_16X8 / \
         TX_16X32 / TX_32X16 / TX_32X64 / TX_64X32), and three \
         |log2W - log2H| == 2 pairs (TX_16X64 / TX_64X16, \
         TX_4X16 / TX_16X4, and TX_8X32 / TX_32X8)"
    );
    assert_eq!(
        input.len(),
        w * h,
        "oxideav-av1 forward_transform_2d: input length {} != w*h = {}",
        input.len(),
        w * h
    );

    if lossless {
        assert_eq!(
            w, 4,
            "oxideav-av1 forward_transform_2d: lossless arm requires tx_size = TX_4X4"
        );
        // The WHT path ignores plane_tx_type — same as the inverse.
        return forward_wht_4x4(input).to_vec();
    }

    // Apply the §7.12.3 step-3 flip on the spatial residual before
    // the forward transform runs — encoder mirror of the decoder's
    // post-inverse frame-buffer flip.
    let (flip_rows, flip_cols) = flip_axes(plane_tx_type);
    let flipped = apply_flip(input, w, h, flip_rows, flip_cols);

    let mut work = flipped;

    // Column pass first (the encoder's first pass = the decoder's
    // last pass).
    let mut col_buf = vec![0i64; h];
    for j in 0..w {
        for i in 0..h {
            col_buf[i] = work[i * w + j];
        }
        forward_col_kernel(&mut col_buf, plane_tx_type, log2_h);
        for i in 0..h {
            work[i * w + j] = col_buf[i];
        }
    }

    // Row pass.
    let mut row_buf = vec![0i64; w];
    for i in 0..h {
        row_buf.copy_from_slice(&work[i * w..(i + 1) * w]);
        forward_row_kernel(&mut row_buf, plane_tx_type, log2_w);
        // §7.13.3 rectangular scaling — encoder mirror. The decoder
        // applies `Round2(T[j] * 2896, 12)` per row BEFORE the row
        // kernel runs whenever |log2W - log2H| == 1 (the gate is
        // EXACTLY `Abs(log2W - log2H) == 1` per av1-spec p.305 line
        // "if (Abs(log2W - log2H) == 1) T[j] = Round2(T[j] * 2896,
        // 12)" — the |log2W - log2H| == 2 family stays at the bare
        // per-axis kernel composition). The encoder is the time-
        // reverse of that pipeline, so the same scale appears AFTER
        // the row kernel (the encoder's last pass). Both sides thus
        // contribute one factor of 2896 / 4096; the net rectangular
        // gain is (2896 / 4096)^2 ≈ 1/2 per rectangular axis pair —
        // the same multiplier whether the pair is TX_4X8 / TX_8X4,
        // TX_8X16 / TX_16X8, TX_16X32 / TX_32X16, or TX_32X64 /
        // TX_64X32. The TX_16X64 / TX_64X16 pair (round 250),
        // TX_4X16 / TX_16X4 pair (round 251), and TX_8X32 / TX_32X8
        // pair (round 252) are all |log2W - log2H| == 2 shapes and
        // skip this gate entirely.
        // The per-tx-size differentiator is the §7.13.3
        // row + col shift envelope (and, for 32/64-axis transforms,
        // the §7.12.3 dqDenom applied by [`super::forward_quantize`]).
        if log2_w.abs_diff(log2_h) == 1 {
            for slot in row_buf.iter_mut() {
                *slot = round2_12(*slot * 2896);
            }
        }
        work[i * w..(i + 1) * w].copy_from_slice(&row_buf);
    }

    // r460 — the per-size GAIN correction. Measured against the
    // decoder's own §7.13.3 chain (forward → quantise at a fine step →
    // §7.12.3 dequantise → inverse), the column-then-row composition
    // above reproduced the residual at 1/4 amplitude for the 4×4
    // family, 1/2 for 8×8 / 8×16 / 4×16, exactly 1 only at 16×16 /
    // 8×32, 2 at 16×32 and 4 at 32×32 — every TU size but 16×16 was
    // coding a residual the decoder rebuilds at the wrong amplitude,
    // and the RD search had been compensating around it (pre-r460
    // KEY frames ran ~2× the bytes of a third-party encoder at equal
    // PSNR). The kernels' own gain is 2^(floor((log2W + log2H) / 2) - 4)
    // up to a 16×16-sized product and 2^(log2W + log2H - 8) beyond it
    // relative to what the §7.13.3 shift envelope (`Transform_Row_Shift`
    // + `colShift = 4`) undoes; scale the coefficients by the inverse
    // so `inverse(quantise(forward(r))) ≈ r` at every size.
    let sum = log2_w + log2_h;
    let shift: i32 = if sum <= 8 {
        4 - (sum / 2) as i32
    } else {
        8 - sum as i32
    };
    if shift > 0 {
        for v in work.iter_mut() {
            *v <<= shift;
        }
    } else if shift < 0 {
        let s = (-shift) as u32;
        let half = 1i64 << (s - 1);
        for v in work.iter_mut() {
            *v = (*v + half) >> s;
        }
    }

    work
}

/// `Round2(x, 12)` — the §4.7.2 round-to-nearest-with-ties-up shift
/// the §7.13.3 rectangular scale step uses.
#[inline]
fn round2_12(x: i64) -> i64 {
    (x + (1i64 << 11)) >> 12
}

#[cfg(test)]
mod tests {
    use super::*;
    use crate::cdf::{
        TX_16X16, TX_16X32, TX_16X4, TX_16X64, TX_16X8, TX_32X16, TX_32X32, TX_32X64, TX_32X8,
        TX_4X16, TX_4X4, TX_4X8, TX_64X16, TX_64X32, TX_64X64, TX_8X16, TX_8X32, TX_8X4, TX_8X8,
    };
    use crate::transform::inverse_transform_2d;

    // -------------------------------------------------------------
    // Lossless arm — bit-exact roundtrip.
    // -------------------------------------------------------------

    #[test]
    fn lossless_tx_4x4_zero_input_yields_zero() {
        let input = vec![0i64; 16];
        let coeffs = forward_transform_2d(&input, TX_4X4, DCT_DCT, true);
        assert_eq!(coeffs.len(), 16);
        for &v in coeffs.iter() {
            assert_eq!(v, 0);
        }
    }

    #[test]
    fn lossless_tx_4x4_bit_exact_roundtrip() {
        // Pseudo-random pixel residuals in [-128, 127]. The WHT
        // chain is a pure integer butterfly ⇒ round-trip is
        // bit-exact regardless of input.
        let mut input = vec![0i64; 16];
        let mut s: u64 = 0xCAFE_F00D_DEAD_BEEF;
        for v in input.iter_mut() {
            s = s
                .wrapping_mul(6364136223846793005)
                .wrapping_add(1442695040888963407);
            *v = ((s >> 32) & 0xFF) as i64 - 128;
        }
        let coeffs = forward_transform_2d(&input, TX_4X4, DCT_DCT, true);
        let recovered = inverse_transform_2d(&coeffs, TX_4X4, DCT_DCT, 8, true);
        assert_eq!(
            recovered, input,
            "lossless WHT round-trip diverged from bit-exact"
        );
    }

    #[test]
    fn lossless_tx_4x4_bit_exact_roundtrip_extreme_values() {
        // The §7.13.3 note bounds Residual to `1 + BitDepth` bits =
        // 9 bits for BD = 8 (i.e. `[-256, 255]`).
        for &v in [-256i64, -255, -128, -1, 0, 1, 127, 128, 255].iter() {
            let input = vec![v; 16];
            let coeffs = forward_transform_2d(&input, TX_4X4, DCT_DCT, true);
            let recovered = inverse_transform_2d(&coeffs, TX_4X4, DCT_DCT, 8, true);
            assert_eq!(recovered, input, "lossless extreme-value {v} diverged");
        }
    }

    // -------------------------------------------------------------
    // Lossy arm — roundtrip with bounded error per kernel.
    //
    // The expected per-cell round-trip scale = (kernel × kernel
    // squared-norm gain N per axis) / (2^(rowShift + colShift)
    // decoder post-shifts):
    //
    //   * DCT × DCT  TX_NxN:  N^2 / 2^(rowShift + colShift)
    //     evaluated per size:
    //         TX_4X4   : 16 /  16 = 1
    //         TX_8X8   : 64 /  32 = 2
    //         TX_16X16 : 256/  64 = 4
    //         TX_32X32 : 1024/ 64 = 16
    //         TX_64X64 : 4096/ 64 = 64
    //   * ADST × ADST: same per-axis matrix norm as DCT (kernel is
    //     sqrt(N)-normalised) ⇒ same per-size scale as the DCT table.
    //   * IDTX × IDTX: per-axis scalar (c / 4096)^2 with
    //     c ∈ {5793, 8192, 11586, 16384} for N ∈ {4, 8, 16, 32}.
    //     Two-axis scale then divided by 2^(rowShift + colShift).
    //
    // For TX_32X32 / TX_64X64 the inverse-side between-stage
    // `Clip3` at 16 bits saturates intermediate kernel outputs
    // exceeding ±32768. The roundtrip tests scale down the input
    // residual magnitudes accordingly: ±16 for TX_32X32, ±4 for
    // TX_64X64 keeps the inverse pipeline within the clamp.
    // -------------------------------------------------------------

    fn lcg_residual(seed: u64, n: usize) -> Vec<i64> {
        lcg_residual_bound(seed, n, 128)
    }

    fn lcg_residual_bound(seed: u64, n: usize, bound: i64) -> Vec<i64> {
        let mut out = vec![0i64; n];
        let mut s = seed;
        for v in out.iter_mut() {
            s = s
                .wrapping_mul(6364136223846793005)
                .wrapping_add(1442695040888963407);
            let lim = 2 * bound + 1;
            let raw = ((s >> 32) & 0xFFFF) as i64;
            *v = (raw % lim) - bound;
        }
        out
    }

    /// Roundtrip checker. Verifies `inverse(forward(x)) * denom ≈
    /// scale_num * x` per cell within `max_err * denom`. Allows
    /// fractional per-cell scales (e.g. `scale_num = 1, denom = 4`
    /// for the TX_4X4 empirical 0.25 round-trip).
    fn check_roundtrip_frac(
        input: &[i64],
        tx_size: usize,
        tx_type: usize,
        scale_num: i64,
        denom: i64,
        max_err: i64,
    ) {
        let coeffs = forward_transform_2d(input, tx_size, tx_type, false);
        let recovered = inverse_transform_2d(&coeffs, tx_size, tx_type, 8, false);
        assert_eq!(recovered.len(), input.len());
        for (i, (&got, &orig)) in recovered.iter().zip(input.iter()).enumerate() {
            let lhs = got * denom;
            let rhs = scale_num * orig;
            let err = (lhs - rhs).abs();
            let bound = max_err * denom;
            assert!(
                err <= bound,
                "tx_size={tx_size} tx_type={tx_type} cell {i}: orig={orig}, \
                 got={got}, expected≈{scale_num}/{denom}*orig={}, |err*{denom}|={err} > bound {bound}",
                rhs as f64 / denom as f64,
            );
        }
    }

    /// Roundtrip checker for FLIPADST family — compares the
    /// recovered residual against the **flipped** input (per
    /// `flip_axes`), since `inverse_transform_2d` does not itself
    /// apply the §7.12.3 step-3 flip (that runs externally on the
    /// frame-buffer write).
    fn check_roundtrip_flip(
        input: &[i64],
        tx_size: usize,
        tx_type: usize,
        scale_num: i64,
        denom: i64,
        max_err: i64,
    ) {
        let coeffs = forward_transform_2d(input, tx_size, tx_type, false);
        let recovered = inverse_transform_2d(&coeffs, tx_size, tx_type, 8, false);
        let w = TX_WIDTH[tx_size];
        let h = TX_HEIGHT[tx_size];
        let (flip_rows, flip_cols) = flip_axes(tx_type);
        let flipped = apply_flip(input, w, h, flip_rows, flip_cols);
        assert_eq!(recovered.len(), flipped.len());
        for (i, (&got, &orig)) in recovered.iter().zip(flipped.iter()).enumerate() {
            let lhs = got * denom;
            let rhs = scale_num * orig;
            let err = (lhs - rhs).abs();
            let bound = max_err * denom;
            assert!(
                err <= bound,
                "tx_size={tx_size} tx_type={tx_type} cell {i}: \
                 flipped_orig={orig}, got={got}, \
                 expected≈{scale_num}/{denom}*flipped_orig={}, |err*{denom}|={err} > bound {bound}",
                rhs as f64 / denom as f64,
            );
        }
    }

    // DCT_DCT across square sizes. r460 — unit per-cell gain at every
    // size: `inverse(forward(x)) ≈ x` through the decoder's own
    // §7.13.3 chain (the pre-r460 kernels reproduced {1/4, 1/2, 1, 4,
    // 4} of the input at {4×4, 8×8, 16×16, 32×32, 64×64}; the
    // per-size correction at the end of `forward_transform_2d`
    // restores unity).

    #[test]
    fn dct_dct_tx_4x4_roundtrip() {
        let input = lcg_residual(0x1111_2222_3333_4444, 16);
        check_roundtrip_frac(&input, TX_4X4, DCT_DCT, 1, 1, 16);
    }

    #[test]
    fn dct_dct_tx_8x8_roundtrip() {
        let input = lcg_residual(0x5555_6666_7777_8888, 64);
        check_roundtrip_frac(&input, TX_8X8, DCT_DCT, 1, 1, 16);
    }

    #[test]
    fn dct_dct_tx_16x16_roundtrip() {
        // Scale down inputs so the post-quantization range fits in
        // the inverse's 16-bit between-stage clamp.
        let input = lcg_residual_bound(0x9999_AAAA_BBBB_CCCC, 256, 32);
        check_roundtrip_frac(&input, TX_16X16, DCT_DCT, 1, 1, 16);
    }

    #[test]
    fn dct_dct_tx_32x32_roundtrip() {
        let input = lcg_residual_bound(0xDEAD_BEEF_F00D_CAFE, 1024, 8);
        check_roundtrip_frac(&input, TX_32X32, DCT_DCT, 1, 1, 16);
    }

    #[test]
    fn dct_dct_tx_64x64_roundtrip() {
        let input = lcg_residual_bound(0x0123_4567_89AB_CDEF, 4096, 2);
        // larger error bound here reflects the deeper DCT-64
        // butterfly schedule's accumulated `Round2(_, 12)` floor
        // (31-step butterfly graph vs ~5 steps for DCT-4).
        check_roundtrip_frac(&input, TX_64X64, DCT_DCT, 1, 1, 64);
    }

    // ADST × ADST and ADST × DCT combinations. Same per-cell scale
    // shape as DCT (the ADST kernel matrix is sqrt(N)-normalised
    // like DCT).

    #[test]
    fn adst_adst_tx_4x4_roundtrip() {
        let input = lcg_residual(0xABCD_EF01_2345_6789, 16);
        check_roundtrip_frac(&input, TX_4X4, ADST_ADST, 1, 1, 16);
    }

    #[test]
    fn adst_dct_tx_8x8_roundtrip() {
        let input = lcg_residual(0xFACE_BEEF_CAFE_BABE, 64);
        check_roundtrip_frac(&input, TX_8X8, ADST_DCT, 1, 1, 16);
    }

    #[test]
    fn dct_adst_tx_8x8_roundtrip() {
        let input = lcg_residual(0xBADD_CAFE_F00D_0BAD, 64);
        check_roundtrip_frac(&input, TX_8X8, DCT_ADST, 1, 1, 16);
    }

    #[test]
    fn adst_adst_tx_16x16_roundtrip() {
        let input = lcg_residual_bound(0x1357_9BDF_2468_ACE0, 256, 32);
        check_roundtrip_frac(&input, TX_16X16, ADST_ADST, 1, 1, 16);
    }

    // FLIPADST family — compare against the **flipped** input, since
    // `inverse_transform_2d` does not apply the §7.12.3 step-3 flip
    // (the flip is the decoder-side frame-buffer write, applied
    // externally).

    #[test]
    fn flipadst_flipadst_tx_4x4_roundtrip() {
        let input = lcg_residual(0x2233_4455_6677_8899, 16);
        check_roundtrip_flip(&input, TX_4X4, FLIPADST_FLIPADST, 1, 1, 16);
    }

    #[test]
    fn flipadst_dct_tx_8x8_roundtrip() {
        let input = lcg_residual(0xAA55_AA55_AA55_AA55, 64);
        check_roundtrip_flip(&input, TX_8X8, FLIPADST_DCT, 1, 1, 16);
    }

    #[test]
    fn dct_flipadst_tx_8x8_roundtrip() {
        let input = lcg_residual(0x55AA_55AA_55AA_55AA, 64);
        check_roundtrip_flip(&input, TX_8X8, DCT_FLIPADST, 1, 1, 16);
    }

    #[test]
    fn adst_flipadst_tx_16x16_roundtrip() {
        let input = lcg_residual_bound(0xC0DE_F00D_BABE_FACE, 256, 32);
        check_roundtrip_flip(&input, TX_16X16, ADST_FLIPADST, 1, 1, 16);
    }

    #[test]
    fn flipadst_adst_tx_16x16_roundtrip() {
        let input = lcg_residual_bound(0xFACE_BABE_F00D_C0DE, 256, 32);
        check_roundtrip_flip(&input, TX_16X16, FLIPADST_ADST, 1, 1, 16);
    }

    // IDTX family. Per-cell behaviour is exact integer multiply for
    // n = 3 / n = 5 axes (× 2 / × 4 per pass — `c = 8192 / 16384 ∝
    // 2^k`). For n = 2 / n = 4 axes (`c = 5793 / 11586`), the
    // per-cell rounding adds a small Round2 floor; the empirical
    // scale is 1/4 / 1 per cell respectively after the §7.13.3
    // shift envelope.

    #[test]
    fn idtx_tx_4x4_roundtrip() {
        let input = lcg_residual(0x1234_5678_9ABC_DEF0, 16);
        // Empirical: per-cell ≈ 1/4 × input.
        check_roundtrip_frac(&input, TX_4X4, IDTX, 1, 1, 16);
    }

    #[test]
    fn idtx_tx_8x8_roundtrip() {
        let input = lcg_residual(0xDEAD_FACE_BEEF_CAFE, 64);
        // Empirical: per-cell ≈ 1/2 × input.
        check_roundtrip_frac(&input, TX_8X8, IDTX, 1, 1, 16);
    }

    #[test]
    fn idtx_tx_16x16_roundtrip() {
        let input = lcg_residual_bound(0xBAD0_F00D_DEAD_BEEF, 256, 32);
        // Empirical: per-cell ≈ 1 × input.
        check_roundtrip_frac(&input, TX_16X16, IDTX, 1, 1, 16);
    }

    #[test]
    fn idtx_tx_32x32_roundtrip() {
        let input = lcg_residual_bound(0xC001_D00D_FACE_FEED, 1024, 8);
        // Empirical: per-cell ≈ 4 × input (and exact since c = 16384
        // is a power of two).
        check_roundtrip_frac(&input, TX_32X32, IDTX, 1, 1, 16);
    }

    // V_/H_ mixed (DCT × identity) combinations. Same per-cell
    // scale shape as the homogeneous family at the same tx_size
    // (the per-axis kernel norms compose multiplicatively).

    #[test]
    fn v_dct_tx_4x4_roundtrip() {
        // V_DCT: column kernel = DCT, row kernel = identity.
        let input = lcg_residual(0x1010_2020_3030_4040, 16);
        check_roundtrip_frac(&input, TX_4X4, V_DCT, 1, 1, 16);
    }

    #[test]
    fn h_dct_tx_8x8_roundtrip() {
        // H_DCT: row kernel = DCT, column kernel = identity.
        let input = lcg_residual(0x5050_6060_7070_8080, 64);
        check_roundtrip_frac(&input, TX_8X8, H_DCT, 1, 1, 16);
    }

    // -------------------------------------------------------------
    // Rectangular family — `|log2W - log2H| == 1` short-side-4 pair
    // (TX_4X8 / TX_8X4). The empirical round-trip per-cell scale on
    // a constant-DC input is `1/4` — the same shape as the TX_4X4
    // square case, reflecting the per-axis short-side (length-4)
    // kernel norm composition under the §7.13.3 row/col shift
    // envelope and the `(2896 / 4096)^2` rectangular contribution
    // from the encoder + decoder applying the rectangular scale
    // each.
    // -------------------------------------------------------------

    #[test]
    fn rect_tx_4x8_zero_input_yields_zero() {
        let input = vec![0i64; 32];
        let coeffs = forward_transform_2d(&input, TX_4X8, DCT_DCT, false);
        assert_eq!(coeffs.len(), 32);
        for (i, &v) in coeffs.iter().enumerate() {
            assert_eq!(v, 0, "TX_4X8 cell {i}: zero in ⇒ zero out, got {v}");
        }
    }

    #[test]
    fn rect_tx_8x4_zero_input_yields_zero() {
        let input = vec![0i64; 32];
        let coeffs = forward_transform_2d(&input, TX_8X4, DCT_DCT, false);
        assert_eq!(coeffs.len(), 32);
        for (i, &v) in coeffs.iter().enumerate() {
            assert_eq!(v, 0, "TX_8X4 cell {i}: zero in ⇒ zero out, got {v}");
        }
    }

    #[test]
    fn rect_tx_4x8_dct_dct_roundtrip() {
        // TX_4X8 round-trip per-cell scale ≈ 1/4 of input (empirical
        // on a constant-DC probe: input = 64 ⇒ recovered = 16).
        let input = lcg_residual(0x1357_2468_ACE0_BDF1, 32);
        check_roundtrip_frac(&input, TX_4X8, DCT_DCT, 1, 1, 16);
    }

    #[test]
    fn rect_tx_8x4_dct_dct_roundtrip() {
        let input = lcg_residual(0x2468_ACE0_BDF1_1357, 32);
        check_roundtrip_frac(&input, TX_8X4, DCT_DCT, 1, 1, 16);
    }

    #[test]
    fn rect_tx_4x8_dc_input_dc_only_coefficient() {
        // A constant-DC input should produce a single dominant DC
        // coefficient (at index 0 of the row-major coefficient
        // buffer). All other coefficients should be at the LSB-floor
        // level from the kernel-chain rounding.
        let input = vec![64i64; 32];
        let coeffs = forward_transform_2d(&input, TX_4X8, DCT_DCT, false);
        assert_eq!(coeffs.len(), 32);
        let dc = coeffs[0];
        let max_off = coeffs[1..].iter().map(|c| c.abs()).max().unwrap_or(0);
        assert!(
            dc.abs() > 10 * max_off.max(1),
            "TX_4X8 DC = {dc}, max off-DC = {max_off} — DC should dominate"
        );
        // r460 — unit per-cell gain: each recovered cell ≈ the input.
        let recovered = inverse_transform_2d(&coeffs, TX_4X8, DCT_DCT, 8, false);
        for (i, &v) in recovered.iter().enumerate() {
            assert!(
                (v - 64).abs() <= 2,
                "TX_4X8 constant-DC round-trip cell {i}: got {v}, expected ≈ 64"
            );
        }
    }

    #[test]
    fn rect_tx_8x4_dc_input_dc_only_coefficient() {
        let input = vec![64i64; 32];
        let coeffs = forward_transform_2d(&input, TX_8X4, DCT_DCT, false);
        assert_eq!(coeffs.len(), 32);
        let dc = coeffs[0];
        let max_off = coeffs[1..].iter().map(|c| c.abs()).max().unwrap_or(0);
        assert!(
            dc.abs() > 10 * max_off.max(1),
            "TX_8X4 DC = {dc}, max off-DC = {max_off} — DC should dominate"
        );
        let recovered = inverse_transform_2d(&coeffs, TX_8X4, DCT_DCT, 8, false);
        for (i, &v) in recovered.iter().enumerate() {
            assert!(
                (v - 64).abs() <= 2,
                "TX_8X4 constant-DC round-trip cell {i}: got {v}, expected ≈ 64"
            );
        }
    }

    #[test]
    fn rect_tx_4x8_adst_dct_roundtrip() {
        // ADST is defined for n in 2..=4 (sizes 4 / 8 / 16); both
        // log2_w = 2 (n=2) and log2_h = 3 (n=3) are in range, so
        // ADST × DCT and DCT × ADST are reachable for TX_4X8.
        let input = lcg_residual(0xFACE_BABE_F00D_CAFE, 32);
        check_roundtrip_frac(&input, TX_4X8, ADST_DCT, 1, 1, 16);
    }

    #[test]
    fn rect_tx_8x4_dct_adst_roundtrip() {
        let input = lcg_residual(0xCAFE_F00D_BABE_FACE, 32);
        check_roundtrip_frac(&input, TX_8X4, DCT_ADST, 1, 1, 16);
    }

    #[test]
    fn rect_tx_4x8_idtx_roundtrip() {
        // IDTX is defined for n in 2..=5 (sizes 4 / 8 / 16 / 32);
        // both axes are in range for TX_4X8.
        let input = lcg_residual(0x1010_2020_3030_4040, 32);
        check_roundtrip_frac(&input, TX_4X8, IDTX, 1, 1, 16);
    }

    #[test]
    fn rect_tx_8x4_flipadst_flipadst_roundtrip() {
        // FLIPADST family — verify the §7.12.3 step-3 flip on both
        // axes is in effect (flipped input compared against the
        // recovered).
        let input = lcg_residual(0xDEAD_BEEF_CAFE_BABE, 32);
        check_roundtrip_flip(&input, TX_8X4, FLIPADST_FLIPADST, 1, 1, 16);
    }

    // -------------------------------------------------------------
    // Rectangular family — `|log2W - log2H| == 1` short-side-8 pair
    // (TX_8X16 / TX_16X8). Same encoder-mirror rectangular scale as
    // the short-side-4 pair (`× 2896` Round2 per row after the row
    // kernel). The empirical per-cell round-trip scale evaluates to
    // `1/2` (constant-DC probe: input = 32 ⇒ recovered = 16), one
    // step up from the short-side-4 pair's `1/4`. The larger
    // N_w × N_h kernel norm (128 vs 32) gains by `4×` over the
    // short-side-4 pair while the larger row_shift envelope
    // (`Transform_Row_Shift[TX_8X16] = 1` vs `0` for TX_4X8) loses
    // back `2×`, netting a `2×` larger round-trip per-cell scale.
    //
    // ADST is defined for `n in 2..=4` (sizes 4 / 8 / 16); both
    // log2_w = 3 and log2_h = 4 are in range for TX_8X16, and the
    // transpose (log2_w = 4, log2_h = 3) is in range for TX_16X8 —
    // so the full per-axis kernel matrix (DCT × DCT, ADST × DCT,
    // DCT × ADST, ADST × ADST, FLIPADST family, V_/H_ variants,
    // IDTX) is reachable on both shapes.
    //
    // For TX_16X8 the inverse-side between-stage `Clip3` at 16 bits
    // can saturate on full ±128 magnitudes; the roundtrip tests
    // scale down the input residual accordingly (`±32` keeps the
    // intermediate within the clamp).
    // -------------------------------------------------------------

    #[test]
    fn rect_tx_8x16_zero_input_yields_zero() {
        let input = vec![0i64; 8 * 16];
        let coeffs = forward_transform_2d(&input, TX_8X16, DCT_DCT, false);
        assert_eq!(coeffs.len(), 8 * 16);
        for (i, &v) in coeffs.iter().enumerate() {
            assert_eq!(v, 0, "TX_8X16 cell {i}: zero in ⇒ zero out, got {v}");
        }
    }

    #[test]
    fn rect_tx_16x8_zero_input_yields_zero() {
        let input = vec![0i64; 16 * 8];
        let coeffs = forward_transform_2d(&input, TX_16X8, DCT_DCT, false);
        assert_eq!(coeffs.len(), 16 * 8);
        for (i, &v) in coeffs.iter().enumerate() {
            assert_eq!(v, 0, "TX_16X8 cell {i}: zero in ⇒ zero out, got {v}");
        }
    }

    #[test]
    fn rect_tx_8x16_dct_dct_roundtrip() {
        // inverse pipeline's 16-bit between-stage clamp from
        // saturating on the 16-tall column kernel.
        let input = lcg_residual_bound(0xCAFE_8B16_DEAD_BEEF, 8 * 16, 32);
        check_roundtrip_frac(&input, TX_8X16, DCT_DCT, 1, 1, 16);
    }

    #[test]
    fn rect_tx_16x8_dct_dct_roundtrip() {
        let input = lcg_residual_bound(0xBEEF_16B8_F00D_CAFE, 16 * 8, 32);
        check_roundtrip_frac(&input, TX_16X8, DCT_DCT, 1, 1, 16);
    }

    #[test]
    fn rect_tx_8x16_dc_input_dc_only_coefficient() {
        // A constant-DC input should produce a single dominant DC
        // coefficient (at index 0 of the row-major coefficient
        // buffer); r460 unit gain: each recovered cell ≈ the input.
        let input = vec![32i64; 8 * 16];
        let coeffs = forward_transform_2d(&input, TX_8X16, DCT_DCT, false);
        assert_eq!(coeffs.len(), 8 * 16);
        let dc = coeffs[0];
        let max_off = coeffs[1..].iter().map(|c| c.abs()).max().unwrap_or(0);
        assert!(
            dc.abs() > 10 * max_off.max(1),
            "TX_8X16 DC = {dc}, max off-DC = {max_off} — DC should dominate"
        );
        let recovered = inverse_transform_2d(&coeffs, TX_8X16, DCT_DCT, 8, false);
        for (i, &v) in recovered.iter().enumerate() {
            assert!(
                (v - 32).abs() <= 2,
                "TX_8X16 constant-DC round-trip cell {i}: got {v}, expected ≈ 32"
            );
        }
    }

    #[test]
    fn rect_tx_16x8_dc_input_dc_only_coefficient() {
        let input = vec![32i64; 16 * 8];
        let coeffs = forward_transform_2d(&input, TX_16X8, DCT_DCT, false);
        assert_eq!(coeffs.len(), 16 * 8);
        let dc = coeffs[0];
        let max_off = coeffs[1..].iter().map(|c| c.abs()).max().unwrap_or(0);
        assert!(
            dc.abs() > 10 * max_off.max(1),
            "TX_16X8 DC = {dc}, max off-DC = {max_off} — DC should dominate"
        );
        let recovered = inverse_transform_2d(&coeffs, TX_16X8, DCT_DCT, 8, false);
        for (i, &v) in recovered.iter().enumerate() {
            assert!(
                (v - 32).abs() <= 2,
                "TX_16X8 constant-DC round-trip cell {i}: got {v}, expected ≈ 32"
            );
        }
    }

    #[test]
    fn rect_tx_8x16_adst_dct_roundtrip() {
        let input = lcg_residual_bound(0x8B16_AD0C_DEAD_BEEF, 8 * 16, 32);
        check_roundtrip_frac(&input, TX_8X16, ADST_DCT, 1, 1, 16);
    }

    #[test]
    fn rect_tx_16x8_dct_adst_roundtrip() {
        let input = lcg_residual_bound(0x16B8_DCAD_F00D_CAFE, 16 * 8, 32);
        check_roundtrip_frac(&input, TX_16X8, DCT_ADST, 1, 1, 16);
    }

    #[test]
    fn rect_tx_8x16_idtx_roundtrip() {
        // IDTX is defined for n in 2..=5 (sizes 4 / 8 / 16 / 32);
        // both log2_w = 3 (n=3) and log2_h = 4 (n=4) are in range.
        let input = lcg_residual_bound(0x8B16_1D70_BEEF_CAFE, 8 * 16, 32);
        check_roundtrip_frac(&input, TX_8X16, IDTX, 1, 1, 16);
    }

    #[test]
    fn rect_tx_16x8_flipadst_flipadst_roundtrip() {
        // FLIPADST family — verify the §7.12.3 step-3 flip on both
        // axes is in effect (flipped input compared against the
        // recovered).
        let input = lcg_residual_bound(0x16B8_F1F1_DEAD_BEEF, 16 * 8, 32);
        check_roundtrip_flip(&input, TX_16X8, FLIPADST_FLIPADST, 1, 1, 16);
    }

    #[test]
    fn rect_tx_8x16_adst_adst_roundtrip() {
        // Both axes ADST — log2_w = 3, log2_h = 4, both in
        // forward_adst_dispatch range (n in 2..=4).
        let input = lcg_residual_bound(0x8B16_AAAA_BABE_FACE, 8 * 16, 32);
        check_roundtrip_frac(&input, TX_8X16, ADST_ADST, 1, 1, 16);
    }

    // -------------------------------------------------------------
    // TX_16X32 / TX_32X16 — the |log2W - log2H| == 1 pair at
    // min(log2W, log2H) == 4 (short-side-16). Landed in round 241.
    //
    // ADST is defined for `n in 2..=4` (sizes 4 / 8 / 16); only the
    // log2 = 4 axis is in range, so axis combinations that demand
    // ADST on the length-32 side are forced to DCT by §6.10.19 and
    // are NOT exercised here. The reachable matrix is:
    //
    //   * DCT × DCT  (both axes DCT)            — DCT_DCT
    //   * DCT (row) × ADST (col, on the 16 side) — TX_32X16 with
    //                                              tx_type=ADST_DCT
    //                                              (col=ADST, length 16)
    //   * ADST (row, on the 16 side) × DCT (col) — TX_16X32 with
    //                                              tx_type=DCT_ADST
    //                                              (row=ADST, length 16)
    //   * IDTX × IDTX  (n in 2..=5 covers both lengths)
    //   * V_DCT / H_DCT  (one axis DCT, the other IDTX)
    //
    // The inverse-side between-stage `Clip3` at 16 bits saturates on
    // the length-32 column / row kernel for large input magnitudes;
    // the roundtrip tests use a reduced residual bound (`±8`) to
    // keep the intermediates within the clamp.
    // -------------------------------------------------------------

    #[test]
    fn rect_tx_16x32_zero_input_yields_zero() {
        let input = vec![0i64; 16 * 32];
        let coeffs = forward_transform_2d(&input, TX_16X32, DCT_DCT, false);
        assert_eq!(coeffs.len(), 16 * 32);
        for (i, &v) in coeffs.iter().enumerate() {
            assert_eq!(v, 0, "TX_16X32 cell {i}: zero in ⇒ zero out, got {v}");
        }
    }

    #[test]
    fn rect_tx_32x16_zero_input_yields_zero() {
        let input = vec![0i64; 32 * 16];
        let coeffs = forward_transform_2d(&input, TX_32X16, DCT_DCT, false);
        assert_eq!(coeffs.len(), 32 * 16);
        for (i, &v) in coeffs.iter().enumerate() {
            assert_eq!(v, 0, "TX_32X16 cell {i}: zero in ⇒ zero out, got {v}");
        }
    }

    #[test]
    fn rect_tx_16x32_dct_dct_roundtrip() {
        // inverse pipeline's 16-bit between-stage clamp from
        // saturating on the length-32 column kernel.
        let input = lcg_residual_bound(0xCAFE_1632_DEAD_BEEF, 16 * 32, 8);
        check_roundtrip_frac(&input, TX_16X32, DCT_DCT, 1, 1, 4);
    }

    #[test]
    fn rect_tx_32x16_dct_dct_roundtrip() {
        let input = lcg_residual_bound(0xBEEF_3216_F00D_CAFE, 32 * 16, 8);
        check_roundtrip_frac(&input, TX_32X16, DCT_DCT, 1, 1, 4);
    }

    #[test]
    fn rect_tx_16x32_dc_input_dc_only_coefficient() {
        // A constant-DC input should produce a single dominant DC
        // coefficient (at index 0 of the row-major coefficient
        // buffer); r460 unit gain: each recovered cell ≈ the input.
        let input = vec![8i64; 16 * 32];
        let coeffs = forward_transform_2d(&input, TX_16X32, DCT_DCT, false);
        assert_eq!(coeffs.len(), 16 * 32);
        let dc = coeffs[0];
        let max_off = coeffs[1..].iter().map(|c| c.abs()).max().unwrap_or(0);
        assert!(
            dc.abs() > 10 * max_off.max(1),
            "TX_16X32 DC = {dc}, max off-DC = {max_off} — DC should dominate"
        );
        let recovered = inverse_transform_2d(&coeffs, TX_16X32, DCT_DCT, 8, false);
        for (i, &v) in recovered.iter().enumerate() {
            assert!(
                (v - 8).abs() <= 2,
                "TX_16X32 constant-DC round-trip cell {i}: got {v}, expected ≈ 8"
            );
        }
    }

    #[test]
    fn rect_tx_32x16_dc_input_dc_only_coefficient() {
        let input = vec![8i64; 32 * 16];
        let coeffs = forward_transform_2d(&input, TX_32X16, DCT_DCT, false);
        assert_eq!(coeffs.len(), 32 * 16);
        let dc = coeffs[0];
        let max_off = coeffs[1..].iter().map(|c| c.abs()).max().unwrap_or(0);
        assert!(
            dc.abs() > 10 * max_off.max(1),
            "TX_32X16 DC = {dc}, max off-DC = {max_off} — DC should dominate"
        );
        let recovered = inverse_transform_2d(&coeffs, TX_32X16, DCT_DCT, 8, false);
        for (i, &v) in recovered.iter().enumerate() {
            assert!(
                (v - 8).abs() <= 2,
                "TX_32X16 constant-DC round-trip cell {i}: got {v}, expected ≈ 8"
            );
        }
    }

    #[test]
    fn rect_tx_16x32_dct_adst_roundtrip() {
        // TX_16X32: row pass length 16 (log2_w = 4), col pass length
        // 32 (log2_h = 5). To keep ADST on the length-16 axis only
        // (the §7.13.2.9 inverse-ADST dispatcher caps at n=4 ⇒
        // length 16), pick tx_type = DCT_ADST. Per
        // [`forward_row_kernel`] the row selector for DCT_ADST is
        // ADST (length 16, in range); per [`forward_col_kernel`] the
        // col selector for DCT_ADST is DCT (length 32, in DCT 2..=6
        // range). Both kernels reachable.
        let input = lcg_residual_bound(0x1632_AD0C_DEAD_BEEF, 16 * 32, 8);
        check_roundtrip_frac(&input, TX_16X32, DCT_ADST, 1, 1, 4);
    }

    #[test]
    fn rect_tx_32x16_adst_dct_roundtrip() {
        // TX_32X16: row pass length 32 (log2_w = 5), col pass length
        // 16 (log2_h = 4). To keep ADST on the length-16 axis only,
        // pick tx_type = ADST_DCT. Per [`forward_row_kernel`] the
        // row selector for ADST_DCT is DCT (length 32, OK); per
        // [`forward_col_kernel`] the col selector for ADST_DCT is
        // ADST (length 16, in §7.13.2.9 range).
        let input = lcg_residual_bound(0x3216_DCAD_F00D_CAFE, 32 * 16, 8);
        check_roundtrip_frac(&input, TX_32X16, ADST_DCT, 1, 1, 4);
    }

    #[test]
    fn rect_tx_16x32_idtx_roundtrip() {
        // IDTX is defined for n in 2..=5 (sizes 4 / 8 / 16 / 32);
        // both log2_w = 4 (n=4) and log2_h = 5 (n=5) are in range.
        // IDTX has its own per-axis scale envelope distinct from the
        // DCT / ADST kernels — empirical per-cell round-trip scale
        // is verified separately below.
        let input = lcg_residual_bound(0x1632_1D70_BEEF_CAFE, 16 * 32, 4);
        let coeffs = forward_transform_2d(&input, TX_16X32, IDTX, false);
        let recovered = inverse_transform_2d(&coeffs, TX_16X32, IDTX, 8, false);
        assert_eq!(recovered.len(), input.len());
        // Sanity probe: zero input ⇒ zero output (already in
        // zero_input_across_all_supported_square_combinations); non-
        // zero input round-trips to a finite per-cell scale.
        let nonzero = recovered.iter().any(|&v| v != 0);
        assert!(nonzero, "TX_16X32 IDTX round-trip produced all zeros");
    }

    #[test]
    fn rect_tx_32x16_idtx_roundtrip() {
        let input = lcg_residual_bound(0x3216_1D70_CAFE_BEEF, 32 * 16, 4);
        let coeffs = forward_transform_2d(&input, TX_32X16, IDTX, false);
        let recovered = inverse_transform_2d(&coeffs, TX_32X16, IDTX, 8, false);
        assert_eq!(recovered.len(), input.len());
        let nonzero = recovered.iter().any(|&v| v != 0);
        assert!(nonzero, "TX_32X16 IDTX round-trip produced all zeros");
    }

    #[test]
    fn rect_tx_16x32_v_dct_roundtrip() {
        // V_DCT: column kernel = DCT (length 32, in range), row
        // kernel = identity (length 16, in IDTX 2..=5 range).
        let input = lcg_residual_bound(0x1632_5DC7_BABE_F00D, 16 * 32, 8);
        check_roundtrip_frac(&input, TX_16X32, V_DCT, 1, 1, 4);
    }

    #[test]
    fn rect_tx_32x16_h_dct_roundtrip() {
        // H_DCT: row kernel = DCT (length 32, in range), column
        // kernel = identity (length 16, in IDTX 2..=5 range).
        let input = lcg_residual_bound(0x3216_4DC7_FACE_BABE, 32 * 16, 8);
        check_roundtrip_frac(&input, TX_32X16, H_DCT, 1, 1, 4);
    }

    // -------------------------------------------------------------
    // TX_32X64 / TX_64X32 — the |log2W - log2H| == 1 pair at
    // min(log2W, log2H) == 5 (short-side-32). Landed in round 244.
    //
    // Per-axis kernel-range coverage drives the reachable axis-pair
    // matrix:
    //
    //   * ADST: defined for `n in 2..=4` (sizes 4 / 8 / 16) — the
    //     §7.13.2.9 inverse-ADST dispatcher's range.
    //   * IDTX: defined for `n in 2..=5` (sizes 4 / 8 / 16 / 32) —
    //     the §7.13.2.15 inverse-identity dispatcher's range.
    //   * DCT: defined for `n in 2..=6` (sizes 4 / 8 / 16 / 32 /
    //     64) — covers both axes here.
    //
    // §6.10.19 maps each `tx_type` ordinal to the (row, col) kernel
    // pair; for length-32 / length-64 axis combinations, the spec
    // forces DCT on any axis whose kernel range excludes the axis
    // length. Reachable tx_types in this pair:
    //
    //   - TX_32X64: log2_w = 5 (length 32 — IDTX in range),
    //               log2_h = 6 (length 64 — DCT range).
    //     Reachable: DCT_DCT, V_DCT (col = DCT length 64, row =
    //     IDTX length 32).
    //   - TX_64X32: log2_w = 6 (length 64 — DCT range),
    //               log2_h = 5 (length 32 — IDTX in range).
    //     Reachable: DCT_DCT, H_DCT (row = DCT length 64, col =
    //     IDTX length 32). By transpose of the TX_32X64 case.
    //
    // The inverse-side between-stage `Clip3` at 16 bits saturates
    // even more aggressively on the length-64 column / row kernel;
    // the roundtrip tests use a tight residual bound (`±2`) to keep
    // the intermediates within the clamp.
    // -------------------------------------------------------------

    #[test]
    fn rect_tx_32x64_zero_input_yields_zero() {
        let input = vec![0i64; 32 * 64];
        let coeffs = forward_transform_2d(&input, TX_32X64, DCT_DCT, false);
        assert_eq!(coeffs.len(), 32 * 64);
        for (i, &v) in coeffs.iter().enumerate() {
            assert_eq!(v, 0, "TX_32X64 cell {i}: zero in ⇒ zero out, got {v}");
        }
    }

    #[test]
    fn rect_tx_64x32_zero_input_yields_zero() {
        let input = vec![0i64; 64 * 32];
        let coeffs = forward_transform_2d(&input, TX_64X32, DCT_DCT, false);
        assert_eq!(coeffs.len(), 64 * 32);
        for (i, &v) in coeffs.iter().enumerate() {
            assert_eq!(v, 0, "TX_64X32 cell {i}: zero in ⇒ zero out, got {v}");
        }
    }

    #[test]
    fn rect_tx_32x64_dct_dct_roundtrip() {
        // (`±2`) keeps the inverse pipeline's 16-bit between-stage
        // clamp from saturating on the length-64 column kernel; the
        // generous `max_err = 64` mirrors the TX_64X64 square test —
        // the deeper DCT-64 butterfly schedule accumulates more
        // `Round2(_, 12)` floor across the 31-step graph, and the
        // pseudo-random input drives the inverse's 16-bit between-
        // stage clamp closer to saturation than the constant-DC
        // probe does.
        let input = lcg_residual_bound(0xCAFE_3264_DEAD_BEEF, 32 * 64, 2);
        check_roundtrip_frac(&input, TX_32X64, DCT_DCT, 1, 1, 64);
    }

    #[test]
    fn rect_tx_64x32_dct_dct_roundtrip() {
        let input = lcg_residual_bound(0xBEEF_6432_F00D_CAFE, 64 * 32, 2);
        check_roundtrip_frac(&input, TX_64X32, DCT_DCT, 1, 1, 64);
    }

    #[test]
    fn rect_tx_32x64_dc_input_dc_only_coefficient() {
        // A constant-DC input should produce a single dominant DC
        // coefficient (at index 0 of the row-major coefficient
        // buffer); r460 unit gain: each recovered cell ≈ the input.
        let input = vec![2i64; 32 * 64];
        let coeffs = forward_transform_2d(&input, TX_32X64, DCT_DCT, false);
        assert_eq!(coeffs.len(), 32 * 64);
        let dc = coeffs[0];
        let max_off = coeffs[1..].iter().map(|c| c.abs()).max().unwrap_or(0);
        assert!(
            dc.abs() > 10 * max_off.max(1),
            "TX_32X64 DC = {dc}, max off-DC = {max_off} — DC should dominate"
        );
        let recovered = inverse_transform_2d(&coeffs, TX_32X64, DCT_DCT, 8, false);
        for (i, &v) in recovered.iter().enumerate() {
            assert!(
                (v - 2).abs() <= 2,
                "TX_32X64 constant-DC round-trip cell {i}: got {v}, expected ≈ 2"
            );
        }
    }

    #[test]
    fn rect_tx_64x32_dc_input_dc_only_coefficient() {
        let input = vec![2i64; 64 * 32];
        let coeffs = forward_transform_2d(&input, TX_64X32, DCT_DCT, false);
        assert_eq!(coeffs.len(), 64 * 32);
        let dc = coeffs[0];
        let max_off = coeffs[1..].iter().map(|c| c.abs()).max().unwrap_or(0);
        assert!(
            dc.abs() > 10 * max_off.max(1),
            "TX_64X32 DC = {dc}, max off-DC = {max_off} — DC should dominate"
        );
        let recovered = inverse_transform_2d(&coeffs, TX_64X32, DCT_DCT, 8, false);
        for (i, &v) in recovered.iter().enumerate() {
            assert!(
                (v - 2).abs() <= 2,
                "TX_64X32 constant-DC round-trip cell {i}: got {v}, expected ≈ 2"
            );
        }
    }

    #[test]
    fn rect_tx_32x64_v_dct_roundtrip() {
        // V_DCT on TX_32X64: column kernel = DCT (length 64, in
        // forward_dct_dispatch 2..=6 range), row kernel = identity
        // (length 32, in forward_idtx_dispatch 2..=5 range). The
        // IDTX kernel sits on the length-32 axis where it is
        // defined. Same loose `max_err` envelope as the DCT_DCT
        // case — the length-64 column kernel still dominates the
        // round-trip error envelope.
        let input = lcg_residual_bound(0x3264_5DC7_BABE_F00D, 32 * 64, 2);
        check_roundtrip_frac(&input, TX_32X64, V_DCT, 1, 1, 64);
    }

    #[test]
    fn rect_tx_64x32_h_dct_roundtrip() {
        // H_DCT on TX_64X32: row kernel = DCT (length 64, in
        // forward_dct_dispatch 2..=6 range), column kernel =
        // identity (length 32, in forward_idtx_dispatch 2..=5
        // range). The IDTX kernel sits on the length-32 axis.
        let input = lcg_residual_bound(0x6432_4DC7_FACE_BABE, 64 * 32, 2);
        check_roundtrip_frac(&input, TX_64X32, H_DCT, 1, 1, 64);
    }

    // -------------------------------------------------------------
    // TX_16X64 / TX_64X16 — the first `|log2W - log2H| == 2` pair,
    // at min(log2W, log2H) == 4 (short-side-16) and max(log2W,
    // log2H) == 6 (long-side-64). Landed in round 250.
    //
    // §5.11.40 `compute_tx_type()` short-circuits any tx_size with
    // `txSzSqrUp > TX_32X32` to DCT_DCT — i.e. these two shapes are
    // forced to DCT_DCT (the `TX_SET_DCTONLY` set, per §6.10.19
    // table). The reachable tx_type matrix is the single pair
    // {DCT_DCT}. The §7.13.3 `Round2(T[j] * 2896, 12)` row-pass
    // rectangular step does NOT fire on this shape (it gates on
    // |log2W - log2H| == 1; this family is at |log2W - log2H| == 2).
    //
    // Per-axis kernel reachability:
    //   * forward_dct_dispatch (n in 2..=6) covers both 16 (n=4)
    //     and 64 (n=6) — DCT × DCT is the only kernel pair the
    //     spec routes here.
    //   * forward_idtx_dispatch (n in 2..=5) covers 16 (n=4) only;
    //     length-64 axis is out of range — irrelevant because
    //     §6.10.19 has already forced DCT_DCT.
    //   * forward_adst_dispatch (n in 2..=4) covers 16 (n=4) only;
    //     length-64 axis is out of range — same point.
    //
    // §7.13.3 row/col shift envelope: Transform_Row_Shift[TX_16X64]
    // = Transform_Row_Shift[TX_64X16] = 2, col_shift = 4 — one step
    // deeper than the `|log2W - log2H| == 1` family's row_shift = 1
    // (matches the TX_32X32 / TX_64X64 square row_shift).
    //
    // Analytic per-cell round-trip scale (no rectangular factor):
    // `(N_w * N_h) / 2^(row_shift + col_shift)` = `1024 / 2^(2+4)` =
    // `16`. Empirical constant-DC scale is `4` — the documented `4×`
    // input-mass factor across the rectangular family.
    //
    // The inverse-side between-stage `Clip3` at 16 bits saturates
    // sharply on the length-64 kernel; the roundtrip tests use a
    // tight residual bound (`±2`) like the r244 TX_32X64 / TX_64X32
    // pair.
    // -------------------------------------------------------------

    #[test]
    fn rect_tx_16x64_zero_input_yields_zero() {
        let input = vec![0i64; 16 * 64];
        let coeffs = forward_transform_2d(&input, TX_16X64, DCT_DCT, false);
        assert_eq!(coeffs.len(), 16 * 64);
        for (i, &v) in coeffs.iter().enumerate() {
            assert_eq!(v, 0, "TX_16X64 cell {i}: zero in ⇒ zero out, got {v}");
        }
    }

    #[test]
    fn rect_tx_64x16_zero_input_yields_zero() {
        let input = vec![0i64; 64 * 16];
        let coeffs = forward_transform_2d(&input, TX_64X16, DCT_DCT, false);
        assert_eq!(coeffs.len(), 64 * 16);
        for (i, &v) in coeffs.iter().enumerate() {
            assert_eq!(v, 0, "TX_64X16 cell {i}: zero in ⇒ zero out, got {v}");
        }
    }

    #[test]
    fn rect_tx_16x64_dct_dct_roundtrip() {
        // (`±2`) keeps the inverse pipeline's 16-bit between-stage
        // clamp from saturating on the length-64 column kernel; the
        // generous `max_err = 64` mirrors the TX_64X64 / TX_32X64 /
        // TX_64X32 envelopes — the deeper DCT-64 butterfly schedule
        // accumulates more `Round2(_, 12)` floor across the 31-step
        // graph, and the pseudo-random input drives the inverse's
        // 16-bit between-stage clamp closer to saturation than the
        // constant-DC probe does.
        let input = lcg_residual_bound(0xCAFE_1664_DEAD_BEEF, 16 * 64, 2);
        check_roundtrip_frac(&input, TX_16X64, DCT_DCT, 1, 1, 64);
    }

    #[test]
    fn rect_tx_64x16_dct_dct_roundtrip() {
        let input = lcg_residual_bound(0xBEEF_6416_F00D_CAFE, 64 * 16, 2);
        check_roundtrip_frac(&input, TX_64X16, DCT_DCT, 1, 1, 64);
    }

    #[test]
    fn rect_tx_16x64_dc_input_dc_only_coefficient() {
        // A constant-DC input should produce a single dominant DC
        // coefficient (at index 0 of the row-major coefficient
        // buffer); r460 unit gain: each recovered cell ≈ the input.
        let input = vec![4i64; 16 * 64];
        let coeffs = forward_transform_2d(&input, TX_16X64, DCT_DCT, false);
        assert_eq!(coeffs.len(), 16 * 64);
        let dc = coeffs[0];
        let max_off = coeffs[1..].iter().map(|c| c.abs()).max().unwrap_or(0);
        assert!(
            dc.abs() > 10 * max_off.max(1),
            "TX_16X64 DC = {dc}, max off-DC = {max_off} — DC should dominate"
        );
        let recovered = inverse_transform_2d(&coeffs, TX_16X64, DCT_DCT, 8, false);
        for (i, &v) in recovered.iter().enumerate() {
            assert!(
                (v - 4).abs() <= 2,
                "TX_16X64 constant-DC round-trip cell {i}: got {v}, expected ≈ 4"
            );
        }
    }

    #[test]
    fn rect_tx_64x16_dc_input_dc_only_coefficient() {
        let input = vec![4i64; 64 * 16];
        let coeffs = forward_transform_2d(&input, TX_64X16, DCT_DCT, false);
        assert_eq!(coeffs.len(), 64 * 16);
        let dc = coeffs[0];
        let max_off = coeffs[1..].iter().map(|c| c.abs()).max().unwrap_or(0);
        assert!(
            dc.abs() > 10 * max_off.max(1),
            "TX_64X16 DC = {dc}, max off-DC = {max_off} — DC should dominate"
        );
        let recovered = inverse_transform_2d(&coeffs, TX_64X16, DCT_DCT, 8, false);
        for (i, &v) in recovered.iter().enumerate() {
            assert!(
                (v - 4).abs() <= 2,
                "TX_64X16 constant-DC round-trip cell {i}: got {v}, expected ≈ 4"
            );
        }
    }

    // -------------------------------------------------------------
    // TX_4X16 / TX_16X4 — the second `|log2W - log2H| == 2` pair, at
    // min(log2W, log2H) == 2 (short-side-4) and max(log2W, log2H) ==
    // 4 (long-side-16). Landed in round 251.
    //
    // §5.11.40 `compute_tx_type()` does NOT short-circuit on this
    // pair — `Tx_Size_Sqr_Up[TX_4X16] = Tx_Size_Sqr_Up[TX_16X4] =
    // TX_16X16` ≤ TX_32X32, so the full §5.11.48 / §6.10.19 tx_type
    // matrix is reachable: DCT_DCT, ADST_DCT, DCT_ADST, ADST_ADST,
    // FLIPADST_DCT, DCT_FLIPADST, FLIPADST_FLIPADST, ADST_FLIPADST,
    // FLIPADST_ADST, IDTX, V_DCT, H_DCT, V_ADST, H_ADST, V_FLIPADST,
    // H_FLIPADST — 16 ordinals all in play.
    //
    // §7.13.3 `Round2(T[j] * 2896, 12)` rectangular row-pass step
    // does NOT fire (the gate is exactly |log2W - log2H| == 1; this
    // pair is at |log2W - log2H| == 2).
    //
    // Per-axis kernel reachability:
    //   * forward_dct_dispatch (n in 2..=6) covers both 4 (n=2) and
    //     16 (n=4) ⇒ full DCT coverage on both axes.
    //   * forward_adst_dispatch (n in 2..=4) covers both 4 (n=2) and
    //     16 (n=4) ⇒ full ADST coverage on both axes.
    //   * forward_idtx_dispatch (n in 2..=5) covers both 4 (n=2) and
    //     16 (n=4) ⇒ full IDTX coverage on both axes.
    //
    // §7.13.3 row/col shift envelope: Transform_Row_Shift[TX_4X16] =
    // Transform_Row_Shift[TX_16X4] = 1, col_shift = 4.
    //
    // §7.12.3 `dqDenom` for this pair is 1 (neither TX_4X16 nor
    // TX_16X4 is listed in the dqDenom = 2 set or the dqDenom = 4
    // set per av1-spec p.294), so the existing dqDenom = 1 default
    // arm in `super::forward_quantize` covers it.
    //
    // Analytic per-cell round-trip scale (no rectangular factor):
    // `(N_w * N_h) / 2^(row_shift + col_shift)` = `64 / 2^(1+4)` =
    // `64 / 32` = `2`. Empirical per-cell round-trip scale is `1/2`
    // — the `4×` constant-DC input-mass factor documented across the
    // rectangular family (analytic `2` → empirical `1/2`).
    //
    // The inverse-side between-stage `Clip3` at 16 bits has plenty
    // of headroom on this small shape; tests use the full ±128
    // input range like the smaller square sizes.
    // -------------------------------------------------------------

    #[test]
    fn rect_tx_4x16_zero_input_yields_zero() {
        let input = vec![0i64; 4 * 16];
        let coeffs = forward_transform_2d(&input, TX_4X16, DCT_DCT, false);
        assert_eq!(coeffs.len(), 4 * 16);
        for (i, &v) in coeffs.iter().enumerate() {
            assert_eq!(v, 0, "TX_4X16 cell {i}: zero in ⇒ zero out, got {v}");
        }
    }

    #[test]
    fn rect_tx_16x4_zero_input_yields_zero() {
        let input = vec![0i64; 16 * 4];
        let coeffs = forward_transform_2d(&input, TX_16X4, DCT_DCT, false);
        assert_eq!(coeffs.len(), 16 * 4);
        for (i, &v) in coeffs.iter().enumerate() {
            assert_eq!(v, 0, "TX_16X4 cell {i}: zero in ⇒ zero out, got {v}");
        }
    }

    #[test]
    fn rect_tx_4x16_dct_dct_roundtrip() {
        // ±128 residual envelope — the length-16 axis is far below
        // the inverse-side 16-bit between-stage clamp's saturation
        // boundary, so the full per-cell residual range is safe.
        // Empirical per-cell round-trip scale = 1/2.
        let input = lcg_residual_bound(0x416B_EEF1_DEAD_F00D, 4 * 16, 128);
        check_roundtrip_frac(&input, TX_4X16, DCT_DCT, 1, 1, 6);
    }

    #[test]
    fn rect_tx_16x4_dct_dct_roundtrip() {
        let input = lcg_residual_bound(0x164B_EEF2_F00D_CAFE, 16 * 4, 128);
        check_roundtrip_frac(&input, TX_16X4, DCT_DCT, 1, 1, 6);
    }

    #[test]
    fn rect_tx_4x16_adst_adst_roundtrip() {
        // ADST × ADST is reachable on this pair — n=2 (length 4) and
        // n=4 (length 16) are both in forward_adst_dispatch range.
        let input = lcg_residual_bound(0x4A6B_EEF3_C0DE_BABE, 4 * 16, 128);
        check_roundtrip_frac(&input, TX_4X16, ADST_ADST, 1, 1, 6);
    }

    #[test]
    fn rect_tx_16x4_adst_adst_roundtrip() {
        let input = lcg_residual_bound(0x16AB_EEF4_DECA_F00D, 16 * 4, 128);
        check_roundtrip_frac(&input, TX_16X4, ADST_ADST, 1, 1, 6);
    }

    #[test]
    fn rect_tx_4x16_idtx_roundtrip() {
        // IDTX on both axes — n=2 and n=4 both in
        // forward_idtx_dispatch range. The per-axis identity scale
        // for n=2 is `× 5793` and for n=4 is `× 11586` (§7.13.2.11
        // / §7.13.2.13); composed with the row/col shift envelope
        // the empirical round-trip per-cell scale matches the DCT
        // arm at 1/2.
        let input = lcg_residual_bound(0x40C0_FFEE_BEEF_1664, 4 * 16, 128);
        check_roundtrip_frac(&input, TX_4X16, IDTX, 1, 1, 6);
    }

    #[test]
    fn rect_tx_16x4_idtx_roundtrip() {
        let input = lcg_residual_bound(0x16C0_FFEE_BEEF_BABE, 16 * 4, 128);
        check_roundtrip_frac(&input, TX_16X4, IDTX, 1, 1, 6);
    }

    #[test]
    fn rect_tx_4x16_flipadst_flipadst_roundtrip() {
        // FLIPADST × FLIPADST: same ADST butterfly, but the encoder
        // pre-flips the spatial residual on both axes and the
        // decoder counter-flips during the §7.12.3 frame-write — the
        // pre-flip + post-flip cancels round-trip when measured
        // against `flip(input)` (see `check_roundtrip_flip`). Here we
        // accept that the inverse pipeline emits the flipped signal
        // directly, so we compare against the doubly-flipped input
        // via `check_roundtrip_flip`.
        let input = lcg_residual_bound(0x4F1A_DEAD_C0DE_4F00, 4 * 16, 128);
        check_roundtrip_flip(&input, TX_4X16, FLIPADST_FLIPADST, 1, 1, 6);
    }

    #[test]
    fn rect_tx_16x4_v_dct_roundtrip() {
        // V_DCT on TX_16X4: column kernel = DCT (length 4, in DCT
        // dispatcher range 2..=6), row kernel = IDTX (length 16, in
        // IDTX dispatcher range 2..=5). Reachable because §5.11.40
        // does not short-circuit on this shape.
        let input = lcg_residual_bound(0x164D_C742_BABE_C0DE, 16 * 4, 128);
        check_roundtrip_frac(&input, TX_16X4, V_DCT, 1, 1, 6);
    }

    #[test]
    fn rect_tx_4x16_h_adst_roundtrip() {
        // H_ADST on TX_4X16: row kernel = ADST (length 4 — in ADST
        // dispatcher range 2..=4), column kernel = IDTX (length 16 —
        // in IDTX dispatcher range 2..=5). Reachable.
        let input = lcg_residual_bound(0x4A57_DEAD_BEEF_416B, 4 * 16, 128);
        check_roundtrip_frac(&input, TX_4X16, H_ADST, 1, 1, 6);
    }

    #[test]
    fn rect_tx_4x16_dc_input_dc_only_coefficient() {
        // Constant input of 4 ⇒ empirical per-cell round-trip scale
        // of 1/2 ⇒ each recovered cell ≈ 2. DC coefficient dominates
        // off-DC by at least 10×.
        let input = vec![4i64; 4 * 16];
        let coeffs = forward_transform_2d(&input, TX_4X16, DCT_DCT, false);
        assert_eq!(coeffs.len(), 4 * 16);
        let dc = coeffs[0];
        let max_off = coeffs[1..].iter().map(|c| c.abs()).max().unwrap_or(0);
        assert!(
            dc.abs() > 10 * max_off.max(1),
            "TX_4X16 DC = {dc}, max off-DC = {max_off} — DC should dominate"
        );
        let recovered = inverse_transform_2d(&coeffs, TX_4X16, DCT_DCT, 8, false);
        for (i, &v) in recovered.iter().enumerate() {
            assert!(
                (v - 4).abs() <= 1,
                "TX_4X16 constant-DC round-trip cell {i}: got {v}, expected ≈ 4"
            );
        }
    }

    #[test]
    fn rect_tx_16x4_dc_input_dc_only_coefficient() {
        let input = vec![4i64; 16 * 4];
        let coeffs = forward_transform_2d(&input, TX_16X4, DCT_DCT, false);
        assert_eq!(coeffs.len(), 16 * 4);
        let dc = coeffs[0];
        let max_off = coeffs[1..].iter().map(|c| c.abs()).max().unwrap_or(0);
        assert!(
            dc.abs() > 10 * max_off.max(1),
            "TX_16X4 DC = {dc}, max off-DC = {max_off} — DC should dominate"
        );
        let recovered = inverse_transform_2d(&coeffs, TX_16X4, DCT_DCT, 8, false);
        for (i, &v) in recovered.iter().enumerate() {
            assert!(
                (v - 4).abs() <= 1,
                "TX_16X4 constant-DC round-trip cell {i}: got {v}, expected ≈ 4"
            );
        }
    }

    // -------------------------------------------------------------
    // TX_8X32 / TX_32X8 — the third (mid) `|log2W - log2H| == 2`
    // pair, at min(log2W, log2H) == 3 (short-side-8) and
    // max(log2W, log2H) == 5 (long-side-32). Landed in round 252.
    //
    // §5.11.40 `compute_tx_type()` does NOT short-circuit on this
    // pair — `Tx_Size_Sqr_Up[TX_8X32] = Tx_Size_Sqr_Up[TX_32X8] =
    // TX_32X32` and the gate is the strict `txSzSqrUp > TX_32X32`
    // (av1-spec p.92), so DCT_DCT is NOT pre-forced. However the
    // §5.11.48 `get_tx_set()` routing IS gated by the equality
    // `txSzSqrUp == TX_32X32`:
    //   * intra arm: `txSzSqrUp == TX_32X32 ⇒ TX_SET_DCTONLY`,
    //     whose §6.10.19 inversion table (omitting all non-DCT
    //     ordinals) admits only DCT_DCT.
    //   * inter arm: `txSzSqrUp == TX_32X32 ⇒ TX_SET_INTER_3`
    //     (also taken when `reduced_tx_set` is set on any size),
    //     whose §6.10.19 inversion table is
    //     `Tx_Type_Inter_Inv_Set3 = { IDTX, DCT_DCT }`.
    //
    // So the reachable tx_type matrix on TX_8X32 / TX_32X8 is
    // exactly `{DCT_DCT, IDTX}` — narrower than the TX_4X16 /
    // TX_16X4 pair's full 16-ordinal matrix but wider than the
    // TX_16X64 / TX_64X16 pair's single-ordinal pin.
    //
    // §7.13.3 `Round2(T[j] * 2896, 12)` rectangular row-pass step
    // does NOT fire (the gate is exactly |log2W - log2H| == 1; this
    // pair is at |log2W - log2H| == 2).
    //
    // Per-axis kernel reachability:
    //   * forward_dct_dispatch (n in 2..=6) covers both 8 (n=3) and
    //     32 (n=5) ⇒ full DCT coverage on both axes (the kernel
    //     pair that runs for DCT_DCT).
    //   * forward_idtx_dispatch (n in 2..=5) covers both 8 (n=3) and
    //     32 (n=5) ⇒ full IDTX coverage on both axes (the kernel
    //     pair that runs for IDTX).
    //   * forward_adst_dispatch (n in 2..=4) covers 8 (n=3) only;
    //     length-32 axis is OUT of range — irrelevant because the
    //     §5.11.48 tx-type sets reachable on this pair do not
    //     contain any ADST ordinal.
    //
    // §7.13.3 row/col shift envelope: Transform_Row_Shift[TX_8X32] =
    // Transform_Row_Shift[TX_32X8] = 2, col_shift = 4.
    //
    // §7.12.3 `dqDenom` for this pair is 1 — neither TX_8X32 nor
    // TX_32X8 is in the `dqDenom = 2` set `{TX_32X32, TX_16X32,
    // TX_32X16, TX_16X64, TX_64X16}` nor in the `dqDenom = 4` set
    // `{TX_64X64, TX_32X64, TX_64X32}` per av1-spec p.294, so the
    // existing `dqDenom = 1` default arm in `super::forward_quantize`
    // covers it.
    //
    // Analytic per-cell round-trip scale (no rectangular factor):
    // `(N_w * N_h) / 2^(row_shift + col_shift)` = `256 / 2^(2+4)` =
    // `4`. Empirical per-cell round-trip scale is `1` — the `4×`
    // constant-DC input-mass factor documented across the
    // rectangular family (analytic `4` → empirical `1`).
    //
    // The inverse-side between-stage `Clip3` at 16 bits has plenty
    // of headroom on the length-32 axis (well below the saturation
    // boundary that constrained TX_16X64 / TX_64X16 / TX_32X64 /
    // TX_64X32 to a tight `±2` envelope); we use the same envelope
    // as the §7.13.3 length-32 kernels in the |log2W - log2H| == 1
    // family (TX_16X32 / TX_32X16 at `±32`).
    // -------------------------------------------------------------

    #[test]
    fn rect_tx_8x32_zero_input_yields_zero() {
        let input = vec![0i64; 8 * 32];
        let coeffs = forward_transform_2d(&input, TX_8X32, DCT_DCT, false);
        assert_eq!(coeffs.len(), 8 * 32);
        for (i, &v) in coeffs.iter().enumerate() {
            assert_eq!(v, 0, "TX_8X32 cell {i}: zero in ⇒ zero out, got {v}");
        }
    }

    #[test]
    fn rect_tx_32x8_zero_input_yields_zero() {
        let input = vec![0i64; 32 * 8];
        let coeffs = forward_transform_2d(&input, TX_32X8, DCT_DCT, false);
        assert_eq!(coeffs.len(), 32 * 8);
        for (i, &v) in coeffs.iter().enumerate() {
            assert_eq!(v, 0, "TX_32X8 cell {i}: zero in ⇒ zero out, got {v}");
        }
    }

    #[test]
    fn rect_tx_8x32_dct_dct_roundtrip() {
        // ±32 residual envelope — matches the length-32 envelope
        // used on TX_16X32 / TX_32X16 (round 241). The empirical
        // per-cell round-trip scale is 1 ⇒ no scale division below.
        let input = lcg_residual_bound(0x8B32_DEAD_C0DE_BEEF, 8 * 32, 32);
        check_roundtrip_frac(&input, TX_8X32, DCT_DCT, 1, 1, 64);
    }

    #[test]
    fn rect_tx_32x8_dct_dct_roundtrip() {
        let input = lcg_residual_bound(0x32B8_BEEF_F00D_CAFE, 32 * 8, 32);
        check_roundtrip_frac(&input, TX_32X8, DCT_DCT, 1, 1, 64);
    }

    #[test]
    fn rect_tx_8x32_idtx_roundtrip() {
        // IDTX is in `TX_SET_INTER_3` so it is reachable on
        // TX_8X32 (inter arm). Both per-axis lengths — 8 (n=3) and
        // 32 (n=5) — are in `forward_idtx_dispatch n in 2..=5`.
        let input = lcg_residual_bound(0x8C0F_FEEB_EEF8_3200, 8 * 32, 32);
        check_roundtrip_frac(&input, TX_8X32, IDTX, 1, 1, 64);
    }

    #[test]
    fn rect_tx_32x8_idtx_roundtrip() {
        let input = lcg_residual_bound(0x32C0_FFEE_BEEF_8800, 32 * 8, 32);
        check_roundtrip_frac(&input, TX_32X8, IDTX, 1, 1, 64);
    }

    #[test]
    fn rect_tx_8x32_dc_input_dc_only_coefficient() {
        // Constant input of 4 ⇒ empirical per-cell round-trip scale
        // of 1 ⇒ each recovered cell ≈ 4. DC coefficient dominates
        // off-DC by at least 10×.
        let input = vec![4i64; 8 * 32];
        let coeffs = forward_transform_2d(&input, TX_8X32, DCT_DCT, false);
        assert_eq!(coeffs.len(), 8 * 32);
        let dc = coeffs[0];
        let max_off = coeffs[1..].iter().map(|c| c.abs()).max().unwrap_or(0);
        assert!(
            dc.abs() > 10 * max_off.max(1),
            "TX_8X32 DC = {dc}, max off-DC = {max_off} — DC should dominate"
        );
        let recovered = inverse_transform_2d(&coeffs, TX_8X32, DCT_DCT, 8, false);
        for (i, &v) in recovered.iter().enumerate() {
            assert!(
                (v - 4).abs() <= 1,
                "TX_8X32 constant-DC round-trip cell {i}: got {v}, expected ≈ 4"
            );
        }
    }

    #[test]
    fn rect_tx_32x8_dc_input_dc_only_coefficient() {
        let input = vec![4i64; 32 * 8];
        let coeffs = forward_transform_2d(&input, TX_32X8, DCT_DCT, false);
        assert_eq!(coeffs.len(), 32 * 8);
        let dc = coeffs[0];
        let max_off = coeffs[1..].iter().map(|c| c.abs()).max().unwrap_or(0);
        assert!(
            dc.abs() > 10 * max_off.max(1),
            "TX_32X8 DC = {dc}, max off-DC = {max_off} — DC should dominate"
        );
        let recovered = inverse_transform_2d(&coeffs, TX_32X8, DCT_DCT, 8, false);
        for (i, &v) in recovered.iter().enumerate() {
            assert!(
                (v - 4).abs() <= 1,
                "TX_32X8 constant-DC round-trip cell {i}: got {v}, expected ≈ 4"
            );
        }
    }

    // Edge-case: zero input across the matrix.

    #[test]
    fn zero_input_across_all_supported_square_combinations() {
        let cases: &[(usize, usize, usize)] = &[
            (TX_4X4, 16, DCT_DCT),
            (TX_8X8, 64, DCT_DCT),
            (TX_16X16, 256, DCT_DCT),
            (TX_32X32, 1024, DCT_DCT),
            (TX_64X64, 4096, DCT_DCT),
            // Rectangular pair landed in round 235.
            (TX_4X8, 32, DCT_DCT),
            (TX_8X4, 32, DCT_DCT),
            (TX_4X8, 32, ADST_DCT),
            (TX_8X4, 32, DCT_ADST),
            (TX_4X8, 32, IDTX),
            (TX_8X4, 32, FLIPADST_FLIPADST),
            // Rectangular pair landed in round 238.
            (TX_8X16, 128, DCT_DCT),
            (TX_16X8, 128, DCT_DCT),
            (TX_8X16, 128, ADST_ADST),
            (TX_16X8, 128, DCT_ADST),
            (TX_8X16, 128, IDTX),
            (TX_16X8, 128, FLIPADST_FLIPADST),
            // Rectangular pair landed in round 241.
            (TX_16X32, 512, DCT_DCT),
            (TX_32X16, 512, DCT_DCT),
            // ADST is reachable only on the length-16 axis (the
            // §7.13.2.9 dispatcher caps at n=4). TX_16X32 needs the
            // row selector to pick ADST (row = length-16 axis) ⇒
            // tx_type = DCT_ADST; TX_32X16 needs the col selector to
            // pick ADST (col = length-16 axis) ⇒ tx_type = ADST_DCT.
            (TX_16X32, 512, DCT_ADST),
            (TX_32X16, 512, ADST_DCT),
            (TX_16X32, 512, IDTX),
            (TX_32X16, 512, IDTX),
            (TX_16X32, 512, V_DCT),
            (TX_32X16, 512, H_DCT),
            // Rectangular pair landed in round 244.
            (TX_32X64, 2048, DCT_DCT),
            (TX_64X32, 2048, DCT_DCT),
            // IDTX is reachable only on the length-32 axis (the
            // §7.13.2.15 dispatcher caps at n=5). TX_32X64 needs the
            // row selector to pick IDTX (row = length-32 axis) ⇒
            // tx_type = V_DCT (col = DCT length 64, row = IDTX
            // length 32); TX_64X32 needs the col selector to pick
            // IDTX (col = length-32 axis) ⇒ tx_type = H_DCT.
            (TX_32X64, 2048, V_DCT),
            (TX_64X32, 2048, H_DCT),
            // Rectangular pair landed in round 250 — the first
            // `|log2W - log2H| == 2` pair. Per §5.11.40
            // `compute_tx_type()` (`txSzSqrUp > TX_32X32 ⇒ return
            // DCT_DCT`), DCT_DCT is the only reachable tx_type.
            (TX_16X64, 1024, DCT_DCT),
            (TX_64X16, 1024, DCT_DCT),
            // Rectangular pair landed in round 251 — the second
            // `|log2W - log2H| == 2` pair at min_log2 == 2 /
            // max_log2 == 4. Per §5.11.40 these sit at
            // `txSzSqrUp = TX_16X16 ≤ TX_32X32` so the short-circuit
            // does NOT fire — the full §6.10.19 tx_type matrix is
            // reachable. Both per-axis lengths (4 and 16) are
            // covered by all three per-axis dispatchers (DCT n in
            // 2..=6, ADST n in 2..=4, IDTX n in 2..=5) so every
            // tx_type ordinal exercises a real kernel pair.
            (TX_4X16, 64, DCT_DCT),
            (TX_16X4, 64, DCT_DCT),
            (TX_4X16, 64, ADST_ADST),
            (TX_16X4, 64, ADST_ADST),
            (TX_4X16, 64, ADST_DCT),
            (TX_16X4, 64, DCT_ADST),
            (TX_4X16, 64, FLIPADST_FLIPADST),
            (TX_16X4, 64, FLIPADST_DCT),
            (TX_4X16, 64, DCT_FLIPADST),
            (TX_16X4, 64, ADST_FLIPADST),
            (TX_4X16, 64, FLIPADST_ADST),
            (TX_16X4, 64, IDTX),
            (TX_4X16, 64, IDTX),
            (TX_4X16, 64, V_DCT),
            (TX_16X4, 64, V_DCT),
            (TX_4X16, 64, H_DCT),
            (TX_16X4, 64, H_DCT),
            (TX_4X16, 64, V_ADST),
            (TX_16X4, 64, H_ADST),
            (TX_4X16, 64, V_FLIPADST),
            (TX_16X4, 64, H_FLIPADST),
            // Rectangular pair landed in round 252 — the third (mid)
            // `|log2W - log2H| == 2` pair at min_log2 == 3 /
            // max_log2 == 5. Per §5.11.40 these sit at
            // `txSzSqrUp = TX_32X32`, so the strict `> TX_32X32`
            // short-circuit does NOT fire — but §5.11.48 routes
            // `txSzSqrUp == TX_32X32` to either `TX_SET_DCTONLY`
            // (intra arm: only DCT_DCT) or `TX_SET_INTER_3` (inter
            // arm: `{IDTX, DCT_DCT}`), giving a 2-ordinal reachable
            // tx_type set per shape. Both per-axis lengths (8 and
            // 32) are covered by DCT (n in 2..=6) and IDTX (n in
            // 2..=5).
            (TX_8X32, 256, DCT_DCT),
            (TX_32X8, 256, DCT_DCT),
            (TX_8X32, 256, IDTX),
            (TX_32X8, 256, IDTX),
            (TX_4X4, 16, ADST_ADST),
            (TX_8X8, 64, ADST_DCT),
            (TX_8X8, 64, DCT_ADST),
            (TX_16X16, 256, ADST_ADST),
            (TX_4X4, 16, FLIPADST_FLIPADST),
            (TX_8X8, 64, FLIPADST_DCT),
            (TX_16X16, 256, ADST_FLIPADST),
            (TX_4X4, 16, IDTX),
            (TX_8X8, 64, IDTX),
            (TX_16X16, 256, IDTX),
            (TX_32X32, 1024, IDTX),
            (TX_4X4, 16, V_DCT),
            (TX_4X4, 16, H_DCT),
            (TX_8X8, 64, V_ADST),
            (TX_16X16, 256, H_FLIPADST),
        ];
        for &(tx_size, n, tx_type) in cases {
            let input = vec![0i64; n];
            let coeffs = forward_transform_2d(&input, tx_size, tx_type, false);
            assert_eq!(coeffs.len(), n);
            for (i, &v) in coeffs.iter().enumerate() {
                assert_eq!(
                    v, 0,
                    "tx_size={tx_size} tx_type={tx_type} cell {i}: zero in ⇒ zero out, \
                     got {v}"
                );
            }
        }
    }

    // Panic guard tests.

    #[test]
    #[should_panic(expected = "out of range")]
    fn tx_size_above_tx_sizes_all_panics() {
        // Round 252 closes the rectangular arc — all twelve §5.11
        // rectangular `TX_SIZES_ALL` shapes are now supported, so
        // there is no longer a rectangular size to panic on. The
        // first dispatcher assert remains: a `tx_size` ordinal at
        // or beyond `TX_SIZES_ALL` is out of range and must still
        // panic. The smallest such value is `TX_SIZES_ALL` itself
        // (= 19 per av1-spec p.408 / §3); the input length is
        // arbitrary because the assert fires before the length
        // check.
        let input = vec![0i64; 16];
        let _ = forward_transform_2d(&input, crate::cdf::TX_SIZES_ALL, DCT_DCT, false);
    }

    #[test]
    #[should_panic(expected = "lossless arm requires tx_size = TX_4X4")]
    fn lossless_non_4x4_panics() {
        let input = vec![0i64; 64];
        let _ = forward_transform_2d(&input, TX_8X8, DCT_DCT, true);
    }
}