tech-userlevel archive
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index][Old Index]
remquo(3) issues
# summary
(Note that unless you already understand remquo or you read the docs I
am pointing to very carefully, this is not going to make a lot of
sense.)
While running tests for PROJ 6.3.0 (wip/proj) I got some failures. PROJ
includes geodesticlib and has replacement code for C99-required math
functions, or can use the system-provided functions. I narrowed the
failure down to using the replacement vs system remquo, and found that
they returned different results. I wrote a test program, but I am not
convinced that my test program provokes the root cause of the
PROJ/geographiclib failure.
The issues are the precise semantics of the returned partial quotient
https://pubs.opengroup.org/onlinepubs/9699919799/functions/remquo.html
and in C++, not in use here, but ~obviously text trying to say the same
thing:
https://en.cppreference.com/w/cpp/numeric/math/remquo
So for the language lawyers, assuming n=3 and thus modulo 8, it seems
that if the quotient is 0, then the answers 0, 8, 16, and so on are all
entirely acceptable.
For arguments x, y of -45.0, 90, the quotient is -0.5, and hence the
sign is negative.
We do not appear to have any atf tests exercising remquo.
It seems obviously in order to look at our remquo implementation and
compare to FreeBSD's.
So I wonder if anyone understands this. Once again, I hope I am not the
expert!
# results from PROJ regression testing
PROJ regression tests (geodtest, specifically) failed on netbsd-8 amd64,
and passed on macos 10.13. They also pass on Linux and FreeBSD, say the
PROJ people and their CI infrastructure.
I get some lines (see my printf below) when input is nan, and n (in this
code, the modulo-2^3 quotient) differs. However, I suspect the partial
quotient is not valid for a nan input.
remquo mismatch delta nan z nan z_c99 nan n 0 n_c99 2146959360 IN nan 90.000000
remquo mismatch delta nan z nan z_c99 nan n 0 n_c99 1 IN nan 90.000000
remquo mismatch delta nan z nan z_c99 nan n 0 n_c99 0 IN nan 90.000000
remquo mismatch delta nan z nan z_c99 nan n 0 n_c99 0 IN nan 90.000000
remquo mismatch delta nan z nan z_c99 nan n 0 n_c99 0 IN nan 90.000000
remquo mismatch delta nan z nan z_c99 nan n 0 n_c99 30716 IN nan 90.000000
remquo mismatch delta nan z nan z_c99 nan n 0 n_c99 30716 IN nan 90.000000
I also get a lot of lines that look like:
remquo mismatch delta 0.000000 z 39.214461 z_c99 39.214461 n 8 n_c99 0 IN 39.214461 90.000000
which appear to be an odd choice of n, but not wrong according to the
spec as I read it.
Here is my modified remquox (with original copyright), which basically is the code in PROJ and
then a second call to the system remquo and comparison.
/*
* This is a C implementation of the geodesic algorithms described in
*
* C. F. F. Karney,
* Algorithms for geodesics,
* J. Geodesy <b>87</b>, 43--55 (2013);
* https://doi.org/10.1007/s00190-012-0578-z
* Addenda: https://geographiclib.sourceforge.io/geod-addenda.html
*
* See the comments in geodesic.h for documentation.
*
* Copyright (c) Charles Karney (2012-2019) <charles%karney.com@localhost> and licensed
* under the MIT/X11 License. For more information, see
* https://geographiclib.sourceforge.io/
*/
static real remquox(real x, real y, int* n) {
real z = remainderx(x, y);
if (n) {
real
a = remainderx(x, 2 * y),
b = remainderx(x, 4 * y),
c = remainderx(x, 8 * y);
*n = (a > z ? 1 : (a < z ? -1 : 0));
*n += (b > a ? 2 : (b < a ? -2 : 0));
*n += (c > b ? 4 : (c < b ? -4 : 0));
if (y < 0) *n *= -1;
if (y != 0) {
if (x/y > 0 && *n <= 0)
*n += 8;
else if (x/y < 0 && *n >= 0)
*n -= 8;
}
}
int n_c99_val;
int *n_c99 = &n_c99_val;
real z_c99 = remquo(x, y, n_c99);
if (z != z_c99 || *n != *n_c99) {
if (z == NAN && z_c99 == NAN) {
printf("remquot both nan\n");
} else {
printf("remquo mismatch delta %f z %f z_c99 %f n %d n_c99 %d IN %f %f\n", z - z_c99, z, z_c99, *n, *n_c99, x, y);
}
} else {
;
/* printf("remquo OK\n"); */
}
return z;
}
# test program
PROJ is far from a minimal reproduction, so I wrote a test program that
compares the the code from PROJ to the system. I ran it on netbsd-8
amd64, i386 and earmv7hf-el with the same results, and also on macos
10.13.
I get output:
q 1.500000 qi 2 z 45.000000 z_sys 45.000000 z_calc 45.000000 n_x 2 n_sys 2 IN -135.000000 -90.000000 ok
q 0.500000 qi 0 z -45.000000 z_sys -45.000000 z_calc -45.000000 n_x 8 n_sys 0 IN -45.000000 -90.000000 BAD_N
q -0.500000 qi 0 z 45.000000 z_sys 45.000000 z_calc 45.000000 n_x -8 n_sys 0 IN 45.000000 -90.000000 BAD_N SIGN:z_sys
q -1.500000 qi -2 z -45.000000 z_sys -45.000000 z_calc -45.000000 n_x -2 n_sys -2 IN 135.000000 -90.000000 ok
q 0.500000 qi 0 z 45.000000 z_sys 45.000000 z_calc 45.000000 n_x 8 n_sys 0 IN 45.000000 90.000000 BAD_N
q 1.500000 qi 2 z -45.000000 z_sys -45.000000 z_calc -45.000000 n_x 2 n_sys 2 IN 135.000000 90.000000 ok
FAILURES 3
There are several cases of the replacement having 8 and the system
having 0. As I see it (integer) 0 has a positive sign, but perhaps that
is unfair since the math integer 0 is not postive.
More concerning is the value of 0 in the third line, compared to -8 in
the geodesic.c code. The sign of the quotient is negative, and 0
isn't.
test_remquo.c follows
----------------------------------------
#include <assert.h>
#include <math.h>
#include <stdio.h>
/* BEGIN extract from proj-6.3.0/src/geodesic.c */
/*
* This is a C implementation of the geodesic algorithms described in
*
* C. F. F. Karney,
* Algorithms for geodesics,
* J. Geodesy <b>87</b>, 43--55 (2013);
* https://doi.org/10.1007/s00190-012-0578-z
* Addenda: https://geographiclib.sourceforge.io/geod-addenda.html
*
* See the comments in geodesic.h for documentation.
*
* Copyright (c) Charles Karney (2012-2019) <charles%karney.com@localhost> and licensed
* under the MIT/X11 License. For more information, see
* https://geographiclib.sourceforge.io/
*/
static double remquox(double x, double y, int* n) {
double z = remainder(x, y);
if (n) {
double
a = remainder(x, 2 * y),
b = remainder(x, 4 * y),
c = remainder(x, 8 * y);
*n = (a > z ? 1 : (a < z ? -1 : 0));
*n += (b > a ? 2 : (b < a ? -2 : 0));
*n += (c > b ? 4 : (c < b ? -4 : 0));
if (y < 0) *n *= -1;
if (y != 0) {
if (x/y > 0 && *n <= 0)
*n += 8;
else if (x/y < 0 && *n >= 0)
*n -= 8;
}
}
return z;
}
/* END extract from proj-6.3.0/src/geodesic.c */
int
test_remquo(double x, double y) {
char *bad;
int ret;
double z_x, z_sys, z_calc;
int n_x_val, n_sys_val;
int *n_x = &n_x_val;
int *n_sys = &n_sys_val;
/* Find the integral quotient per remainder(3). */
double q = x / y;
int qi = rint(q);
n_x_val = -99;
z_x = remquox(x, y, n_x);
assert(n_x == &n_x_val);
n_sys_val = -99;
z_sys = remquo(x, y, n_sys);
assert(n_sys == &n_sys_val);
ret = 0;
bad = "ok";
if (z_x != z_sys) {
ret = 1;
bad = "BAD_Z";
}
if (*n_x != *n_sys) {
ret = 1;
bad = "BAD_N";
}
/* Calculate remainder from computed quotient. */
z_calc = x - qi * y;
printf("q %f qi %d \tz %f z_sys %f z_calc %f\tn_x %d n_sys %d IN %f %f \t%s",
q, qi,
z_x, z_sys, z_calc,
*n_x, *n_sys, x, y,
bad);
if (q < 0 && *n_x >= 0) {
printf(" SIGN:z_x");
}
if (q >= 0 && *n_x < 0) {
printf(" SIGN:z_x");
}
if (q < 0 && *n_sys >= 0) {
ret = 1;
printf(" SIGN:z_sys");
}
if (q >= 0 && *n_sys < 0) {
ret = 1;
printf(" SIGN:z_sys");
}
printf("\n");
return ret;
}
int
main(int argc, char **argv) {
int fail = 0;
fail += test_remquo(-135.0, -90.0);
fail += test_remquo(-45.0, -90.0);
fail += test_remquo(45.0, -90.0);
fail += test_remquo(135.0, -90.0);
fail += test_remquo(45.0, 90.0);
fail += test_remquo(135.0, 90.0);
printf("FAILURES %d\n", fail);
return 0;
}
Home |
Main Index |
Thread Index |
Old Index