proj geodesic-signtest Rounding Error

When PROJ4 version 9.3.1 is installed compiled and built from source code, one of the 62 tests fails: 
ctest --rerun-failed --output-on-failure
Test project proj-9.3.1/build
    Start 2: geodesic-signtest
1/1 Test #2: geodesic-signtest ................***Failed    0.00 sec
Line 140: cos(-810) != 0 (6.12323e-17)
Line 161: cos(810) != 0 (6.12323e-17)
Line 173 : sincos accuracy fail
3 failures

Work Around

Ignore the ctest failure. It is due to rounding error.

remquo bug

The PROJ documentation says there was a problem with rounding error in the GNU C runtime library which was reported in 2014 and fixed in 2015, (See https://sourceware.org/bugzilla/show_bug.cgi?id=17569). However The current GNU version of remquo for 810 degrees (expressed as double) with second argument 90 degrees returns 8 and 90 degrees, rather than 9 and 0. This is due to using floating point arithmetic.

PROJ src/tests/geodsigntest.c has a conditional compilation switch OLD_BUGGY_REMQUO which disables test #2

PROJ sincosdx

The source code for sincosdx src/geodesic.c says it uses remquo to reduce the input angle (in degrees) to the range -45 to +45 before converting to radians. This appears to be wrong and it actually converts it to 0 to 90 before converting to radians.

The impact of the remquo bug when given 810 degrees is to convert it to 90. This is converted to radians and the GNU cos function returns 6.12323e-17 rather than exactly 0.

Possible fixes for sincosdx

  1. The tests should be rewritten to allow for floating point imprecision and the documentation updated.
  2. The comment could be corrected to replace -45 to +45 to 0 to 90.
  3. The input angle could be reduced to 0 to 45 degrees before conversion to radians. Meaning the logic for decided which quadrant the input was in would have to replaced by code for deciding which octant instead and complicating the conversion from sin/cos in first octant back to required output.
  4. remquo could be replaced by (possibly 64 bit) integer division and calculation of remainder. This would entail an ugly conversion from double to fixed point.
W.B.Langdon Back Started 20 Feb 2024.