C fmod function: Floating point error and optimization

Question

I'm trying to calculate the true course from one point to anoter on the surface of the earth in as few CPU cycles as possible. The result should be a double 0 <= tc < 360, however in a few special cases i get the result 360 (should be reported as 0). I realize that this is due to machine precision when working with fmod and floating point number, but what will be the most efficient workaround of the problem?

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

#define EPS 1e-15 // EPS a small number ~ machine precision
#define R2D 57.295779513082320876798154814105   //multiply radian with R2D to get degrees
#define D2R 0.01745329251994329576923690768489  //multiply degrees with D2R to get radians
#define TWO_PI 6.283185307179586476925286766559 //2*Pi

/*----------------------------------------------------------------------------
* Course between points

* We obtain the initial course, tc1, (at point 1) from point 1 to point 2
* by the following. The formula fails if the initial point is a pole. We can
* special case this with as IF statement
----------------------------------------------------------------------------
  Implementation
  Argument 1: INPUT - Pointer to double containing Latitude  of point 1 in degrees
  Argument 2: INPUT - Pointer to double containing Longitude of point 1 in degrees
  Argument 3: INPUT - Pointer to double containing Latitude  of point 2 in degrees
  Argument 4: INPUT - Pointer to double containing Longitude of point 2 in degrees

  RETURN: Double containing initial course in degrees from point1 to point 2
--------------------------------------------------------------------------*/
double _stdcall CourseInitial (double *lat1, double *lon1, double *lat2, double *lon2)
{
    double radLat1 = D2R * *lat1;
    double radLon1 = D2R * *lon1;
    double radLat2 = D2R * *lat2;
    double radLon2 = D2R * *lon2;
    double tc = 0;

    if (cos(radLat1) < EPS) {     // EPS a small number ~ machine precision
      if (radLat1 > 0) {
          tc = 180;  //  starting from N pole
      } else {
          tc = 0;    //  starting from S pole
      }
    } else {
      // Calculate true course [180, 540]
      tc = R2D * atan2(sin(radLon2-radLon1),
                       cos(radLat1) * tan(radLat2) - sin(radLat1) * cos(radLon2-radLon1)
                      ) + 360;
    }

    //Return 0 <= true course < 360
    return fmod(tc, 360);
}

int main(void)
{
    double lat1 = 89;
    double lon1 = 17;
    double lat2 = 68;
    double lon2 = -163;
    double tc = 0;
    tc = CourseInitial(&lat1, &lon1, &lat2, &lon2);
    printf("The course from point 1 to 2 is: %.5f", tc);

    return 0;

}

Output:

The course from point 1 to 2 is: 360.00000

The calculation looks fine, `fmod` shouldn't be having trouble with a representable integer value like 360. You're asking printf to round off the value at 5 decimals though, so anything between `359.999995 <= n < 360.0000` will look a bit iffy in the trace. — doynax, Mar 03 '15 at 17:04
Use `printf("%.*e", DECIMAL_DIG - 1, tc);` to print `tc` with sufficient precision. See http://stackoverflow.com/q/16839658/2410359 — chux - Reinstate Monica, Mar 03 '15 at 17:25
"... this is due to machine precision when working with fmod and floating point number" is not so. Expect `fmod()` to be exact [Is fmod() exact when y is an integer?](http://stackoverflow.com/a/20929050/2410359) — chux - Reinstate Monica, Mar 03 '15 at 18:25
All I can tell you is this code sample prints 0.00000 on my Linux box w gcc 4.1.2 and also prints 0.0000 on my windows box using Visual Studio 2010 (yeah, I'm behind on my version), as you would expect. What compiler are you using? — JohnH, Mar 03 '15 at 19:32
I use MinGW together with Eclipse Luna on Win 7. When I compile without optimization I get 0.00000, however when I compile with -O1 or higher optimization I get 360.00000. — Arild M Johannessen, Mar 03 '15 at 20:14
The problem still exist even though I replace the printf line with printf("The course from point 1 to 2 is: %.17g", tc); — Arild M Johannessen, Mar 03 '15 at 20:17
Why 17? Unless you do a test `if (tc == 360.0)` using 5, 17 or what ever is guessing the precision of the calculation. What is the result of `printf("%.*e", DECIMAL_DIG - 1, tc);`? If all else fails, add `if (tc >= 360.0) tc = 0.0;` — chux - Reinstate Monica, Mar 03 '15 at 20:43
I'm trying to use DECIMAL_DIG but I get an error. 'DECIMAL_DIG' undeclared (first use in this function). Haven't figured out how to use it yet. — Arild M Johannessen, Mar 03 '15 at 21:21
Add `#include `. If that does not work, use `DBL_DIG + 3`. If that does not work, use `printf("%.*e", 40, tc);`. The point is, when code has trouble, test if `tc` is <, == or > 360.0. — chux - Reinstate Monica, Mar 03 '15 at 21:55
The problem of formatting floats is a dead end. See answer below and comments. It seems that the problem is not in the formatting but in the return value of fmod. — Arild M Johannessen, Mar 03 '15 at 22:17

Arild M Johannessen · Accepted Answer · 2015-03-10T09:06:54.327

Location of problem

The difference in the result between the two different levels of optimization occurs when calculating radLon2-radLon1. The result of this calculation is shown here.

-O0 radLon2-radLon1 = 0xC00921FB54442D19
-O1 radLon2-radLon1 = 0xC00921FB54442D18

The difference is in the least significant bit bringing the -O0 result past the real value of pi.

-O0 radLon2-radLon1 = -3.14159265358979356008717331861
Pi with 50 decimals =  3.14159265358979323846264338327950288419716939937510
-O1 radLon2-radLon1 = -3,14159265358979311599796346854

-O0 calculation loads one value on the Floating Point Stack and subtracts the other with the FSUB command (line 004014a8:). The -O1 calculation loads both values onto the Floating Point Stack and subtracts them with the FSUBP command (line 00401480:)

Disassembly after compilation with no optimization -O0

(...)
004014a5:   fld QWORD PTR [ebp-0x30]     // ST(0) = radLon2         (from [ebp-0x30])
004014a8:   fsub QWORD PTR [ebp-0x20]    // ST(0) = ST(0) - radLon1 (from [ebp-0x20])
004014ab:   fstp QWORD PTR [esp]         // [esp] = (radLon2-radLon1) = 0xC00921FB54442D19
004014ae:   call 0x4080d0 <sin>          // ST(0) = sin(radLon2-radLon1)
(...)
-------------------------------------------------------------------------------------
Significant values:

radLon2 =              0xC006C253F2B53437  (-2.84488668075075734620327239099     )     
radLon1 =              0x3FD2FD3B0C77C70D  ( 0.296705972839036047350447233839    )    
radLon2-radLon1 =      0xC00921FB54442D19  (-3.14159265358979356008717331861     ) 
sin(radLon2-radLon1) = 0x3CB72D0000000000  ( 3.21628574467824890348310873378e-16 )                       

Later atan2(y, x) is calculated with these values
x =           0x3FF0B04ED1755590  ( 1.04304391688978981278523860965     )
y =           0x3CB72D0000000000  ( 3.21628574467824890348310873378e-16 )
atan2(y, x) = 0x3CB63828CAA39659  ( 3.08355735803412799607014393888e-16 )

Disassembly after compilation with optimization -O1

(...)
29            double radLon2 = D2R * *lon2;
0040146c:   fld QWORD PTR ds:0x40a0c0                   // ST(0) = D2R          (from [ds:0x40a0c0]) 
00401472:   fmul QWORD PTR [esp+0x30]                   // ST(0) = ST(0) * lon2 (from [ESP+0x30])
                                                        // ST(0) = -2.8448866807507573
27            double radLon1 = D2R * *lon1;
00401476:   fld QWORD PTR ds:0x40a0c0                   // ST(0) = D2R          (from [ds:0x40a0c0]) 
0040147c:   fmul QWORD PTR [esp+0x20]                   // ST(0) = ST(0) * lon1 (from [ESP+0x20])
                                                        // ST(0) = 0.29670597283903605 (radLon1)
                                                        // ST(1) = -2.8448866807507573 (radLon2)
00401480:   fsubp st(1),st                              // ST(1) = ST(1)-ST(0) then POP stack 
00401482:   fst QWORD PTR [esp+0x20]                    // [esp+0x20] = (radLon2-radLon1) = 0xC00921FB54442D18 
(...)

40              tc = R2D * atan2(sin(radLon2-radLon1),
00401492:   fld QWORD PTR [esp+0x20]                    // ST(0)=(radLon2-radLon1)
00401496:   fstp QWORD PTR [esp]                        // [esp]=(radLon2-radLon1)
00401499:   call 0x4080e0 <sin>                         // ST(0)=sin(radLon2-radLon1)
(...)
-------------------------------------------------------------------------------------
Significant values

radLon2 =              0xC006C253F2B53437  (-2.84488668075075734620327239099     )     
radLon1 =              0x3FD2FD3B0C77C70D  ( 0.296705972839036047350447233839    )    
radLon2-radLon1 =      0xC00921FB54442D18  (-3,14159265358979311599796346854     ) 
sin(radLon2-radLon1) = 0xBCA1A60000000000  (-1.22460635382237725821141793858e-16 )                       

Later atan2(y, x) is calculated with these values
x =           0x3FF0B04ED1755590  ( 1.04304391688978981278523860965     )
y =           0xBCA1A60000000000  (-1.22460635382237725821141793858e-16 )
atan2(y, x) = 0xBCA0EB8D90F27437  (-1.1740697913027295863036855646E-16 ) 
=====================================================================================

Attempted solution

In CourseInitial() radLon1 and radLon2 is not used independently. So I tried the following.

double radDeltaLon = D2R * (*lon2-*lon1);
(...)
tc = R2D * atan2(sin(radDeltaLon),
                 cos(radLat1) * tan(radLat2) - sin(radLat1) * cos(radDeltaLon)
                ) + 360;

This did not work. Debugging showed that the problematic value close to Pi showed up another place in the code and the end result was the same.

One solution

At the end of each of the defined constants I added an L and with this converting them to Long Doubles (80-bit floating point numbers). This is the same precision that the CPU has in it's Floating Point Registers and solved the problem in some cases.

#define R2D 57.295779513082320876798154814105L   //multiply radian with R2D to get degrees
#define D2R 0.01745329251994329576923690768489L  //multiply degrees with D2R to get radians
#define TWO_PI 6.283185307179586476925286766559L //2*Pi

Final solution

// Calculate true course [-180, 180)
tc = atan2(sin(radDeltaLon),
           cos(radLat1) * tan(radLat2) - sin(radLat1) * cos(radDeltaLon)
          );

if (fabs(tc) < EPS) {
    tc = 0;  //Prevents fmod(tc, 360) from returning 360 due to rounding error
} else {
    tc *= R2D; //Convert to degrees after tc has been checked for machine precision
    tc += 360; //tc [180, 540)
}
return fmod(tc, 360); // returns tc [0, 360)

score 0 · Answer 2 · answered Mar 07 '15 at 09:34

A comparison of the given value of the constant D2R and the closest 64 bit and 80 bit floating point number.

80 bit 0x3FF98EFA351294E9C8AF = 1.745329251994329577083321213687439055206596094649285078048706054e-2
D2R                           = 1.745329251994329576923690768489e-2
80 bit 0x3FF98EFA351294E9C8AE = 1.74532925199432957691391462423657898739293159451335668563842773e-2

64 bit 0x3F91DF46A2529D3A     = 1.745329251994329894381863255148346070200204849243164e-2
D2R                           = 1.745329251994329576923690768489e-2
64 bit 0x3F91DF46A2529D39     = 1.7453292519943295474371680597869271878153085708618164e-2

These are the values chosen by my compiler:

0x3FF98EFA351294E9C8AE = 0 011111111111001 1000111011111010001101010001001010010100111010011100100010101110
0x3F91DF46A2529D39     = 0     01111111001  0001110111110100011010100010010100101001110100111001

The conversion was performed with tools and information on the web pages listed below:
http://www.exploringbinary.com/binary-converter
http://apfloat.appspot.com
http://en.wikipedia.org/wiki/Extended_precision
http://en.wikipedia.org/wiki/Double-precision_floating-point_format