mysterious remquo bug

Greg Troxel <[email protected]>
Newsgroups gmane.os.netbsd.devel.userlevel
Message-ID <[email protected]>
C99/POSIX define remquo:

  https://pubs.opengroup.org/onlinepubs/9799919799/functions/remquo.html

The basic idea is that return value r is a remainder, and in q is stored
the low-order k bits of the quotient, but with q having the same sign as
the quotient.

The point is to reduce possibly large angles before using sin/cos.

I'm doing all this on netbsd-10, amd64, system gcc 10, intel 9th
generation i7.  The following test program calls remquo with (-90, 90).
This is from a test in geodesiclib in proj of a high-accuracy sincos
implementation, which reduces angles to [45,-45] and then adjusts sign
by quadrant.  But, it can be understood from the remquo specification.

The remainder is 0, and that's ok (-0, but not worrying about that).

The true quotient value is -1, and the low order bits should be
not only congruent to -1, but < 0.

The test program calls remquo twice, one with the -90 in a double, and
once with it as a literal.  The first time gets a positive q -- which
messes up the sin calculation, wwhich is how I got here.  The second is
correct.

Adding -ffloat-store didn't change the output.

I got the same behavior on netbsd-9 earmv7hf-el.

src/lib/libm/src/s_remquo.c is somewhat hard to understand....

Does anyone consider themselves the guardian of numerical correctness?
You are hereby nerdsniped!

----------------------------------------
#include <math.h>
#include <stdio.h>

int
main()
{
	double x = -90.0;
	double r; int q = 0;

	r = remquo(x, 90, &q);
	printf("remquo(x == -90.0, 90) r %f q %d", r, q);
	if (q > 0)
		printf("\tERROR: q positive");
	printf("\n");

	r = remquo(-90.0, 90, &q);
	printf("remquo(-90.0, 90) r %f q %d", r, q);
	if (q > 0)
		printf("\tERROR: q positive");
	printf("\n");
}
----------------------------------------
remquo(x == -90.0, 90) r -0.000000 q 1  ERROR: q positive
remquo(-90.0, 90) r -0.000000 q -1
----------------------------------------
lmpx.com only provides a reader for public news (NNTP) servers. It is not affiliated with the servers or forums shown here and is not responsible for the content of articles, which is written by their respective authors.