Reimplement php_round_helper() using modf()
#12220
Merged
Conversation
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
This change makes the implementation much easier to understand, by explicitly handling the various cases. It fixes rounding for `0.49999999999999994`, because no loss of precision happens by adding / subtracing `0.5` before turning the result into an integral float. Instead the fractional parts are explicitly compared. see phpGH-12143 (this fixes one of the reported cases) Closes phpGH-12159 which was an alternative attempt to fix the rounding issue for `0.49999999999999994`
TimWolla
commented
Sep 15, 2023
ext/standard/math.c
Outdated
Comment on lines
102
to
106
| if (fractional >= 0.5) { | ||
| return integral + copysign(1.0, integral); | ||
| } | ||
|
|
||
| return integral; |
There was a problem hiding this comment.
clang is able to compile this branchless as per: https://godbolt.org/z/5cb8W1416, which is likely a good thing.
|
@TimWolla Since |
This is for visual consistency with the other modes.
|
Are you able to do some perf tests to see if there's any difference? |
This makes the code even clearer to understand and also improves the assembly, allowing the compiler to use an actual jump table for the switch cases.
Using #include <math.h>
#include <stdbool.h>
#include <stdio.h>
#include <stdint.h>
#ifndef PHP_ROUND_HALF_UP
#define PHP_ROUND_HALF_UP 0x01 /* Arithmetic rounding, up == away from zero */
#endif
#ifndef PHP_ROUND_HALF_DOWN
#define PHP_ROUND_HALF_DOWN 0x02 /* Down == towards zero */
#endif
#ifndef PHP_ROUND_HALF_EVEN
#define PHP_ROUND_HALF_EVEN 0x03 /* Banker's rounding */
#endif
#ifndef PHP_ROUND_HALF_ODD
#define PHP_ROUND_HALF_ODD 0x04
#endif
double old(double value, int mode) {
double tmp_value;
if (value >= 0.0) {
tmp_value = floor(value + 0.5);
if ((mode == PHP_ROUND_HALF_DOWN && value == (-0.5 + tmp_value)) ||
(mode == PHP_ROUND_HALF_EVEN && value == (0.5 + 2 * floor(tmp_value/2.0))) ||
(mode == PHP_ROUND_HALF_ODD && value == (0.5 + 2 * floor(tmp_value/2.0) - 1.0)))
{
tmp_value = tmp_value - 1.0;
}
} else {
tmp_value = ceil(value - 0.5);
if ((mode == PHP_ROUND_HALF_DOWN && value == (0.5 + tmp_value)) ||
(mode == PHP_ROUND_HALF_EVEN && value == (-0.5 + 2 * ceil(tmp_value/2.0))) ||
(mode == PHP_ROUND_HALF_ODD && value == (-0.5 + 2 * ceil(tmp_value/2.0) + 1.0)))
{
tmp_value = tmp_value + 1.0;
}
}
return tmp_value;
}
double new_jorg(double value, int mode) {
double tmp_value;
if (value >= 0.0) {
tmp_value = floor(value + 0.5);
if ((mode == PHP_ROUND_HALF_DOWN && value == (-0.5 + tmp_value)) ||
(mode == PHP_ROUND_HALF_EVEN && value == (0.5 + 2 * floor(tmp_value/2.0))) ||
(mode == PHP_ROUND_HALF_ODD && value == (0.5 + 2 * floor(tmp_value/2.0) - 1.0)) ||
value < (-0.5 + tmp_value))
{
tmp_value = tmp_value - 1.0;
}
} else {
tmp_value = ceil(value - 0.5);
if ((mode == PHP_ROUND_HALF_DOWN && value == (0.5 + tmp_value)) ||
(mode == PHP_ROUND_HALF_EVEN && value == (-0.5 + 2 * ceil(tmp_value/2.0))) ||
(mode == PHP_ROUND_HALF_ODD && value == (-0.5 + 2 * ceil(tmp_value/2.0) + 1.0)) ||
value > (0.5 + tmp_value))
{
tmp_value = tmp_value + 1.0;
}
}
return tmp_value;
}
double new_tim(double value, int mode) {
double integral, fractional;
fractional = fabs(modf(value, &integral));
switch (mode) {
case PHP_ROUND_HALF_UP:
if (fractional >= 0.5) {
return integral + copysign(1.0, integral);
}
return integral;
case PHP_ROUND_HALF_DOWN:
if (fractional > 0.5) {
return integral + copysign(1.0, integral);
}
return integral;
case PHP_ROUND_HALF_EVEN:
if (fractional > 0.5) {
return integral + copysign(1.0, integral);
}
if (fractional == 0.5) {
bool even = !fmod(integral, 2.0);
if (!even) {
return integral + copysign(1.0, integral);
}
}
return integral;
case PHP_ROUND_HALF_ODD:
if (fractional > 0.5) {
return integral + copysign(1.0, integral);
}
if (fractional == 0.5) {
bool even = !fmod(integral, 2.0);
if (even) {
return integral + copysign(1.0, integral);
}
}
return integral;
}
__builtin_unreachable();
}
static inline uint64_t rotl(const uint64_t x, int k) {
return (x << k) | (x >> (64 - k));
}
static uint64_t s[4];
uint64_t next(void) {
const uint64_t result = s[0] + s[3];
const uint64_t t = s[1] << 17;
s[2] ^= s[0];
s[3] ^= s[1];
s[1] ^= s[2];
s[0] ^= s[3];
s[2] ^= t;
s[3] = rotl(s[3], 45);
return result;
}
int
main() {
s[0] = 0xbe0abf86eeacfd2d;
s[1] = 0x5212a180ba6c1136;
s[2] = 0xbb1f87b46572ab77;
s[3] = 0xab0fde1b8ab187da;
const double step_size = 1.0 / (1ULL << 53);
double sum = 0;
for (size_t i = 0; i < 100000000; i++) {
sum += FUNC(step_size * (next() >> 11), PHP_ROUND_HALF_UP);
}
printf("%.17g\n", sum);
}with old being the unfixed implementation that returns incorrect results for Within a Debian Sid container: |
This comment was marked as resolved.
This comment was marked as resolved.
|
This was on a Intel(R) Core(TM) i5-2430M. The implementation in this PR is indeed the slowest of all of them, with the majority of the time spent in |
|
Differences become much smaller when changing Then all versions are competitive with gcc and Jorg's version is much slower with clang, because the branch misses skyrocket. |
|
With main being: (i.e. |
|
My conclusion from these tests is that the performance heavily depends on the compiler used, the input distribution, and the rounding mode (with my version benefiting from the I'd argue that the |
|
For reference, LLVM seems to have recently merged a faster version of |
|
@TimWolla 's implementation appears to be numerically stable. Regarding the case of slowness, isn't this within an acceptable range considering accuracy? |
@Girgias Also the clang version is tested for the non |
mvorisek
reviewed
Sep 19, 2023
bukka
approved these changes
Sep 19, 2023
|
Now merged, thank you. |
This was referenced Sep 19, 2023
TimWolla
added a commit
to TimWolla/php-src
that referenced
this pull request
Sep 23, 2023
Since phpGH-12220 the implementation of `php_round_helper()`, which performs rounding to an integral value, is easy to verify for correctness up to the floating point precision. If rounding to 0 places is desired, i.e. the userland `round()` function is called with `$precision = 0`, we bypass all logic for the decimal point adjustment and instead directly call `php_round_helper()`. This change fixes the remaining two cases of phpGH-12143 and likely guarantees correct rounding for all possible inputs and `$precision = 0`.
This was referenced Sep 23, 2023
Closed
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.
This change makes the implementation much easier to understand, by explicitly handling the various cases.
It fixes rounding for
0.49999999999999994, because no loss of precision happens by adding / subtracing0.5before turning the result into an integral float. Instead the fractional parts are explicitly compared.see GH-12143 (this fixes one of the reported cases)
Closes GH-12159 which was an alternative attempt to fix the rounding issue for
0.49999999999999994/cc @jorgsowa