JuliaMath / openlibm

High quality system independent, portable, open source libm implementation
https://openlibm.org
Other
515 stars 141 forks source link

C versions of remquo deliver wrong quotient if x == -y #283

Open statham-arm opened 1 year ago

statham-arm commented 1 year ago

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.)