mysterious remquo bug

[Greg Troxel <gdt%lexort.com@localhost>]

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