The PSP has a vector unit, the VFPU, which is set up as a traditional MIPS co-processor that shares the instruction stream with the main CPU, but has its own pipeline. Much like modern GPUs it has a special function unit, for functions like sin/cos, exp2/log2, square root, reciprocals, etc. It has long been know exactly what these instructions do, but not how they do it - they all return approximations made in an unknown way, none of them return the correctly rounded IEEE result.
vrcp(reciprocal)vrsq(reciprocal of square root)vsqrt(square root)vexp2(pow(2.0f, input))vsin(sin, with 4.0 = one lap instead of 2*PI)vcos(cos, with 4.0 = one lap instead of 2*PI)vasin(arc-sine, same 4.0 circle system)vlog2
There's also vrot which is a convenience function for sin+cos for building rotation matrices.
fp64, a long-term contributor to the project, has already implemented bit-exact versions of these (#16984, etc), by dumping the full result tables from a real PSP and adding correction tables on top of close approximations. This costs multiple megabytes of memory though, and some performance - it would be much neater if we knew the underlying functions. fp64 has also implemented a bit-exact version of the VFPU's dot product instruction, which really deserves its own article.
Encouraged by skmp's work reversing similar instructions on the Sega Dreamcast, Reversing the SH-4 fpu with the help of three AI models, I thought I'd simply set Claude on it.
And lo and behold, it worked. It took Claude Opus 5.5 about an hour.
This will save PPSSPP a few megabytes of shipped tables in future versions, and we will enable accurate emulation of these functions by default due to the somewhat faster performance than the previous method.
The pull request: VFPU: Replace fp64's correction tables with accurate implementation of special functions.
Below is Claude's own writeup (very lightly edited) about how it did it.
Claude says #
PPSSPP has had the functions bit-exact for a while, thanks to fp64's work in issue #16946. The method was careful measurement: an easy-to-compute base approximation, plus per-64-input correction deltas, plus lists of exceptions. It shipped as about 4.9 MB of .dat binary files. That's the right way to get correct answers when you don't know the algorithm. It also leaves an itch: the hardware certainly doesn't have 5 MB of ROM for this. What does it actually do?
This post is about finding that out, one function at a time, without ever looking inside the chip. The answer turned out to be a single, rather textbook circuit. All eight functions now come from about 10 KB of coefficients, bit-exact over every one of the 2^32 float inputs.
The setup
fp64's code is an exact oracle. It agrees with the hardware on every input, and it runs in a second on a laptop. So I could dump every output of every function, for example all 2^23 mantissas of 1/x over [1, 2), and interrogate the data as much as I liked. Before relying on that, a new hardware test (pspautotests cpu/vfpu/exact) ran special values and exponent sweeps on a real PSP. It agreed with the oracle, and showed a few mismatches in the other CPU backends along the way.
The other tool that mattered was exhaustive checking. A hypothesis about a 23-bit function isn't right until it reproduces all 2^32 inputs. That takes 30 seconds, so there's no excuse.
First attempts: good enough isn't the answer
The obvious guess is a piecewise polynomial. Fit a quadratic to each chunk of 1/x and see how many outputs come out exactly right. The answer was about 93%, whatever chunk size or polynomial degree I tried. Fitting coefficients that are floats and evaluating them with the VFPU's own dot product, which has its own exact rounding rules, got about 92%. One suggestion from Henrik was that the functions are built out of vdot operations, since on the Dreamcast the FIPR inner-product unit turned out to be reused by the function approximator. That was a good lead, but a single vdot of float coefficients also stalled at 92%.
The failures had a consistent shape. The best possible real-valued quadratic missed by about one 24-bit ulp, in a band exactly one ulp wider than the output's truncation window. Something discrete was going on underneath, and least-squares fitting would never find it. I switched to feasibility instead of fitting: for a candidate formula, can any coefficients reproduce every output exactly? Each output then turns into an interval constraint on the unknowns, and intervals intersect cheaply. This question, rather than "how close can I get", drove everything after it.
Finding the segments
The first structural question was how many pieces the function is made of. Within every run of 64 consecutive inputs, the output is exactly a straight line. Past that, I checked how many of those 64-input intervals can share one slope. The answer was 1,024 intervals, which is 2^16 inputs, and never 2,048. So the top 7 bits of the mantissa pick one of 128 segments. Within a segment the slope is shared, and it equals the true derivative at the segment's centre. That's the textbook layout of a hardware quadratic interpolator: a small table indexed by the top bits, plus arithmetic on the rest.
The fingerprint
With the slope pinned, each 64-input interval still carried its own offset. Subtracting a smooth quadratic from those offsets left this:
It's a sawtooth, one ulp tall, and the same in every function. The rate at which each ramp climbs matched something specific: the local slope of that function's squared term. For rcp that's −0.273 ulps per interval, and the ramp climbs 0.28. For exp2 it's −0.151, and the ramp climbs 0.15. rsqrt's ramp looks different, but its slope of −0.958 aliases to a small drift, exactly as a floor would make it.
So the squared term is floored to whole ulps on its own before it is added in. The linear term isn't floored with it, and neither is the total. That's also why no polynomial could fit: a floor with its own phase is invisible to least squares.
A long detour, and some bugs of my own
Knowing that the squared term is floored didn't immediately say how. I tried floors with a phase and floors at quarter and half ulps. I tried truncated squarers that drop low partial products, split squarers, squarers keeping only 12 significant bits, and two-step products. All of them came out at the same wall: 99.65% of outputs right, with the rest one step off, at points where the squared term sat within about 0.015 ulp of an integer.
Several of the dead ends were my own bugs. One search stepped a coefficient so coarsely that it could never hit the true value. One extraction subtracted a constant twice. A zsh quirk made one parameter sweep silently test nothing. When a hypothesis fails cleanly, the right move is to check the checker first.
The break: identical margins
What broke it open wasn't a new hypothesis but a coincidence in the output of a sweep. I fitted every segment of rcp separately and printed the best-fit squared coefficient and its error margin. Neighbouring segments often showed identical values to four decimals, for example segments 115 to 120 of rcp all at 143.9636. Identical decimals meant identical integer data.
Two things followed. The squared coefficients are all multiples of 8/2^20, so the coefficient is a small integer n, with a squared term of n·t²/2^17. And the per-interval integer sequences were byte-for-byte identical across functions: exp2's segment 32 matches rcp's segment 116, and rsqrt's segment 100 matches rcp's segment 54. The squared-term circuit is one unit shared by all of them, and its output depends only on n and t.
That turned the squarer from a guess into a measurement. With 109 different values of n in the data, I could solve for the squarer's output T(t) directly. For each t, each n constrains T(t) to an interval, and 109 intervals intersect very tightly:
The squarer rounds t² up to a multiple of 256. With that, the squared term is (n · ⌈t²/256⌉) >> 9, and it matched every n and every t with zero mismatches.
The model for a segment is then:
t = |(x2 >> 6) − 512|
v = c0 + ((m · x2) >> 17) + ((n · ⌈t²/256⌉) >> 9)
result = v & ~3 (in ulps of the segment's exponent)
Here c0 is a whole number of ulps, m has about 18 bits and n has 7 to 8. For each segment and each candidate (m, n), every output pins c0 to an interval, so fitting a segment is a scan over a few thousand m values. All 512 segments of rcp, rsqrt, sqrt and exp2 fit. The resulting code matched fp64's over all 2^32 inputs on the first run, apart from one bug in the wrapper arithmetic.
sin: counting backwards
sin didn't fit, at any segment size. The clue came from its few non-monotonic outputs, places where the output steps the wrong way. For rcp those only happen between inputs 64k + 63 and 64k + 64, at the interval boundaries. For sin they happen one input later, between 64k and 64k + 1. That's what you get if the hardware counts from the other end. It indexes the quarter wave with y = 2^23 − x, which amounts to computing sin as the cosine of the complementary angle. Re-indexed that way, the same interpolator fits.
What remained was scale. Most of sin's range is below 0.5, and some segments cross into a lower binade partway through. The rule turned out to be simple. Each segment works in the ulps of its first, largest output, and truncates to 4 of those even for results that drop into the next binade down. asin confirmed it from the other side. Its values rise within a segment, so results that climb into the next binade keep one extra bit. Those are exactly the 1.25% of asin outputs with 23 significant bits instead of 22.
log2: the hardest one
log2 for inputs in [1, 4) fitted immediately. For larger inputs, and for inputs below 1, it didn't. The result is exponent + log2(1.m), and as the exponent grows it takes up bits the fraction can no longer use:
Every hypothesis about how the hardware reaches the coarser step failed for a while. The outputs below 1 looked like a different, less accurate computation altogether. The resolution was that the datapath cuts the coefficients down to match the output's precision. At level d, where the step is 2^(d+2) units:
- the slope m loses its low d + 2 bits;
- the squared coefficient loses the low d bits of its magnitude;
- c0 absorbs the part of the squared term that was dropped, as evaluated at the segment's edge.
For negative exponents that is level 7, with a step of 2^-15. There the squared term disappears entirely, which is why that path looked like plain linear interpolation. The sum is truncated toward zero, which rounds negative results up. That single rule also covers the region just below 1.0 that fp64's code handled as a special case, including the sign of −0.
What came out
| Tables | Interpolator | |
|---|---|---|
| Data | 4.9 MB loaded from assets/vfpu |
10.5 KB of static const coefficients |
| rcp | 4.2 ns | 3.6 ns |
| rsqrt | 5.4 ns | 3.6 ns |
| exp2 | 6.3 ns | 3.9 ns |
| log2 | 7.1 ns | 5.1 ns |
| sin | 14.6 ns | 5.9 ns |
| asin | 6.0 ns | 3.6 ns |
All are bit-exact over every 32-bit input. vrot needs sine and cosine together, and sharing the argument reduction makes the pair about 10% faster.
What does this say about the chip? It's the design from the literature on hardware function evaluation, for example Piñeiro, Oberman, Muller and Bruguera's minimax quadratic interpolator, or NVIDIA's multifunction interpolator (Oberman and Siu, 2005), which covers much the same set of functions. It has a 128-entry ROM per function, a squarer that sees only the top 10 bits of the offset, and coefficient widths trimmed to what the output needs. Whether the final adder is shared with vdot, as on the Dreamcast, can't be told from outputs alone. The separately truncated products are certainly compatible with it.
Lessons
- An exact oracle plus exhaustive checking beats cleverness. Every idea here was cheap to test against all 2^32 inputs, so wrong ideas died quickly.
- Ask whether it's possible, not how close you can get. Least squares hid the discrete structure. Interval feasibility exposed it.
- Watch for suspicious coincidences. Identical margins across segments was the single most useful observation. It came from output I wasn't looking for.
- Suspect your own tools. Several "impossible" results were bugs in my search code.