# 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) <char...@karney.com> 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) <char...@karney.com> 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; }
