Fix sinh() losing precision on 32-bit x86 - #23582
Open
marc-mabe wants to merge 2 commits into
Open
Conversation
On i386 the x87 FPU is the platform default and zend_init_fpu() clamps it to double precision (53-bit mantissa) for the whole request. libm's sinh() computes on the x87 unit and relies on the extended range to round correctly to double, so while that clamp is in place sinh() is off by one or more ULP. This matters because it makes sinh() silently architecture-dependent: the same script prints different floats on i386 than on x86-64 or arm64, comparisons against precomputed constants fail, and the error propagates into anything built on sinh(). Every expected value is the correctly rounded double, verified against a high-precision reference rather than against another platform, and written as the shortest decimal that round-trips back to it. All arguments are exact binary fractions. The first group is inputs that are already correct everywhere, so a failure points at the affected inputs rather than at the whole function. This test fails on i386 and passes wherever XPFPA_HAVE_CW is 0, which includes x86-64 and arm64. The following commit fixes i386.
Restore extended FPU precision around the libm sinh() call and truncate the result back to double. Why this matters: zend_init_fpu() clamps the x87 FPU to a 53-bit mantissa for the whole request so that PHP's own double arithmetic behaves like IEEE-754 binary64. glibc's sinh(), however, computes on the x87 unit and relies on the extended range to round correctly to double, so the clamp costs sinh() accuracy. The loss is observable from userland and makes sinh() silently architecture-dependent, which breaks float comparisons, cached or serialized results, and any computation that accumulates the error. Measured against a high-precision reference rounded to nearest double, over 720 sampled inputs: i386 before this change 75/720 correctly rounded (10.4%) i386 after this change 715/720 correctly rounded (99.3%) arm64, unchanged 601/720 correctly rounded (83.5%) Note that arm64 is itself not a correctly rounded baseline for sinh(), so these figures are stated against the high-precision reference rather than against another platform. This is an improvement, not a guarantee: sinh(67) regresses by 1 ULP, because the clamped result there is genuinely the correctly rounded one. The accompanying test therefore asserts only inputs that this change fixes and that are also correct on x86-64 and arm64. zend_sinh() in zend_float.h compiles down to a plain sinh() wherever XPFPA_HAVE_CW is 0, which includes x86-64 and arm64, so no other platform is affected.
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Problem
On 32-bit x86,
sinh()returns results that are off by many ULP compared to other platforms, and that are not correctly rounded:Measured against (ex - e-x) / 2 evaluated to 80 decimal digits and rounded to nearest double, over 720 sampled inputs, the i386 result is correctly rounded in only 75 cases (10.4%).
This makes
sinh()silently architecture-dependent, so comparisons against precomputed constants fail, cached or serialized values differ between machines, and the error propagates into anything built onsinh().Cause
init_executor()callszend_init_fpu(), which clamps the x87 FPU to double precision (a 53-bit mantissa) for the whole request so that PHP's owndoublearithmetic behaves like IEEE-754 binary64 rather than carrying x87's 64-bit excess precision.libm's
sinh(), however, computes on that same x87 unit and relies on the extended range internally in order to round correctly todouble. While the clamp is in place it loses that headroom, sosinh()inside PHP is markedly less accurate than the very same libmsinh()is outside PHP.Fix
zend_sinh()inZend/zend_float.hrestores extended precision for the duration of the libm call and truncates the result back todoubleon return, using the XPFPA macros already in that header.PHP_FUNCTION(sinh)calls it.Wherever
XPFPA_HAVE_CWis 0 — x86-64, arm64, and every other target without x87 precision control —zend_sinh()compiles down to a plainsinh()call, so no other platform is affected in behaviour or code generation.Result on i386, same 720 inputs and same reference:
The accuracy figures are stated against a high-precision reference rather than against another platform, because for
sinh()no mainstream platform is a correct baseline: glibc on arm64 is itself only correctly rounded for 601/720 (83.5%) of the same inputs. Comparing i386 against arm64 would therefore understate the improvement and mislabel many fixed cases as regressions — against arm64 the change looks like 534 improvements and 26 regressions, whereas against the true value it is 641 improvements and a single regression.Known limitations
This is a large improvement, not a guarantee.
sinh(67)regresses by 1 ULP, because the clamped result there is genuinely the correctly rounded one. It is the only regression in the sample, and the test asserts only inputs that this change fixes.Tests
ext/standard/tests/math/sinh_fpu_precision.phpt.All arguments are exact binary fractions and every expected value is the correctly rounded double, written as the shortest decimal that round-trips back to it. Inputs were chosen from the intersection of cases that are wrong on clamped i386, correct on i386 with this change, and correct on x86-64/arm64.
Verified locally with a full
make clean+./configure+makeper configuration, on i386 and arm64:Relationship to the other XPFPA PRs
This is the same class of bug as the
pow()/fpow(),exp()andexpm1()fixes and adds a helper to the same part ofZend/zend_float.h.Notes
Zend/zend_float.h