diff --git a/README.md b/README.md index 75f0b7e..358e66e 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.228 s - 2609x real time - 39.8 MB/s of bitstream, 230.1 M samples/s - 2503 ns per granule and channel + 595.59 s of audio in 0.222 s + 2688x real time + 41.0 MB/s of bitstream, 237.1 M samples/s + 2430 ns per granule and channel per granule and channel: - polyphase synthesis 756 ns 30.2% - IMDCT and windowing 744 ns 29.7% - alias reduction 51 ns 2.1% + polyphase synthesis 763 ns 31.4% + IMDCT and windowing 762 ns 31.4% + alias reduction 44 ns 1.8% frequency inversion 22 ns 0.9% - everything else 930 ns 37.1% + everything else 838 ns 34.5% ``` "Everything else" is Huffman decoding and requantization. It used to be two @@ -332,66 +332,105 @@ processor time. ### How it got there -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. +It started out ten times slower, and measurement — not intuition — found every +change 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. -The two filterbanks were the first target, and both had the same two problems. -Each was computing outputs the standard's own symmetries give away — writing -the synthesis matrixing as D[u] = Σ s[j]·cos(u(2j+1)π/64), the row the standard +**The filterbanks** were the first target, and both had the same two problems. +Each computed outputs that the standard's own symmetries give away: writing the +synthesis matrixing as D[u] = Σ s[j]·cos(u(2j+1)π/64), the row the standard calls N[i] is D[16 + i] and D[64 − u] = −D[u], so half the rows are free and D[32] is zero; the inverse MDCT's 36 outputs likewise obey x[17 − i] = −x[i] over the first half and x[53 − i] = x[i] over the second. And both were written as dot products, which end in a horizontal sum across a vector that costs more -than the multiplying did. Transposing the coefficient tables so the results -accumulate in separate lanes and are never summed across removed that: 16.8 µs -to 0.87 for the filterbank, 4.8 to 0.68 for the transform. - -Then the bit-level half. Requantization was branching on the sign of every -spectral value, which is as good as random and mispredicted half the time; -moving the sign across as a bit instead took it from 19% of the profile to 8%. -The gain term is a power of two in quarter steps, so a table replaced a libm -call. And the spectrum is no longer cleared before it is filled, nor -requantized above the last line the Huffman data actually coded — between them -`memset` fell from 5.5% of the profile to 0.6%. - -Four things that seemed obvious and were not, all of them measured and backed -out: - -- **Caching the bit reader's window** so that consecutive reads share one load. - Slightly slower: the branch to decide whether the window still covers the - read costs more than an L1 load. -- **The same branchless sign trick in `spectralValue`.** It helped in - requantization because it unblocked vectorising a 576-sample loop; here there - is no loop to vectorise and it only lengthened a serial dependency. 5% slower. -- **Hoisting the Huffman region and table lookup** out of the per-pair loop. - Exactly neutral. -- **Inlining the Huffman code word lookup.** 5% slower, from code bloat. - -A later pass took another tenth off, and mostly by removing work rather than -doing it faster. The synthesis history had been kept in a buffer of twice the -length it needed, every value written to both halves so that the window's -reads never straddled the wrap; keeping it as a ring of sixteen 64-value -blocks instead makes that impossible by construction, because every read the -standard asks for falls inside one block. Half the memory, and no mirrored -write. And the inverse MDCT now looks at a subband before transforming it: an -encoder writes nothing above its lowpass, so an eighth of them are silent on -broadband material and a fifth on anything resembling music, and the transform -of nothing is nothing. - -Two more things were tried and backed out. Taking the code word and both sign -bits from a single 32-bit peek, rather than three separate loads, was 8% -slower — the accumulated bit count lengthens the dependency chain before the -next iteration can start, and on this path latency costs more than loads do. -Widening the Huffman lookup table was not worth measuring once the tree walk -behind it was counted: it is taken for 2.7% of code words. +than the multiplying did — transposing the coefficient tables so that results +accumulate in separate lanes removed that entirely. + +Two more followed. The synthesis history had been kept in a buffer of twice the +length it needed, every value written to both halves so the window's reads +never straddled the wrap; a ring of sixteen 64-value blocks makes straddling +impossible by construction, because every read the standard asks for falls +inside one block. And building the 64 matrix rows out of the 32 computed +values was sixty-odd scalar moves per slot, where two shuffles and two sign +patterns do it in six vector operations. + +**Requantization** was branching on the sign of every spectral value, which is +as good as random and mispredicted half the time; moving the sign across as a +bit took it from 19% of the profile to 8%. Its gain term is a power of two in +quarter steps, so a table replaced a libm call. Almost every spectral line is +small enough to have been coded without an escape, so a second small table +carries those values with their sign already in them. And the spectrum is no +longer cleared before being filled, nor requantized above the last line the +Huffman data actually coded — `memset` fell from 5.5% of the profile to 0.6%. + +**The inverse MDCT** now looks at a subband before transforming it. An encoder +writes nothing above its lowpass, so an eighth of them are silent on the +broadband noise this is benchmarked with and a fifth on anything resembling +music, and the transform of nothing is nothing. + +**Huffman decoding** is a chain of dependent reads — each one's position +depends on the width of the last — so an eight-byte load sat on the critical +path of every one. It now runs against a cache of the next 57 bits held in a +register, refilled once per code word rather than touched three times. +Counting rather than guessing found two more: 39.5% of pairs come from a table +with an escape, so the test for it is a coin flip taken sixteen million times a +minute — and it is fixed for a whole region, so it now sits outside the loop. +And the count1 partition was working out four values from four magnitude bits +and a run of sign bits with four dependent tests, three million times a minute; +there are only 256 combinations, so they are all worked out at compile time. + +### What did not work + +Rather more than what did, and each was quicker to measure than to argue about. + +- **Caching the bit reader's window** inside `Reader`, so consecutive reads + share a load. Slower: the branch deciding whether the window still covers the + read costs more than the L1 load it saves. The same idea as a local buffer + over the hot loop, where the compiler can keep it in registers, is the win + described above. +- **The branchless sign trick in `spectralValue`.** It helped in requantization + because it unblocked vectorising a 576-sample loop; here there is no loop to + vectorise and it only lengthened a serial dependency. 5% slower. +- **Taking the code word and both sign bits from one 32-bit peek** rather than + three loads. 8% slower against the reader, because the peek was a load and + accumulating the bit count lengthened the chain. Against the register cache + the same shape is free, and is how the loop is written now. +- **Splitting the synthesis accumulator four ways** to break the multiply-add + dependency chain. 6% slower: a 32-lane multiply-add is already four + independent 8-lane ones, so the parallelism was there and four accumulators + only spilled registers. Two accumulators measured identically to one. +- **Matrixing two slots per pass** over the four-kilobyte coefficient table, on + the theory that the stage was bound by loads. No change at all. +- **Padding the inverse MDCT's vectors to 24 lanes** rather than 32, to waste + fewer on the 18 rows it needs. Five times slower — LLVM legalises a vector + that is not a power of two into pieces costing far more than the lanes they + save, the same trap an 18-lane vector fell into earlier. +- **Refilling the bit cache once per two code words** instead of once per word. + Slightly slower. +- **Widening the Huffman lookup table**, which was not worth measuring once the + tree walk behind it was counted: it is taken for 2.7% of code words. Both filterbanks are also written out longhand in `reference.zig`, the way the -standard's flowcharts read, and tests hold the fast versions against them. -That is what made all of this safe: every change was checked against a second -implementation, against ffmpeg's decoder, and against the fuzz suite, and the -output never moved by a bit. +standard's flowcharts read, with tests holding the fast versions against them. +That is what made all of this safe to attempt: every change was checked against +a second implementation, against ffmpeg's decoder, and against the fuzz suite. + +### On "the same output" + +The standard does not define one correct set of samples. ISO/IEC 11172-4 +defines compliance as staying within a tolerance of a reference decoder, and +every implementation in floating point lands somewhere inside it. What is +checked here — and what "agrees to within one least significant bit" means +throughout — is agreement with ffmpeg's decoder, not that this decoder's own +output never moves. + +It does move, occasionally and by one. The order a sum is accumulated in +changes its rounding, so splitting the inverse MDCT's accumulator shifted 386 +samples in 441,000 by a single count. That is far inside the tolerance and is +not a defect, but it is worth being exact about: these changes preserve +agreement with a reference, not byte-for-byte reproducibility against an +earlier build. ## Limitations diff --git a/src/layer3.zig b/src/layer3.zig index 1dfdb43..385360c 100644 --- a/src/layer3.zig +++ b/src/layer3.zig @@ -389,6 +389,23 @@ fn gainFor(exponent: i32) f32 { return gains[@intCast(at)]; } +/// v raised to 4/3, keeping v's sign, for the values that fit in a Huffman +/// code word without an escape. +/// +/// Almost every spectral line is one of these, and having the sign already +/// in the table saves taking the absolute value and putting the sign back +/// afterwards on each of the 576 lines of every granule. +const small_pow43 = blk: { + @setEvalBranchQuota(20_000); + var t: [31]f32 = undefined; + for (0..31) |n| { + const v: f64 = @floatFromInt(@as(i32, @intCast(n)) - 15); + const size = std.math.pow(f64, @abs(v), 4.0 / 3.0); + t[n] = @floatCast(if (v < 0) -size else size); + } + break :blk t; +}; + const pow43 = blk: { @setEvalBranchQuota(200_000); var t: [1024]f32 = undefined; @@ -439,13 +456,20 @@ fn requantize( for (line..line + count) |k| { const value = quantized[k]; - const size = magnitude(@abs(value)); - // The sign moves across as a bit rather than a branch. Spectral - // values change sign about as often as not and in no pattern, so - // branching on it is a coin flip the processor has to pay for - // every line. - const sign = @as(u32, @bitCast(value)) & 0x8000_0000; - out[k] = @as(f32, @bitCast(@as(u32, @bitCast(size)) | sign)) * gain; + // One unsigned comparison covers both ends of the small table: + // anything outside -15...15 wraps to something far above 30. + const biased: u32 = @bitCast(value +% 15); + if (biased <= 30) { + out[k] = small_pow43[biased] * gain; + } else { + // Only an escape coded value gets here. The sign moves + // across as a bit rather than a branch, since spectral values + // change sign about as often as not and in no pattern. + @branchHint(.unlikely); + const size = magnitude(@abs(value)); + const sign = @as(u32, @bitCast(value)) & 0x8000_0000; + out[k] = @as(f32, @bitCast(@as(u32, @bitCast(size)) | sign)) * gain; + } } line += count; } diff --git a/src/main.zig b/src/main.zig index 5686bfd..ff5e3ff 100644 --- a/src/main.zig +++ b/src/main.zig @@ -6,6 +6,8 @@ const std = @import("std"); const mp3 = @import("mp3"); +const native_endian = @import("builtin").cpu.arch.endian(); + const usage = \\usage: zig-mp3 [--raw] [output.wav] \\ @@ -245,7 +247,7 @@ fn convert( var decoder = mp3.Decoder.init(); var buffer: mp3.FrameBuffer = undefined; - var encoded: [mp3.max_frame_samples * 2]u8 = undefined; + var encoded: [mp3.max_frame_samples]i16 = undefined; var channels: u8 = 0; var sample_rate: u32 = 0; @@ -292,11 +294,15 @@ fn convert( } if (samples.len == 0) continue; - for (samples, 0..) |sample, n| { - const scaled = std.math.clamp(sample * 32768.0, -32768.0, 32767.0); - std.mem.writeInt(i16, encoded[n * 2 ..][0..2], @intFromFloat(scaled), .little); + // The library's conversion, which the compiler vectorises; + // writing the samples one at a time through `writeInt` did + // not, and cost a twentieth of the whole program. + mp3.toInt16(samples, &encoded); + const pcm = encoded[0..samples.len]; + if (native_endian != .little) { + for (pcm) |*sample| sample.* = @byteSwap(sample.*); } - const bytes = encoded[0 .. samples.len * 2]; + const bytes = std.mem.sliceAsBytes(pcm); if (to_stdout) { try held.appendSlice(gpa, bytes); } else {