openlibm
openlibm copied to clipboard
C versions of remquo deliver wrong quotient if x == -y
The following test program uses all three of remquof, remquo and remquol to compute the remainder of −3 mod +3. The remainder should be zero, and the quotient should be (some low bits of) −1.
#include <stdio.h>
#include <math.h>
int main(void)
{
int quo;
float remf = remquof(-3.0F, 3.0F, &quo);
printf("remquof(-3,+3) = %g, quo=%08x\n", remf, quo);
double rem = remquo(-3.0, 3.0, &quo);
printf("remquo (-3,+3) = %g, quo=%08x\n", rem, quo);
long double reml = remquol(-3.0L, 3.0L, &quo);
printf("remquol(-3,+3) = %Lg, quo=%08x\n", reml, quo);
}
Compiling openlibm for amd64, so that these functions are provided by amd64/s_remquo.S and friends, this behaves as expected:
$ cc test.c libopenlibm.a && ./a.out
remquof(-3,+3) = -0, quo=ffffffff
remquo (-3,+3) = -0, quo=ffffffff
remquol(-3,+3) = -0, quo=ffffffff
But on aarch64, where the C versions in src/s_remquo.c and friends are used, the wrong quotient is delivered:
$ cc test.c libopenlibm.a && ./a.out
remquof(-3,+3) = -0, quo=00000001
remquo (-3,+3) = -0, quo=00000001
remquol(-3,+3) = -0, quo=00000001
Looking at the source code, it appears that there's a special case for |x| = |y|, which always sets *quo = 1 without checking signs.
(Tested against the current HEAD, commit 12f5ffcc990e16f4120d4bf607185243f5affcb8, in each case.)