Floating point on 32-bit x86: x87 precision (XPFPA) and fpmath=sse View as Markdown

Hi internals,

While working on 64-bit integers for 32-bit platforms [1], I found floating point results that differ only on 32-bit x86. I would like to reduce the platform differences where PHP is running as much as possible.

Why i386 differs from amd64, arm64 and arm32

On amd64 (SSE2), arm64 and arm32 (VFP), the compiler computes double arithmetic in IEEE-754 binary64: 53-bit mantissa, 11-bit exponent, one rounding per operation.
On i386, compilers use the x87 FPU by default. Its registers have a 64-bit mantissa and a 15-bit exponent, and values are rounded to double only when they are stored to memory.

To get closer to binary64, zend_init_fpu() sets the x87 precision to 53 bits for every request (the XPFPA macros in Zend/zend_float.h). That leaves two differences from the other three platforms:

1. libm runs with the reduced precision too. glibc's i386 pow(), exp(), expm1(), sinh() etc. compute on the x87 unit and need its extended precision to round their result correctly. Inside PHP they lose it, so they are much less accurate than the same functions outside PHP.
2. Only the mantissa is reduced, not the exponent. A result in the subnormal range (below PHP_FLOAT_MIN) is rounded once in the wide register and again when stored, and can end up one ULP off.

In both cases i386 is less accurate compared to the other platforms.

Part 1: libm accuracy -- XPFPA wrappers

#23578 pow(), ** and fpow() https://github.com/php/php-src/pull/23578
#23579 exp() https://github.com/php/php-src/pull/23579
#23581 expm1() https://github.com/php/php-src/pull/23581
#23582 sinh() https://github.com/php/php-src/pull/23582
Examples on i386 today:
var_dump(10 ** -2); // float(0.010000000000000002) vs float(0.01)
var_dump(10 ** -2 === 0.01); // false on i386, true on amd64/arm64
var_dump(exp(10)); // float(22026.465794806725) vs ...718

Each PR adds a small inline helper to Zend/zend_float.h (zend_pow(), zend_exp(), zend_expm1(), zend_sinh()). It switches the FPU to extended precision for the libm call and rounds the result back to double. Where XPFPA_HAVE_CW is 0 -- amd64, arm64, arm32 -- the helper compiles to the plain libm call, so nothing changes there.

With the PRs, nearly all results I checked on i386 are correctly rounded, compared to roughly one in ten before. It is not a guarantee:
a few inputs stay wrong, and one gets 1 ULP worse. Each PR's test contains inputs that are wrong on i386 today and correct with the change.

Cost: two fnstcw/fldcw pairs per call. Measured on an x86_64 Linux host running the i386 build (-m32 -msse2), JIT off. Each function was called 1,000,000 times in a loop, repeated 15 times, and the median time per call is used:
function cost added per call
exp() ~5-7 ns (~5% of the libm call)
expm1() ~5-7 ns
sinh() ~4 ns
**, fpow() ~3-8 ns

JIT should not need any changes: as far as I can see it has no implementation of its own for these functions or for **. It calls the regular function handlers, and ** goes through pow_function() and from there zend_pow(). Compile-time evaluation calls the same handlers. Is this assumption correct, or is there a path I missed?

Part 2: subnormal double rounding -- SSE2 instead of x87

The wrappers can't fix the second problem, as it happens in PHP's own arithmetic:
$a = 5.678650331314532e-158;
$b = 2.3715405351663569e-153;
var_dump($a * $b);
// i386 (x87): 1.34671494457485E-310
// amd64, arm64, arm32: 1.3467149445748E-310 (correctly rounded)

Normal-range results are identical; only some results below PHP_FLOAT_MIN differ.
Compiling with -msse2 -mfpmath=sse makes i386 compute in binary64 like the other platforms and removes the difference.

Unlike the other SIMD code in php-src, which falls back to scalar code or picks an implementation at run time, the floating point unit is fixed at compile time and crashes with
SIGILL on a CPU without SSE2.

This doesn't make Part 1 unnecessary: zend_init_fpu() still reduces the x87 precision, and the system's i386 libm still runs on x87.

The PRs are open for some time already without any comments - comments welcome!
Can we enable -msse2 -mfpmath=sse by default on x86 (require SSE2).
WDYT?

Thanks,
Marc Bennewitz