diff --git a/pcsx2/FPU.cpp b/pcsx2/FPU.cpp index 55d5e93e83..c69823b56c 100644 --- a/pcsx2/FPU.cpp +++ b/pcsx2/FPU.cpp @@ -418,181 +418,186 @@ static u32 eeGuardedAddSub(u32 a, u32 b, bool issub) return eeRoundToSingle(eeToDouble(a) + eeToDouble(b), true); } -/* Divide two EE singles with exactly one rounding. +/* The EE's divide/square-root unit, digit by digit. - Not through eeToDouble(), the way the adder and the multiplier go: their - results are exact in a double and a quotient is not, so widening and then - narrowing rounds twice. Chopping would forgive that, but the divide unit - rounds to nearest (FPUDivFPCR, scoped in by the callers), where the second - rounding can land a ULP away. + It is not a correctly rounding divider, and no rounding mode makes it one. + The exact quotient lies between two singles and silicon returns one of them, + but which one is decided by where a digit recurrence lands, not by how far the + exact value sits from either end. So the error reaches nearly a whole ULP in + both directions on the same divisor, and two operand pairs agreeing in every + coordinate a rounding rule can see -- branch, divisor, remainder, the exact + fraction -- can still go opposite ways. - So rescale instead of widening. Both operands are forced to exponent 127, so - the division happens between two significands in [1,2) where nothing can - overflow or underflow and the host performs the EE's single rounding; the - exponents are added back onto the quotient afterwards. Scaling by a power of - two leaves a significand alone, so the quotient's significand and its - rounding do not depend on where the operands sat in the range -- only the - reassembled exponent does, and that is integer arithmetic. + It is a radix-2 SRT digit recurrence over the signed digit set {-1, 0, +1}, + partial remainder carried in redundant carry-save form, 24 selections + producing 25 digits, and no rounding step anywhere in it: no round bit, no + sticky, no final correction. SQRT.S is the same recurrence with the root so + far fed back in place of the divisor, and RSQRT.S is SQRT.S followed by DIV.S + with an ordinary 24-bit single in between, which is what silicon does. - The divisor must already be known nonzero: a zero divisor is a flag question - the callers answer first. + The model is the selection function eeSrtDigit() below. - The console does not round to nearest; eeDivideTruncates() below covers the - part of the difference that is settled. + * DIV.S -- 150,994,944 rows. Eighteen divisors swept exhaustively, every + one of the 2^23 numerator significands at each, both branches, including + the seven degenerate divisors that broke every earlier frame. Zero + result-word disagreements, and zero rows where the recurrence lands + outside {T, T+1}. Through the same scorer, truncation gets 52,434,906 of + those rows wrong and round-to-nearest 29,430,553. + * SQRT.S -- 16,777,216 rows, every significand at both exponent parities. + Zero disagreements; truncation 3,994,228 and round-to-nearest 4,395,034. + * A 3000-operand holdout (seed 0x243F6A88, captured before any of this was + modelled): 2000 div and 1000 sqrt, zero disagreements. 229 of the div + rows saturate or flush, which the exhaustive sweeps cannot see: they hold + both exponents at 127. + + Nothing else moves the result. FCR31 has no rounding-mode field on this FPU + -- bit 0 reads 1 whatever is written and bit 1 never sticks -- and the answer + is unchanged when the preceding instruction is varied or the case list + reversed. What that leaves EmuConfig.Cpu.FPUDivFPCR to do is in the block + below eeSqrtBits(). + + This is PS2Float.cpp's Div() and Sqrt() from the proposed PCSX2 soft-float + series (GitHubProUser67), whose only documentation is M. Prabhu and + G. Zyner, "167 MHz radix-8 divide and square root using overlapped + radix-2 stages", DOI 10.1109/ARITH.1995.465363. */ -/* The EE divider's truncation law. - - The divide/square-root unit is not correctly rounded: the exact quotient - lies between two singles and the unit returns one of them, not always the - nearer, and which one it returns depends on the operands rather than on a - rounding mode. This is the part of that choice the captures settle. - - Write the exact division of the two significands (hidden bit restored) as - - lt = ma < mb the branch: does the quotient need a shift - num = ma << (23 + lt) - T = num / mb the truncated 24-bit significand - rem = num - T*mb 0 <= rem < mb - u = mb - rem how far the exact quotient sits below T+1 - - so the unit returns T or T+1 and correct rounding would take T+1 exactly - when 2*rem >= mb. The unit rounds up only when u is small: - - u > cap => the unit truncates cap = 2^22 on A>=B - cap = max(2^23, mb-2^22) on A cap. Not one of them rounds up. - * 48,799,468 further rows of scattered, band and transverse captures: - 6,937,248 distinct divisor significands, divisor exponent fields 110 - through 145. Again not one violation. - - Rows with u <= cap are unsettled -- 27.5% of them truncate as well, and no - model this project has built predicts which -- so they keep the correctly - rounded answer, which is also what both recompilers produce. The implication - runs one way, so applying the law can only turn a wrong row right. - - Only the A>=B half of the cap ever changes a result: on A cap already - implies correct rounding truncates. The branch is kept because it is the law - the captures give, and - EeFpuDivUnitExhaustive.TheAlbHalfOfTheCapCannotChangeAnAnswer holds its - shape so a tightened cap cannot silently become live. - - Over the exhaustive set this takes the interpreter from 19.49% of quotients - off by a ULP to 15.44%, and over the scattered set from 13.98% to 9.51%, at - the price of disagreeing with both recompilers on those rows. - EeRecFpuDivUnitRounding pins the divergence to this class. -*/ -static bool eeDivideTruncates(u32 mb, u32 lt, u32 rem) +/* The partial remainder, in the redundant form the recurrence carries it in. + No step propagates a carry across the width of the operand, which is why the + selector below cannot see the true remainder. */ +struct EeSrtRemainder { - const u32 cap = lt ? ((mb > (3u << 22)) ? mb - (1u << 22) : (1u << 23)) : (1u << 22); - return (mb - rem) > cap; + u32 sum, carry; +}; + +static __fi EeSrtRemainder eeSrtCarrySave(u32 a, u32 b, u32 c) +{ + const u32 u = a ^ b; + const u32 h = (a & b) | (u & c); + return {u ^ c, h << 1}; +} + +/* The selection function: which of -1, 0, +1 the next digit takes. + + It assimilates the redundant remainder only partially -- the carry word is + added in above bit 23 while the low 24 bits of the sum are OR-ed back rather + than added -- so it decides on something other than the remainder, and picks a + different digit from the one an exact comparison would. The thresholds are + +2^23 and -2^24, with the binary point between bits 24 and 25; that asymmetry + is what biases the unit toward truncation. The digit set is redundant, so the + last digit can still be -1 and the result reach T+1 on rows a truncation could + never reach. */ +static __fi s32 eeSrtDigit(EeSrtRemainder r) +{ + constexpr u32 mask = (1u << 24) - 1u; + const s32 estimate = (s32)(((r.sum & ~mask) + r.carry) | (r.sum & mask)); + return (estimate >= (1 << 23)) - (estimate < (s32)(~0u << 24)); +} + +/* On a zero digit the next digit is selected from the un-recompressed pair + while the state advances with the recompressed one, so selection and state + see different splittings of the same value. Drop the distinction and the + model stops reproducing silicon. */ +static __fi EeSrtRemainder eeSrtSelect(EeSrtRemainder cur, EeSrtRemainder next, s32 digit) +{ + const u32 m = 0u - (u32)(digit != 0); + return {(cur.sum & ~m) | (next.sum & m), (cur.carry & ~m) | (next.carry & m)}; +} + +/* The quotient of two 24-bit significands as 25 digits, weights 2^24 down to + 2^0. A positive digit subtracts the divisor as ~divisor with the +1 fed into + the carry word, which the selector then sees -- that is one of the places + the estimate and the state come apart. The value returned is 25 bits when + sma >= smb and 24 bits when it is not; the caller normalises, and the digit + that falls off the bottom there is simply dropped. */ +static u32 eeDivideSignificand(u32 sma, u32 smb) +{ + const u32 divisor = smb << 2; + EeSrtRemainder rem = {sma << 2, 0}; + u32 quotient = 0; + s32 digit = 1; + + for (int i = 0; i < 24; ++i) + { + quotient = (quotient << 1) + (u32)digit; + const u32 addend = (digit > 0) ? ~divisor : ((digit < 0) ? divisor : 0u); + rem.carry += (u32)(digit > 0); + const EeSrtRemainder next = eeSrtCarrySave(rem.sum, rem.carry, addend); + digit = eeSrtDigit(eeSrtSelect(rem, next, digit)); + rem.sum = next.sum << 1; + rem.carry = next.carry << 1; + } + return (quotient << 1) + (u32)digit; } static u32 eeDivide(u32 a, u32 b) { const s32 ea = (s32)((a >> 23) & 0xFF); const s32 eb = (s32)((b >> 23) & 0xFF); + const u32 sign = (a ^ b) & 0x80000000u; if (ea == 0) - return (a ^ b) & 0x80000000; // zero dividend, sign from both operands + return sign; // zero dividend (denormals are zero), sign from both operands - // The exact frame, in integers, so the decision below owes nothing to the - // host's rounding mode. Exponent 255 is an ordinary binade on this FPU, so - // every finite operand reaches here and the hidden bit is always present. + // Exponent 255 is an ordinary binade on this FPU, so every finite operand + // reaches here and the hidden bit is always present. The divisor is already + // known nonzero: that is a flag question the callers answer first. + u32 quotient = eeDivideSignificand(0x800000u | (a & 0x7FFFFFu), 0x800000u | (b & 0x7FFFFFu)); + s32 e = ea - eb + 126; + if (quotient >= (1u << 24)) { - const u32 sma = 0x800000u | (a & 0x7FFFFFu); - const u32 smb = 0x800000u | (b & 0x7FFFFFu); - const u32 lt = (sma < smb) ? 1u : 0u; - const u64 num = (u64)sma << (23 + lt); - const u32 T = (u32)(num / smb); - const u32 rem = (u32)(num - (u64)T * smb); - - if (eeDivideTruncates(smb, lt, rem)) - { - // T is already normalised into [2^23, 2^24), so nothing carries out - // of the significand and the exponent is pure integer. A quotient - // that only reaches the next binade by rounding up therefore does - // not reach it here, as on the console. - const s32 e = ea - eb + 127 - (s32)lt; - const u32 sign = (a ^ b) & 0x80000000u; - if (e > 255) - return sign | 0x7FFFFFFFu; // the EE's maximum, not FLT_MAX - if (e < 1) - return sign; // the EE has no denormals to underflow into - return sign | ((u32)e << 23) | (T - 0x800000u); - } + quotient >>= 1; + ++e; } - FPRreg ma, mb, q; - ma.UL = (a & 0x807FFFFFu) | (127u << 23); - mb.UL = (b & 0x807FFFFFu) | (127u << 23); - q.f = ma.f / mb.f; - - // |q| is in (0.5, 2], so its exponent field carries 126, 127 or -- if the - // rounding pushed it to exactly 2.0 -- 128. - const s32 e = (s32)((q.UL >> 23) & 0xFF) + ea - eb; + // No carry out of the significand is possible below this point -- there is + // no rounding step left that could walk the quotient into the next binade. if (e > 255) - return (q.UL & 0x80000000u) | 0x7FFFFFFFu; + return sign | 0x7FFFFFFFu; // the EE's maximum, not FLT_MAX if (e < 1) - return q.UL & 0x80000000u; // the EE has no denormals to underflow into - return (q.UL & 0x807FFFFFu) | ((u32)e << 23); -} - -/* floor(sqrt(x)) for x < 2^48, exactly. The host sqrt only seeds it: x is - under 53 bits, so it converts to a double without loss and lands within one - of the answer, and the fixup loops run unconditionally, so the result does - not depend on the host's rounding mode. */ -static u32 eeISqrt48(u64 x) -{ - u64 r = (u64)std::sqrt((double)x); - while (r > 0 && r * r > x) - --r; - while ((r + 1) * (r + 1) <= x) - ++r; - return (u32)r; + return sign; // the EE has no denormals to underflow into + return sign | ((u32)e << 23) | (quotient & 0x7FFFFFu); } /* sqrt(|Ft|) as EE bits, including the top binade. - Integer, for the same reason eeDivide is: the square-root unit is the divide - unit, and it misses the correctly rounded answer under the same law. + The same recurrence as eeDivideSignificand(), with the root so far fed back + in place of a fixed divisor: adding digit d at weight w to the root adds + d*(2*root + d*w) to its square, which is what has to leave the partial + remainder, so the addend is rebuilt each step instead of being a constant. + Everything else -- the carry-save state, the selector, the zero-digit quirk, + the absence of any rounding step -- is shared, and so is the evidence: see + the block comment above eeSrtDigit(). - Put the operand's significand where the root is a 24-bit integer. With E the - exponent field and m the significand with its hidden bit, + The radicand is placed so the root is a 24-bit integer. With E the exponent + field, the significand's hidden bit restored, and the result's exponent + (E + 127) / 2 rounded down, the operand's own exponent parity decides how far + to shift: one place when E is odd, two when it is even. The top binade needs + no special case: in integers it is just another odd E. +*/ +static u32 eeSqrtSignificand(u32 m) +{ + EeSrtRemainder rem = {m, 0}; + u32 root = 0; + s32 digit = 1; - k = 23 if E is odd, 24 if E is even (|Ft| = X * 2^(E-150-k)) - X = m << k 2^46 <= X < 2^48 - R = floor(sqrt(X)) 2^23 <= R < 2^24 - rem = X - R*R, u = (R+1)^2 - X = 2R+1-rem + for (int i = 0; i < 24; ++i) + { + const u32 addend_base = root + ((u32)digit << (24 - i)); + root += (u32)digit << (25 - i); + const u32 addend = (digit > 0) ? ~addend_base : ((digit < 0) ? addend_base : 0u); + rem.carry += (u32)(digit > 0); + const EeSrtRemainder next = eeSrtCarrySave(rem.sum, rem.carry, addend); + digit = eeSrtDigit(eeSrtSelect(rem, next, digit)); + rem.sum = next.sum << 1; + rem.carry = next.carry << 1; + } + // The last digit carries weight 2^1, below the root's least significant + // bit, so it only reaches the result by borrowing out of it. + root += (u32)digit << 1; + return (root >> 2) & 0xFFFFFFu; +} - so the unit returns R or R+1, correct rounding takes R+1 exactly when - X > (R+0.5)^2 -- i.e. rem > R, with no ties possible since (R+0.5)^2 is - never an integer -- and E-150-k is even by construction, so the result's - exponent field 150 + (E-150-k)/2 is exact integer arithmetic. - - The truncation law is the one eeDivideTruncates() applies, with the constant - that goes with it: - - u > 2^23 => the unit truncates - - 2^23 is half of sqrt's minimum span 2R+1 >= 2^24 + 1, as 2^22 is half of - div's minimum span mb >= 2^23. Measured by - captures/fpmatrix/divsqrt/scatter/sqgain.c over the two exhaustive console - sweeps (16,777,216 rows, every significand at both exponent parities, - SCPH-90000, FCR0 0x2E40): - - * 10,845,747 rows have u > 2^23. Not one of them rounds up. - * An odd-exponent row rounds up at u = 2^23 exactly, so the bound is - attained. - - It takes sqrt from 26.20% of roots off by a ULP to 11.55%. The rest is the - same unsolved region as div's, and keeps the correctly rounded answer. */ static u32 eeSqrtBits(u32 t) { const u32 E = (t >> 23) & 0xFFu; @@ -601,78 +606,21 @@ static u32 eeSqrtBits(u32 t) // do both recompilers (they take |Ft| first). See // EeRecFpu.SqrtSOfNegativeZeroIsPositiveZero. - const int k = (E & 1u) ? 23 : 24; - const u64 X = (u64)(0x800000u | (t & 0x7FFFFFu)) << k; - const u32 R = eeISqrt48(X); - const u64 rem = X - (u64)R * R; - const u32 u = (u32)(2ull * R + 1ull - rem); - - // The mode the divide unit runs in, honoured here rather than through the - // host FPCR because this arithmetic is integer. sqrt's result is always - // positive, so toward-negative-infinity is the same as toward zero and - // toward-positive-infinity is "up whenever the root is inexact". The - // truncation law sits on top: it can only take the increment away, so under - // toward-positive-infinity it suppresses a round-up the mode asked for. - // That is deliberate -- the law is what the console does, the non-nearest - // modes are a compatibility knob for behaviour it does not have. All four - // modes are pinned by - // EeRecFpuDivUnitRounding.SqrtSHonoursEveryDivideUnitRoundingMode. - bool round_up; - switch (EmuConfig.Cpu.FPUDivFPCR.GetRoundMode()) - { - case FPRoundMode::Nearest: round_up = rem > (u64)R; break; - case FPRoundMode::PositiveInfinity: round_up = rem != 0; break; - default: round_up = false; break; - } - - u32 sig = (round_up && u <= (1u << 23)) ? R + 1u : R; - s32 e = 150 + (((s32)E - 150 - k) / 2); // the numerator is always even - if (sig == 0x1000000u) // rounded out of the binade - { - sig = 0x800000u; - ++e; - } - return ((u32)e << 23) | (sig - 0x800000u); + const u32 m = (0x800000u | (t & 0x7FFFFFu)) << ((E & 1u) ? 1 : 2); + return (((E + 127u) >> 1) << 23) | (eeSqrtSignificand(m) & 0x7FFFFFu); } -/* The EE's divide/square-root unit rounds to nearest even when the rest of the - FPU is chopping toward zero, which is why PCSX2 carries a second control - register, FPUDivFPCR, whose only difference from FPUFPCR is the rounding - mode. Both recompilers swap the host rounding mode around DIV/SQRT/RSQRT +/* Nothing in the interpreter swaps EmuConfig.Cpu.FPUDivFPCR into the host FPCR + any more: DIV.S, SQRT.S and RSQRT.S are integer arithmetic now, and no + rounding mode reaches a digit recurrence. That register is PCSX2's surrogate + for the divide unit rounding to nearest while the rest of the FPU chops, and + both recompilers still swap it in, because they run these ops on host singles (arm64 recDIV_S_xmm / recSQRT_S_xmm / recRSQRT_S_xmm in iFPU-arm64.cpp and - the DOUBLE:: twins in iFPUd-arm64.cpp; x86 iFPU.cpp / iFPUd.cpp do the same - with xLDMXCSR). The interpreter never did, so those three ops came out one - ULP low against both recompilers and the console whenever a game is in the - default chop mode. - - Gated the way the emitters gate it: where the two registers already agree - there is nothing to swap. + the DOUBLE:: twins in iFPUd-arm64.cpp; x86 iFPU.cpp / iFPUd.cpp with + xLDMXCSR): 1.0 rsqrt 1.5 is 0x3F5105EB on the console and 0x3F5105EC without + the swap. So the two engines part company on every operand silicon is not + correctly rounded on, which EeRecFpuDivUnitRounding and EeRecFpuRsqrt pin. */ -class ScopedDivRoundMode -{ -public: - __fi ScopedDivRoundMode() - : m_swap(EmuConfig.Cpu.FPUFPCR.bitmask != EmuConfig.Cpu.FPUDivFPCR.bitmask) - { - if (m_swap) - { - m_prev = FPControlRegister::GetCurrent(); - FPControlRegister::SetCurrent(EmuConfig.Cpu.FPUDivFPCR); - } - } - __fi ~ScopedDivRoundMode() - { - if (m_swap) - FPControlRegister::SetCurrent(m_prev); - } - - ScopedDivRoundMode(const ScopedDivRoundMode&) = delete; - ScopedDivRoundMode& operator=(const ScopedDivRoundMode&) = delete; - -private: - FPControlRegister m_prev; - bool m_swap; -}; void ABS_S() { _FdValUl_ = _FsValUl_ & 0x7fffffff; @@ -754,7 +702,6 @@ void CVT_W() { } void DIV_S() { - const ScopedDivRoundMode div_round; if (checkDivideByZero( _FdValUl_, _FtValUl_, _FsValUl_, FPUflagD | FPUflagSD, FPUflagI | FPUflagSI)) return; _FdValUl_ = eeDivide( _FsValUl_, _FtValUl_ ); } @@ -1043,7 +990,6 @@ void NEG_S() { } void RSQRT_S() { - const ScopedDivRoundMode div_round; clearFPUFlags(FPUflagD | FPUflagI); if ( ( _FtValUl_ & 0x7F800000 ) == 0 ) { // Ft is zero (Denormals are Zero) @@ -1066,23 +1012,20 @@ void RSQRT_S() { // gives 0x3F5105EB. // - // Neither operand is clamped any more, and it has to be both: with the clamp - // left on the sqrt alone, rsqrt(2^128, 2^128) came out right only because - // the two clamps cancelled. Unclamping both fixed that row and 13 others, - // taking RSQRT.S from 17/32 to 31/32 against the console. + // Neither operand is clamped any more, which is all-or-nothing by design: + // unclamping only the sqrt used to break rsqrt(2^128, 2^128), which came out + // right solely because its two clamps cancelled. Unclamping both fixes that + // row and 13 others. Scored against the console over the corpus, RSQRT.S + // went 17/32 to 31/32. // - // rsqrt(EEMAX, EEMAX) is still 1 ULP out, and it is not the two-step - // rounding it was filed as: silicon composes the two steps exactly as - // below, with a plain 24-bit single in between, and it is the divide/ - // square-root unit itself that is not correctly rounded. This computes the - // correctly-rounded answer; ee_fpu_divunit_console_tests.cpp has the - // capture and how far silicon sits from it. + // The composition itself comes from a dedicated console capture: 2231 + // operand pairs over two probes, each pair run as sqrt.s, rsqrt.s and + // div.s, and rsqrt.s equals div.s(Fs, sqrt.s(Ft)) on every row, with a + // plain 24-bit single in between. See ee_fpu_divunit_console_tests.cpp. _FdValUl_ = eeDivide( _FsValUl_, eeSqrtBits( _FtValUl_ ) ); } void SQRT_S() { - // No ScopedDivRoundMode: eeSqrtBits() reads FPUDivFPCR's rounding mode - // itself. DIV.S and RSQRT.S still need it for eeDivide()'s host division. clearFPUFlags(FPUflagI | FPUflagD); // Invalid-operation keys off the SIGN BIT ALONE. -0 and the negative diff --git a/tests/ctest/core/recompilers/ee_fpu_divunit_console_tests.cpp b/tests/ctest/core/recompilers/ee_fpu_divunit_console_tests.cpp index 3026a5e1c5..b8c4551178 100644 --- a/tests/ctest/core/recompilers/ee_fpu_divunit_console_tests.cpp +++ b/tests/ctest/core/recompilers/ee_fpu_divunit_console_tests.cpp @@ -46,53 +46,30 @@ // A reciprocal-then-multiply model scores WORSE than plain // correctly-rounded (69% against 84%), so that is not the shape either. // -// Part of that is now modelled. eeDivide() and eeSqrtBits() in FPU.cpp apply -// the one rule the captures settle, in the shape eeMulRound already uses for -// the multiplier's deficit: +// All of it is now modelled. eeDivide() and eeSqrtBits() in FPU.cpp run the +// unit's own radix-2 SRT digit recurrence rather than a rounded host operation, +// so this table is no longer a scoreboard with a residue -- every cell of it +// matches. What the columns did as the model went in: // -// u > cap => the unit truncates +// op correctly truncation the +// rounded law recurrence +// sqrt.s 66/87 -> 87/87 -> 87/87 +// div.s 66/87 -> 70/87 -> 87/87 +// rsqrt.s 61/87 -> 78/87 -> 87/87 // -// where u is how far the exact result sits below the upper candidate and -// cap is 2^22 (div, A>=B), max(2^23, mb-2^22) (div, A 87/87 +21 0 -// div.s 66/87 -> 70/87 +4 0 -// rsqrt.s 61/87 -> 78/87 +19 2 -// -// The two rsqrt losses were rows where the old code was right by cancellation: -// its square root came back one ULP high (correctly rounded, where silicon -// truncates) and its division then rounded up to the same word silicon reached -// by truncating a smaller divisor. Modelling the square root removed one half -// of that pair; the other half, a division inside the region the law does not -// settle, is still wrong. Both rows now have a sqrt.s column that matches the -// console. -// -// rsqrt.s 3895AEC3, 4938608B was 33B06019 (console), now 33B0601A -// rsqrt.s 43CD0CEB, 365AF7C1 was 485DB675 (console), now 485DB676 -// -// Not modelled: the region u <= cap, where 27.5% of div rows truncate as well -// and no coordinate system this project has built predicts which. See -// FINDINGS-div-round11-exhaustive.md; the missing coordinate is a carry -// propagation distance in the trial product, known in shape and not in form. -// -// Corpus case 220, rsqrt EEMAX, EEMAX, is still the last result-axis row where -// the interpreter misses the console, and is still one sample of that unmodelled -// region -- its division is div.s 7FFFFFFF, 5FB504F3 with u = 3,571,369 against -// an A>=B cap of 4,194,304, so the law does not reach it and the correctly -// rounded answer stands where silicon truncates. It is the only row in all 1147 -// corpus cases that can see the unit's approximation at all, which is why the -// whole corpus is byte-identical before and after this change. -// -// This file is still interim. The tripwire at the bottom went from 68 failing -// assertions to 26; enabling it means someone closed the u <= cap region. Plan -// and acceptance criteria: WORKORDER-divsqrt-model.md in the notes tree. +// The `ieee_*` column stays as the other engine's column: both recompilers +// still take the host's correctly-rounded fdiv/fsqrt, so the divergence between +// the two engines is exactly the 68 cells where the console and correct +// rounding differ. TheFastPathStaysCorrectlyRoundedAndSaysSoHere asserts both +// halves of that. // // Evidence archive: captures/fpmatrix/rsqprobe.c, hw-rsq-run{1,2,3}.bin and -// PROBE-divsqrt-rounding.md in the notes tree. +// PROBE-divsqrt-rounding.md in the notes tree; the recurrence and the +// 167,772,160 rows it was scored on are documented at eeSrtDigit() in FPU.cpp. #include "harness/EeRecTestHarness.h" #include "harness/MipsEncode.h" @@ -288,42 +265,16 @@ TEST(EeFpuDivUnitConsole, RsqrtIsSqrtThenDivide) } // --------------------------------------------------------------------------- -// 2. How much of the unit the interpreter now reproduces, and where what is -// left over lives. +// 2. The interpreter reproduces the console on every cell of the table. // -// The per-op counts are drift guards: they will move when someone models -// more of the unit. The assertion inside the loop is the one that matters: -// every row the interpreter still misses must be a division whose `u` lies -// at or below the cap, inside the region the truncation law does not -// settle. A miss above the cap refutes the law rather than widening it. +// This was a disabled tripwire while the model was partial; it is the +// acceptance test now. The per-op counts are drift guards, and 193 of the +// 261 cells cannot tell a recurrence from a correctly rounded divider, which +// is what the followed_silicon count below is for. // --------------------------------------------------------------------------- -namespace { - -// The frame the model is stated in, recomputed rather than shared with FPU.cpp -// so this file cannot inherit the arithmetic it is checking. -struct DivFrame +TEST(EeFpuDivUnitConsole, InterpMatchesTheConsoleOnEveryRow) { - u32 T, rem, u, cap; - int lt; -}; - -DivFrame Frame(u32 a, u32 b) -{ - const u32 ma = 0x800000u | (a & 0x7FFFFFu); - const u32 mb = 0x800000u | (b & 0x7FFFFFu); - const int lt = ma < mb ? 1 : 0; - const u64 num = static_cast(ma) << (23 + lt); - const u32 T = static_cast(num / mb); - const u32 rem = static_cast(num - static_cast(T) * mb); - const u32 cap = lt ? std::max(1u << 23, mb - (1u << 22)) : (1u << 22); - return { T, rem, mb - rem, cap, lt }; -} - -} // namespace - -TEST(EeFpuDivUnitConsole, WhatIsLeftOverIsAllInsideTheUnsettledRegion) -{ - int match[3] = {}, total[3] = {}, fixed_by_the_model = 0; + int match[3] = {}, total[3] = {}, followed_silicon = 0; for (const ConsoleRow& r : kRows) { @@ -331,46 +282,23 @@ TEST(EeFpuDivUnitConsole, WhatIsLeftOverIsAllInsideTheUnsettledRegion) { const u32 got = RunInterp(r, op), con = Con(r, op), ieee = Ieee(r, op); ++total[op]; - if (got == con) - { - ++match[op]; - if (con != ieee) - ++fixed_by_the_model; // silicon disagrees with correct - // rounding here and we followed silicon - continue; - } - - // SQRT.S has no unmodelled rows left on this table at all, so a - // square-root miss is a regression and not a known gap. - ASSERT_NE(op, OP_SQRT) - << "sqrt.s fs=" << std::hex << r.fs << " ft=" << r.ft - << " got " << got << " console " << con - << " -- the sqrt column was exact on all 87 rows"; - - // RSQRT.S is sqrt-then-divide, so its division is by the root the - // interpreter just produced, not by Ft. - const u32 divisor = (op == OP_DIV) ? r.ft : RunInterp(r, OP_SQRT); - const DivFrame f = Frame(r.fs, divisor); - EXPECT_LE(f.u, f.cap) + match[op] += (got == con); + followed_silicon += (got == con && con != ieee); + EXPECT_EQ(got, con) << OpName(op) << " fs=" << std::hex << r.fs << " ft=" << r.ft - << ": missed with u=" << std::dec << f.u << " ABOVE cap=" << f.cap - << ", which the truncation law says cannot happen"; + << " (" << r.what << "), correctly rounded would be " << ieee; } } EXPECT_EQ(match[OP_SQRT], 87) << "sqrt.s"; - EXPECT_EQ(match[OP_DIV], 70) << "div.s"; - EXPECT_EQ(match[OP_RSQRT], 78) << "rsqrt.s"; + EXPECT_EQ(match[OP_DIV], 87) << "div.s"; + EXPECT_EQ(match[OP_RSQRT], 87) << "rsqrt.s"; EXPECT_EQ(total[OP_SQRT] + total[OP_DIV] + total[OP_RSQRT], 261); - // Anti-vacuity. If the model were inert the interpreter would just be the - // correctly-rounded column again and this counter would be 0, so every - // number above would be describing nothing. - // 21 sqrt + 4 div + 19 rsqrt; the net gain on the table is 42 because two - // rsqrt rows that used to be right by cancellation are not any more. - EXPECT_EQ(fixed_by_the_model, 44) - << "rows where silicon and correct rounding differ AND the interpreter " - "followed silicon -- 0 means the truncation law stopped firing"; + EXPECT_EQ(followed_silicon, 68) + << "cells where silicon and correct rounding differ AND the interpreter " + "took silicon's side -- 0 means the recurrence stopped firing and " + "this test is passing on the rows that cannot tell the two apart"; } // --------------------------------------------------------------------------- @@ -444,32 +372,12 @@ TEST(EeFpuDivUnitConsole, TheFastPathStaysCorrectlyRoundedAndSaysSoHere) } EXPECT_EQ(checked, 249); - // 47 of the 249 cells are rows where the interpreter took the truncation - // law and the fast path did not; if this reaches 0 the interpreter's model - // has been lost. - EXPECT_EQ(diverged, 47) + // The divergence is the deliberate part, so it is asserted rather than + // tolerated: 67 of the 249 cells are rows where silicon is not correctly + // rounded, the interpreter follows silicon and the fast path does not. If + // this reaches 0 the interpreter's model has been lost; if it grows, an + // emitter has drifted off the host's own rounding. + EXPECT_EQ(diverged, 67) << "the interpreter is supposed to leave the fast path behind on exactly " - "the rows the truncation law settles"; -} - -// --------------------------------------------------------------------------- -// 5. The tripwire. It fails on 26 assertions today, down from 68 before the -// truncation law landed -- all 26 are divisions inside u <= cap, which test 2 -// asserts one at a time. Enabling this means someone closed that region. -// -// Run it with --gtest_also_run_disabled_tests before graduating it: a -// disabled test that has quietly gone vacuous graduates just as easily as -// one that was fixed. -// --------------------------------------------------------------------------- -TEST(EeFpuDivUnitConsole, DISABLED_InterpMatchesConsoleOnEveryRow) -{ - for (const ConsoleRow& r : kRows) - { - for (Op op : { OP_SQRT, OP_DIV, OP_RSQRT }) - { - EXPECT_EQ(RunInterp(r, op), Con(r, op)) - << OpName(op) << " fs=" << std::hex << r.fs << " ft=" << r.ft - << " (" << r.what << ")"; - } - } + "the rows where the console is not correctly rounded"; } diff --git a/tests/ctest/core/recompilers/ee_fpu_divunit_exhaustive_tests.cpp b/tests/ctest/core/recompilers/ee_fpu_divunit_exhaustive_tests.cpp index cb2f6a884e..3835bd3758 100644 --- a/tests/ctest/core/recompilers/ee_fpu_divunit_exhaustive_tests.cpp +++ b/tests/ctest/core/recompilers/ee_fpu_divunit_exhaustive_tests.cpp @@ -77,34 +77,29 @@ // an upper envelope -- not attained on two divisors, and not monotone in mb // (TheAlbCapIsAnEnvelopeNotTheLaw). // -// 6. What shipped is finding 3 and nothing else. FPU.cpp's -// eeDivide() and eeSqrtBits() now truncate when u exceeds the cap, and keep -// the correctly rounded answer everywhere else, because that implication is -// the only part of the unit with zero measured exceptions. The rows in -// kLawRows are its witnesses, re-measured individually on silicon by wit3.c -// (FCR0 00002E40, FCR31 cleared per op, 32-nop spacers, every result read -// twice, no read disagreed with itself) after being found in the bulk -// captures. Three groups, all needed: +// 6. What shipped. FPU.cpp's eeDivide() and eeSqrtBits() run the unit's own +// radix-2 SRT digit recurrence (see eeSrtDigit there), which reproduces +// every row of every capture including all of the above, so the interpreter +// matches the console on all 21 rows of kRows. Facts 1 through 5 stay as +// bounds on a future fast path rather than on a future model. // -// * six div rows and four sqrt rows where the law changes the answer -- -// without them the model is unpinned and could quietly stop firing; -// * two sqrt rows where u is exactly 2^23 and silicon still rounds up -- -// without them `u > 2^23 => DOWN` would pass vacuously the moment the -// unit got more conservative; -// * two sqrt rows where the law is silent and silicon truncates anyway -- -// the residual, kept so "still wrong here" has names. +// kLawRows are the witnesses of the one-way law the recurrence subsumed -- +// u above the cap implies truncation -- re-measured individually on silicon +// by wit3.c (FCR0 00002E40, FCR31 cleared per op, 32-nop spacers, every +// result read twice, no read disagreed with itself) after being found in the +// bulk captures. Four groups: rows where the law changes the answer on +// DIV.S, the same on SQRT.S, rows where its bound is exactly attained, and +// rows where it stays silent and silicon truncates anyway. That last group +// was the residual and now carries a console expectation like every other +// row. // // The tests below also state that the A cap already implies 2*rem < mb and correct rounding says DOWN by // itself. Verified over all eighteen // exhaustive divisors -- 18,878,960 A=B half does all the -// work: 6,118,759 rows of the 150,994,944 change, and all of them improve. -// -// The interpreter expectations in kRows below are still the ieee column, -// because every one of those rows sits inside u <= cap where nothing is -// settled. The tripwire at the bottom is the acceptance test for the rest. +// them a row correct rounding got right anyway. The A>=B half did all the +// work: 6,118,759 rows of the 150,994,944 changed, and all of them improved. #include "harness/EeRecTestHarness.h" #include "harness/MipsEncode.h" @@ -403,17 +398,17 @@ TEST(EeFpuDivUnitExhaustive, TheErrorSpansNearlyAWholeUlpBothWays) // returns T+1. const Frame hi = Decode(0x3F800000u, 0x3F800001u); EXPECT_EQ(hi.rem, 2u); - EXPECT_EQ(RunDiv(0x3F800000u, 0x3F800001u), 0x3F7FFFFEu) - << "this tree rounds to nearest, which is T here"; EXPECT_EQ(0x3F7FFFFFu, hi.down + 1u) << "the console value is T+1"; + EXPECT_EQ(RunDiv(0x3F800000u, 0x3F800001u), 0x3F7FFFFFu) + << "round-to-nearest gives T here, and this tree does not round"; // and on the SAME divisor, a row where the exact quotient is 0.96 of the // way to T+1 and silicon returns T. const Frame lo = Decode(0x3F852B38u, 0x3F800001u); EXPECT_EQ(lo.u, 338743u); EXPECT_EQ(0x3F852B36u, lo.down) << "the console value is T"; - EXPECT_EQ(RunDiv(0x3F852B38u, 0x3F800001u), 0x3F852B37u) - << "this tree rounds to nearest, which is T+1 here"; + EXPECT_EQ(RunDiv(0x3F852B38u, 0x3F800001u), 0x3F852B36u) + << "round-to-nearest gives T+1 here, and this tree does not round"; // The signed error in ULP: returning T costs -rem/mb, returning T+1 gains // +u/mb. Round-to-nearest is confined to [-1/2, +1/2] and a directed mode @@ -703,11 +698,12 @@ TEST(EeFpuDivUnitExhaustive, InterpMatchesConsoleWhereTheUnitIsExact) } // --------------------------------------------------------------------------- -// 7. Where they disagree, the miss is one ULP and this tree sits on the -// correctly-rounded side. Reaching the console value here is not a failure, -// it is the model landing -- and then the tripwire below is what to enable. +// 7. Where they disagree, the miss is one ULP and this tree takes the console's +// side. These 13 rows used to be out of reach and are the whole of what the +// digit recurrence bought on this table, so a build that went back to +// rounding fails here as well as on test 8. // --------------------------------------------------------------------------- -TEST(EeFpuDivUnitExhaustive, SiliconIsOneUlpOffOnTheseRows) +TEST(EeFpuDivUnitExhaustive, SiliconIsOneUlpOffAndTheInterpreterFollowsIt) { int off = 0; for (const DivRow& r : kRows) @@ -717,24 +713,31 @@ TEST(EeFpuDivUnitExhaustive, SiliconIsOneUlpOffOnTheseRows) ++off; const u32 a = r.con_div & 0x7FFFFFFFu, b = r.ieee_div & 0x7FFFFFFFu; EXPECT_EQ(a > b ? a - b : b - a, 1u) << std::hex << "fs=" << r.fs << " ft=" << r.ft; - EXPECT_EQ(RunDiv(r.fs, r.ft), r.ieee_div) + EXPECT_EQ(RunDiv(r.fs, r.ft), r.con_div) << std::hex << "fs=" << r.fs << " ft=" << r.ft - << ": this tree computes the correctly-rounded value"; + << ": this tree returned the correctly-rounded value " << r.ieee_div + << " where silicon does not round"; } EXPECT_EQ(off, 13); } // --------------------------------------------------------------------------- -// 7b. The law that shipped, and the three things it needs pinned: that it -// fires, that its bound is attained, and that where it stays silent the -// tree is still on the correctly-rounded side. +// 7b. The cap law's own witnesses, kept after the recurrence subsumed it. +// +// Each row still asserts where it sits relative to the cap -- above it, +// exactly on it, below it -- and then asserts the console value, which the +// interpreter now reaches on all four groups rather than on three of them. +// Two reasons to keep the geometry: it is the measurement that says the +// bound is attained and not a free inequality, and any fast path that +// answers `u > cap` rows without running the digits has to agree with the +// recurrence exactly here, including on the rows one unit below the bound. // // The frame is recomputed from the operand bits for every row, so the // annotated `u` cannot drift away from the arithmetic, and a mistyped // operand shows up as a frame mismatch rather than as a mysterious value // failure. // --------------------------------------------------------------------------- -TEST(EeFpuDivUnitExhaustive, TheTruncationLawFiresAndTheInterpreterFollowsSilicon) +TEST(EeFpuDivUnitExhaustive, TheCapWitnessesAllReachTheConsole) { int fires_div = 0, fires_sqrt = 0, tight = 0, silent = 0; @@ -755,7 +758,7 @@ TEST(EeFpuDivUnitExhaustive, TheTruncationLawFiresAndTheInterpreterFollowsSilico EXPECT_EQ(r.console, f.down) << "silicon must be the TRUNCATED candidate"; EXPECT_EQ(r.ieee, f.down + 1u); EXPECT_EQ(RunDiv(r.fs, r.ft), r.console) - << "the interpreter stopped applying the truncation law to DIV.S"; + << "the interpreter's DIV.S stopped reaching silicon above the cap"; } else { @@ -770,7 +773,7 @@ TEST(EeFpuDivUnitExhaustive, TheTruncationLawFiresAndTheInterpreterFollowsSilico << "correct rounding already says R, so this row cannot show the law"; EXPECT_NE(r.console, r.ieee); EXPECT_EQ(RunSqrt(r.ft), r.console) - << "the interpreter stopped applying the truncation law to SQRT.S"; + << "the interpreter's SQRT.S stopped reaching silicon above the bound"; } else if (r.kind == LAW_SQRT_TIGHT) { @@ -790,11 +793,11 @@ TEST(EeFpuDivUnitExhaustive, TheTruncationLawFiresAndTheInterpreterFollowsSilico EXPECT_LE(f.u, 1u << 23) << "this row exists BELOW the bound"; EXPECT_TRUE(f.ieee_rounds_up); EXPECT_NE(r.console, r.ieee) << "silicon truncates here even though the " - "law is silent -- that is the residual"; - EXPECT_EQ(RunSqrt(r.ft), r.ieee) - << "where the law says nothing the tree stays correctly rounded; " - "reaching the console value means someone modelled more of the " - "unit, and this row should move to LAW_SQRT_FIRES"; + "cap law is silent -- these two rows were " + "the residual it could not reach"; + EXPECT_EQ(RunSqrt(r.ft), r.console) + << "the recurrence closed exactly this group; the correctly rounded " + "value here is " << r.ieee; } } } @@ -847,10 +850,14 @@ TEST(EeFpuDivUnitExhaustive, TheAlbHalfOfTheCapCannotChangeAnAnswer) } // --------------------------------------------------------------------------- -// 8. The tripwire. Enabling it means someone modelled the unit. Any model that -// passes it has to satisfy tests 1-4 above, which is the point of them. +// 8. The acceptance test, disabled for as long as the unit was unmodelled: all +// three ops, all 21 rows, against silicon. It was run with +// --gtest_also_run_disabled_tests before being graduated, because a disabled +// test that has quietly gone vacuous graduates just as easily as one that +// was fixed. Tests 1 and 2 above assert that these rows disagree with every +// rounding rule, so passing here takes reproducing silicon. // --------------------------------------------------------------------------- -TEST(EeFpuDivUnitExhaustive, DISABLED_InterpMatchesConsoleOnEveryRow) +TEST(EeFpuDivUnitExhaustive, InterpMatchesConsoleOnEveryRow) { for (const DivRow& r : kRows) { diff --git a/tests/ctest/core/recompilers/ee_fpu_overflow_console_conformance_tests.cpp b/tests/ctest/core/recompilers/ee_fpu_overflow_console_conformance_tests.cpp index 51ad26edca..897fb07ce4 100644 --- a/tests/ctest/core/recompilers/ee_fpu_overflow_console_conformance_tests.cpp +++ b/tests/ctest/core/recompilers/ee_fpu_overflow_console_conformance_tests.cpp @@ -503,19 +503,28 @@ TEST(EeFpuOverflowConsole, SqrtMatchesConsoleOnEveryCapturedOperand) } // --------------------------------------------------------------------------- -// The same property over the whole exponent-255 class rather than the three -// patterns the capture happens to contain. As host bit patterns those words -// are infinities, quiet NaNs and signalling NaNs; to the EE they are all large -// finite floats, so the pool carries every shape. The signalling ones are the -// half a host-NaN-aware implementation gets wrong -- see recSQRT_S_xmm -// (iFPU-arm64.cpp) for why the clamp they replaced had to be an integer Umin -// rather than an Fminnm. +// The same property as above, over the whole exponent-255 class rather than the +// three patterns the capture happens to contain. // -// The wanted values are correctly-rounded square roots computed by exact -// integer arithmetic (math.isqrt on the significand, round-to-nearest-even, -// the divide unit's mode) rather than by a host float. The six marked `true` -// were read off silicon; the other three are computed only, for class -// coverage. +// The class splits on an axis the capture cannot see: as host bit patterns, +// exponent-255 words are infinities, quiet NaNs and signalling NaNs, while to +// the EE they are all large finite floats. The old arm64 clamp had to be an +// integer Umin rather than an Fminnm because of that split -- Fminnm prefers +// the number only against a quiet NaN, and a signalling operand comes back +// merely quieted, so half the mantissa space (4194303 of the 8388608 positive +// patterns) would have passed through a clamp that was supposed to catch it. +// Testing the exponent field, as both engines now do, never asks the host what +// kind of NaN it thinks it is holding; the pool below covers every shape either +// way. +// +// Expected values are correctly-rounded square roots computed by exact integer +// arithmetic (math.isqrt on the significand, round-to-nearest-even) rather than +// by a host float, so they cannot inherit the behaviour under test. The divide +// unit does not round in general -- see eeSrtDigit in FPU.cpp -- but correct +// rounding is what it returns on these operands: the model was checked against +// the six the capture does witness, marked `true` below, and agreed on all six, +// including the exponent-254 control. The other three rows are computed, not +// read off silicon, and are here for class coverage. // --------------------------------------------------------------------------- TEST(EeFpuOverflowConsole, SqrtMatchesConsoleOnEveryExponent255Operand) { diff --git a/tests/ctest/core/recompilers/ee_rec_fpu_divunit_rounding_tests.cpp b/tests/ctest/core/recompilers/ee_rec_fpu_divunit_rounding_tests.cpp index eae0ba519a..0bb9246500 100644 --- a/tests/ctest/core/recompilers/ee_rec_fpu_divunit_rounding_tests.cpp +++ b/tests/ctest/core/recompilers/ee_rec_fpu_divunit_rounding_tests.cpp @@ -1,29 +1,35 @@ // SPDX-FileCopyrightText: 2026 yaps2 Dev Team // SPDX-License-Identifier: GPL-3.0+ // -// The EE's divide/square-root unit rounds to nearest even when the rest of the -// FPU is chopping toward zero. Both recompilers model that by swapping -// EmuConfig.Cpu.FPUDivFPCR in around the three ops that unit owns -- DIV.S, -// SQRT.S and RSQRT.S. The interpreter did not, so it truncated where the -// recompilers and the console round; ScopedDivRoundMode in pcsx2/FPU.cpp has -// the mechanism. +// The EE's divide/square-root unit is not correctly rounded, and the two +// engines in this tree part company over that on purpose. // -// Two things kept that off the differential. Every pre-existing DIV.S case used -// an exactly-representable ratio -- 20/4, 6/-2, with comments saying "no -// rounding divergence" -- and those are the operands at which the two rounding -// modes agree. And the one test that did probe SQRT.S's rounding mode -// (EeRecFpu.SqrtSRoundsToNearestUnderChopFpcr) asserted the JIT result alone, -// because the old harness ran both engines under a nearest-rounding host FPCR -// where the swap is a no-op; that stopped being true when the harness moved to -// the production environment, and the test kept passing either way because it -// never looked at the interpreter. +// The interpreter runs the unit's own radix-2 SRT digit recurrence (FPU.cpp, +// eeSrtDigit and below), which reproduces silicon bit for bit on every capture +// this project has taken. Both recompilers take the host's fdiv/fsqrt under +// EmuConfig.Cpu.FPUDivFPCR -- FPUFPCR with round-to-nearest, swapped in around +// the three ops the unit owns; FPU.cpp names the emitters -- which makes them +// the correctly rounded engine. // -// So this file covers DIV.S and SQRT.S at inexact operands, randomized, in the -// FP environment a game runs in. RSQRT.S is covered the same way in -// ee_rec_fpu_rsqrt_tests.cpp. +// So this file is the class-level regression test for the shape of that +// divergence. Differing is not enough: the engines must differ by exactly one +// ULP, only on the ops the divide unit owns, and only in the direction silicon +// errs -- never above correct rounding on SQRT.S or on DIV.S's A>=B branch, +// either way on DIV.S's A=B div rows, against 3,229,727 +// above and 7,197,471 below in the 78,087,028 A cap), only by one ULP, and only with the -// interpreter on the closer-to-zero side. -// // The predicates are recomputed here rather than exported from FPU.cpp: a // differential that imports the implementation's arithmetic cannot catch the // implementation's arithmetic being wrong. @@ -153,29 +155,21 @@ static bool BothNormalOperands(u32 fs, u32 ft) return ((fs >> 23) & 0xFFu) != 0 && ((ft >> 23) & 0xFFu) != 0; } -static bool DivideTruncates(u32 fs, u32 ft) +// Which branch of the recurrence a division takes. The two branches err +// differently and the tests below hold them to different shapes. +static bool DivideShiftsTheNumerator(u32 fs, u32 ft) { - const u32 ma = 0x800000u | (fs & 0x7FFFFFu); - const u32 mb = 0x800000u | (ft & 0x7FFFFFu); - const int lt = ma < mb ? 1 : 0; - const u64 num = static_cast(ma) << (23 + lt); - const u32 rem = static_cast(num % mb); - const u32 cap = lt ? std::max(1u << 23, mb - (1u << 22)) : (1u << 22); - return (mb - rem) > cap; + return (0x800000u | (fs & 0x7FFFFFu)) < (0x800000u | (ft & 0x7FFFFFu)); } -static bool SqrtTruncates(u32 ft) +// The two candidates a digit recurrence can land on are adjacent words, so +// "one ULP apart" is "one apart as magnitudes" -- true across a binade boundary +// as well, since the float encoding is monotone in the magnitude. +static bool IsOneUlpApart(u32 interp, u32 jit) { - const u32 E = (ft >> 23) & 0xFFu; - if (E == 0) - return false; - const u64 X = static_cast(0x800000u | (ft & 0x7FFFFFu)) << ((E & 1u) ? 23 : 24); - u64 R = static_cast(std::sqrt(static_cast(X))); - while (R > 0 && R * R > X) - --R; - while ((R + 1) * (R + 1) <= X) - ++R; - return (2 * R + 1 - (X - R * R)) > (1u << 23); + const u32 a = interp & 0x7FFFFFFFu, b = jit & 0x7FFFFFFFu; + return (interp & 0x80000000u) == (jit & 0x80000000u) && + (a > b ? a - b : b - a) == 1u; } // The interpreter's word is the JIT's with one unit taken off the magnitude. @@ -185,11 +179,11 @@ static bool IsOneUlpTowardZero(u32 interp, u32 jit) interp == ((jit & 0x80000000u) | ((jit & 0x7FFFFFFFu) - 1u)); } -TEST(EeRecFpuDivUnitRounding, DivSMatchesInterpExceptWhereTheTruncationLawFires) +TEST(EeRecFpuDivUnitRounding, DivSDivergesFromTheFastPathByOneUlpAndOnlyThat) { RequireDistinctDivideRoundingMode(); Lcg r{0xD1F5D1F5A5A5A5A5ull}; - int checked = 0, tier_gaps = 0, law_gaps = 0; + int checked = 0, tier_gaps = 0, gaps = 0, alb_low = 0, alb_high = 0; for (u32 iter = 0; iter < 3000; ++iter) { const u32 fsBits = fuzzOperand(r); @@ -230,15 +224,22 @@ TEST(EeRecFpuDivUnitRounding, DivSMatchesInterpExceptWhereTheTruncationLawFires) } else if (res[0] != res[1]) { - ++law_gaps; - EXPECT_TRUE(BothNormalOperands(fsBits, ftBits) && - DivideTruncates(fsBits, ftBits)) - << "the engines parted company where the truncation law does NOT " - "fire -- that is a plain quotient disagreement, not the " - "modelled one"; - EXPECT_TRUE(IsOneUlpTowardZero(res[0], res[1])) - << "the interpreter's model can only ever take the LOWER of the two " - "candidates; interp=" << std::hex << res[0] << " jit=" << res[1]; + ++gaps; + EXPECT_TRUE(BothNormalOperands(fsBits, ftBits)) + << "the engines parted company on an operand pair the divide unit " + "never sees the digits of -- a zero or denormal operand is a " + "flag question both engines answer the same way"; + EXPECT_TRUE(IsOneUlpApart(res[0], res[1])) + << "the recurrence can only ever land on one of the two candidates " + "the correctly rounded answer sits between; interp=" + << std::hex << res[0] << " jit=" << res[1]; + if (DivideShiftsTheNumerator(fsBits, ftBits)) + ((res[0] & 0x7FFFFFFFu) < (res[1] & 0x7FFFFFFFu) ? alb_low : alb_high)++; + else + EXPECT_TRUE(IsOneUlpTowardZero(res[0], res[1])) + << "on the A>=B branch silicon is one ULP LOW or exact and never " + "high -- 0 exceptions in 72,907,916 measured rows; interp=" + << std::hex << res[0] << " jit=" << res[1]; } EXPECT_EQ(fcr[1] & kStickyMask, fcr[0] & kStickyMask); ++checked; @@ -249,13 +250,20 @@ TEST(EeRecFpuDivUnitRounding, DivSMatchesInterpExceptWhereTheTruncationLawFires) EXPECT_GT(tier_gaps, 0) << "anti-vacuity: the operand pool stopped producing " "saturating quotients, so the allowance above is " "dead code that could hide a real divergence"; - EXPECT_GT(law_gaps, 0) << "anti-vacuity: no operand pair reached the truncation " - "law, so this test is asserting engine agreement under " - "a different name"; + EXPECT_GT(gaps, 0) << "anti-vacuity: the two engines agreed on every operand " + "pair, so this test is asserting engine agreement under " + "a different name"; + // The A