diff --git a/README.md b/README.md index f444a50..75f0b7e 100644 --- a/README.md +++ b/README.md @@ -298,16 +298,16 @@ stages on their own: ``` $ zig build bench decoding 44100 Hz, 2 channels - 595.59 s of audio in 0.232 s - 2563x real time - 39.1 MB/s of bitstream, 226.1 M samples/s - 2548 ns per granule and channel + 595.59 s of audio in 0.228 s + 2609x real time + 39.8 MB/s of bitstream, 230.1 M samples/s + 2503 ns per granule and channel per granule and channel: - polyphase synthesis 858 ns 33.7% - IMDCT and windowing 772 ns 30.3% - alias reduction 52 ns 2.1% + polyphase synthesis 756 ns 30.2% + IMDCT and windowing 744 ns 29.7% + alias reduction 51 ns 2.1% frequency inversion 22 ns 0.9% - everything else 844 ns 33.1% + everything else 930 ns 37.1% ``` "Everything else" is Huffman decoding and requantization. It used to be two @@ -332,7 +332,7 @@ processor time. ### How it got there -It started out nine times slower than this, and measurement — not intuition — +It started out ten times slower than this, and measurement — not intuition — found every one of the changes that mattered. `perf` is in the devshell; because the whole granule pipeline inlines into one symbol, marking the stages `noinline` for a profiling run is what makes it legible. diff --git a/src/imdct.zig b/src/imdct.zig index 6fa1337..16bacb5 100644 --- a/src/imdct.zig +++ b/src/imdct.zig @@ -149,11 +149,14 @@ fn transform(in: *const [subband_len]f32, out: *[36]f32, block_type: BlockType) // The three windows' inputs are interleaved: window `w` holds the // values at 3k + w, which is how reordering left them. for (0..3) |w| { - var accumulator: ShortVector = @splat(0); - for (0..6) |k| { - accumulator += short_cos[k] * @as(ShortVector, @splat(in[3 * k + w])); + var sums: [2]ShortVector = @splat(@splat(0)); + for (0..3) |k| { + inline for (0..2) |lane| { + sums[lane] += short_cos[2 * k + lane] * + @as(ShortVector, @splat(in[3 * (2 * k + lane) + w])); + } } - const x: [short_lanes]f32 = accumulator; + const x: [short_lanes]f32 = sums[0] + sums[1]; const at = 6 * w + 6; for (0..3) |i| { out[at + i] += x[i] * window[2][i]; @@ -168,11 +171,18 @@ fn transform(in: *const [subband_len]f32, out: *[36]f32, block_type: BlockType) return; } - var accumulator: LongVector = @splat(0); - for (0..subband_len) |k| { - accumulator += long_cos[k] * @as(LongVector, @splat(in[k])); + // Three running totals rather than one: a multiply-add cannot start + // until the one before it has finished, so a single accumulator makes + // the chain decide how long this takes rather than the arithmetic. + // Eighteen inputs divide evenly by three. + var sums: [3]LongVector = @splat(@splat(0)); + var k: usize = 0; + while (k < subband_len) : (k += 3) { + inline for (0..3) |lane| { + sums[lane] += long_cos[k + lane] * @as(LongVector, @splat(in[k + lane])); + } } - const x: [long_lanes]f32 = accumulator; + const x: [long_lanes]f32 = (sums[0] + sums[1]) + sums[2]; const w = window[@intFromEnum(block_type)]; for (0..9) |i| { diff --git a/src/synthesis.zig b/src/synthesis.zig index 31af6a4..a2f40fc 100644 --- a/src/synthesis.zig +++ b/src/synthesis.zig @@ -37,10 +37,31 @@ const block_len = 64; const Block = @Vector(subbands, f32); -/// Lane order that reverses a `Block`. -const reverse = blk: { +/// The 64 rows of the standard's matrix are a rearrangement of the 32 values +/// actually computed: row i is D[16 + i], D[32] is zero, and above that D[u] +/// is -D[64 - u]. Written as two shuffles and two sign patterns, the whole +/// rearrangement is six vector operations instead of sixty-four moves. +const rows_low = blk: { var m: [subbands]i32 = undefined; - for (0..subbands) |i| m[i] = subbands - 1 - i; + var signs: [subbands]f32 = undefined; + for (0..subbands) |i| { + if (i < 16) { + m[i] = @intCast(16 + i); // D[16 + i] + signs[i] = 1; + } else if (i == 16) { + m[i] = 0; // D[32], which is zero; the sign pattern makes it so + signs[i] = 0; + } else { + m[i] = @intCast(48 - i); // -D[64 - (16 + i)] + signs[i] = -1; + } + } + break :blk .{ .order = m, .signs = @as(Block, signs) }; +}; + +const rows_high = blk: { + var m: [subbands]i32 = undefined; + for (0..subbands) |i| m[i] = @intCast(if (i < 16) 16 - i else i - 16); break :blk m; }; @@ -105,27 +126,15 @@ pub fn synthesize( const sample: Block = @splat(spectrum[j * slots + slot]); accumulator += cosines[j] * sample; } - const d: [32]f32 = accumulator; - state.head = (state.head + blocks - 1) & (blocks - 1); - // Row i of the standard's matrix is D[16 + i]. D[32] is zero, and - // above it D[u] is -D[64 - u], so the 64 rows come out of the 32 - // values just computed with nothing but a reflection and a sign. - // Rows 17 to 47 are the first 31 of D reversed and negated, which is - // a shuffle rather than 31 separate loads. - const flipped: [32]f32 = -@shuffle(f32, accumulator, undefined, reverse); - + // The 64 rows, as two shuffles of the 32 values just computed. const rows = state.v[state.head * block_len ..][0..block_len]; - @memcpy(rows[0..16], d[16..32]); - rows[16] = 0; // cos of an odd multiple of pi/2 - @memcpy(rows[17..48], flipped[0..31]); - for (0..16) |k| rows[48 + k] = -d[k]; + rows[0..subbands].* = @shuffle(f32, accumulator, undefined, rows_low.order) * + rows_low.signs; + rows[subbands..][0..subbands].* = + -@shuffle(f32, accumulator, undefined, rows_high); - // The standard builds a 512 element vector out of alternating half - // blocks of the history, multiplies it by the window and sums every - // 32nd element. Folding those together leaves sixteen multiply-adds - // of 32 lanes each. // The standard builds a 512 element vector out of alternating half // blocks of the history, multiplies it by the window and sums every // 32nd element. Folding those together leaves sixteen multiply-adds @@ -143,7 +152,7 @@ pub fn synthesize( out[at..][0..subbands].* = total; } else { const lanes: [subbands]f32 = total; - for (lanes, 0..) |sample, j| out[(at + j) * stride] = sample; + for (lanes, 0..) |sample, lane| out[(at + lane) * stride] = sample; } } }