//! The Emotion Engine's FPU, as the tfrag lighting pass uses it. //! //! Spec: `docs/formats/ps2-float.md`. //! //! **This is not `f32`.** The EE has no NaN, no infinity and no denormals: an //! exponent of 0 *is* zero and an exponent of 255 is an ordinary value. Every //! result truncates toward zero rather than rounding to nearest, and the adder //! discards the low bits of the smaller operand that sit more than one bit //! below the larger's LSB. //! //! That last rule is the one that matters for fidelity. `(a + b)` and //! `a + b` in Rust can differ in the last bit, and the lighting accumulates //! those differences across every vertex until the colours visibly drift. //! //! Arithmetic is therefore done on raw `u32` bit patterns. Passing `f32` in //! would mean `to_bits`, an operation, `from_bits` on every step -- and the //! temptation to "just use f32" that follows from it. /// Sign bit. pub const SIGN: u32 = 0x8000_0000; /// Largest representable magnitude. pub const MAX: u32 = 0x7fff_ffff; /// 1.0. pub const ONE: u32 = 0x3f80_0000; fn exp(x: u32) -> i32 { ((x >> 23) & 0xff) as i32 } fn man(x: u32) -> u64 { ((x & 0x7f_ffff) | 0x80_0000) as u64 } /// Pack `sign * mag * 2^(lsb_exp)`, truncating to 24 significant bits. /// /// Overflow clamps to ±`MAX` and underflow flushes to ±0, because the hardware /// has neither an infinity nor a denormal to land on. fn pack(sign: u32, mag: u64, lsb_exp: i32) -> u32 { if mag == 0 { return sign; } let k = 63 - mag.leading_zeros() as i32; let e = lsb_exp + (k - 23); let m = if k >= 23 { mag >> (k - 23) } else { mag << (23 - k) }; if e > 255 { sign | MAX } else if e < 1 { sign } else { sign | ((e as u32) << 23) | ((m as u32) & 0x7f_ffff) } } /// `MUL`. Truncating. pub fn mul(a: u32, b: u32) -> u32 { let s = (a ^ b) & SIGN; if exp(a) == 0 || exp(b) == 0 { return s; } // man(a) * man(b) has its binary point at bit 46, so >> 23 leaves bit 23 // worth 2^(ea + eb - 254). pack(s, (man(a) * man(b)) >> 23, exp(a) + exp(b) - 127) } /// `ADD`, truncating, with the EE's operand alignment. pub fn add(a: u32, b: u32) -> u32 { let (ea, eb) = (exp(a), exp(b)); if ea == 0 || eb == 0 { if ea != 0 { return a; } if eb != 0 { return b; } // Both zero: 0 unless both are negative zero, which is -0. return if a & b & SIGN != 0 { SIGN } else { 0 }; } let (hi, lo) = if ea >= eb { (a, b) } else { (b, a) }; let d = exp(hi) - exp(lo); // A 25-bit exponent gap leaves nothing of the smaller operand to add. if d >= 25 { return hi; } let lo = if d >= 1 { lo & (u32::MAX << (d - 1)) } else { lo }; let sv = |x: u32, m: u64| if x & SIGN != 0 { -(m as i64) } else { m as i64 }; let sum = sv(hi, man(hi) << d) + sv(lo, man(lo)); if sum == 0 { return 0; } pack(if sum < 0 { SIGN } else { 0 }, sum.unsigned_abs(), exp(lo)) } /// `SUB`. pub fn sub(a: u32, b: u32) -> u32 { add(a, b ^ SIGN) } /// `DIV`. Truncated quotient; division by zero gives ±`MAX`, not an infinity, /// because the EE has none. pub fn div(a: u32, b: u32) -> u32 { let s = (a ^ b) & SIGN; if exp(b) == 0 { return s | MAX; } if exp(a) == 0 { return s; } let q = (man(a) << 24) / man(b); pack(s, q, exp(a) - exp(b) + 127 - 1) } /// `SQRT` of the magnitude, truncated. pub fn sqrt(x: u32) -> u32 { if exp(x) == 0 { return 0; } let e = exp(x) - 127; let odd = e.rem_euclid(2); let r = isqrt(man(x) << (23 + odd)); pack(0, r, (e - odd) / 2 + 127) } /// `RSQRT a, b`: the square root is truncated *first*, then divided. Doing it /// as one step gives a different answer. pub fn rsqrt(a: u32, b: u32) -> u32 { div(a, sqrt(b)) } /// `ITOF12`: signed integer scaled by 1/4096, truncating integers wider than 24 /// bits. pub fn itof12(i: i32) -> u32 { if i == 0 { return 0; } pack( if i < 0 { SIGN } else { 0 }, i.unsigned_abs() as u64, 127 + 23 - 12, ) } /// Ordering key. The sign bit has to be folded in or every negative compares /// above every positive. fn key(x: u32) -> i64 { if x & SIGN != 0 { -((x & MAX) as i64) } else { x as i64 } } /// `MAX`, where -0 < 0. pub fn max(a: u32, b: u32) -> u32 { if key(b) > key(a) { b } else { a } } /// `MIN`, where -0 < 0. pub fn min(a: u32, b: u32) -> u32 { if key(b) < key(a) { b } else { a } } /// Integer square root, exact. fn isqrt(n: u64) -> u64 { let mut r = (n as f64).sqrt() as u64; while r * r > n { r -= 1; } while (r + 1) * (r + 1) <= n { r += 1; } r } #[cfg(test)] mod tests { use super::*; fn f(x: f32) -> u32 { x.to_bits() } fn v(x: u32) -> f32 { f32::from_bits(x) } #[test] fn one_is_one() { assert_eq!(v(mul(ONE, ONE)), 1.0); assert_eq!(v(add(ONE, 0)), 1.0); assert_eq!(v(div(ONE, ONE)), 1.0); } #[test] fn there_is_no_denormal() { // The smallest f32 denormal would be exp 0 with a nonzero mantissa. // On the EE that bit pattern is simply zero. let denormal = 1u32; assert_eq!(exp(denormal), 0); assert_eq!(v(mul(denormal, ONE)), 0.0); } #[test] fn there_is_no_infinity() { // Dividing by zero saturates rather than producing an infinity. The // comparison is on bits: `MAX` is exponent 255 with a non-zero // mantissa, which Rust reads back as NaN precisely because the EE has // no infinity for it to be. That asymmetry is the whole reason this // module exists. assert_eq!(div(ONE, 0), MAX); assert!(v(MAX).is_nan()); assert_eq!(sqrt(f(-1.0)), f(1.0)); } #[test] fn overflow_saturates() { let big = f(3.0e38); assert_eq!(mul(big, big), MAX); } #[test] fn results_truncate_toward_zero() { // 1/3 is 0.333...; the EE truncates rather than rounding to nearest. let third = div(ONE, f(3.0)); assert!(v(third) < 1.0 / 3.0, "got {}", v(third)); assert!(v(third) > 1.0 / 3.0 - 1e-6); } #[test] fn a_large_operand_absorbs_a_small_one() { let big = f(1.0e10); let tiny = f(1.0e-10); // A 40-bit exponent gap is past the adder's 25-bit window. assert_eq!(v(add(big, tiny)), v(big)); } #[test] fn sqrt_and_rsqrt_are_truncated_at_each_step() { let two = f(2.0); let s = sqrt(two); // Truncated, so the root never overshoots. assert!(v(s) * v(s) <= 2.0); // RSQRT is sqrt-then-divide, with truncation between the two. That is // the contract to preserve: computing 1/sqrt(x) in one go agrees only // by luck, and disagrees as soon as either step lands on a different // bit. for x in [2.0f32, 3.0, 7.0, 11.5, 0.25, 1000.0] { assert_eq!(rsqrt(ONE, f(x)), div(ONE, sqrt(f(x))), "x={x}"); } } #[test] fn negatives_order_below_positives() { let neg = f(-5.0); let pos = f(5.0); assert_eq!(max(neg, pos), pos); assert_eq!(min(neg, pos), neg); assert_eq!(max(pos, neg), pos); // -0 and +0 compare equal, and the hardware keeps the first operand, // so both min and max return it. An unsigned compare would get this // backwards for every negative number. assert_eq!(min(SIGN, 0), SIGN); assert_eq!(max(SIGN, 0), SIGN); } #[test] fn itof12_is_a_twelve_bit_fraction() { assert_eq!(v(itof12(4096)), 1.0); assert_eq!(v(itof12(-4096)), -1.0); assert_eq!(v(itof12(2048)), 0.5); assert_eq!(itof12(0), 0); } #[test] fn subtraction_is_addition_of_the_negated_operand() { assert_eq!(sub(f(5.0), f(3.0)), add(f(5.0), f(-3.0))); assert!((v(sub(f(5.0), f(3.0))) - 2.0).abs() < 1e-6); } }