Fixes issue #3, test = GeodSolve99

cffk · cffk · commit 9416cf72c9fc · 2024-07-31T12:02:23.000-04:00
Problem was that output of sincosd(+/-45) was inconsistent because of
directed rounding by Math.round.  Fix by porting special treatment of
angles +/-45 and +/-30 from C++ library.
diff --git a/geodesic/Math.js.in b/geodesic/Math.js.in
@@ -296,14 +296,21 @@ geodesic.Accumulator = {};
   m.sincosd = function(x) {
     // In order to minimize round-off errors, this function exactly reduces
     // the argument to the range [-45, 45] before converting it to radians.
-    var r, q, s, c, sinx, cosx;
-    r = x % 360;
-    q = Math.round(r / 90);     // If r is NaN this returns NaN
-    r -= 90 * q;
-    // now abs(r) <= 45
-    r *= this.degree;
+    var d, r, q, s, c, sinx, cosx;
+    d = x % 360;
+    q = Math.round(d / 90);     // If d is NaN this returns NaN
+    d -= 90 * q;
+    // now abs(d) <= 45
+    r = d * this.degree;
     // Possibly could call the gnu extension sincos
     s = Math.sin(r); c = Math.cos(r);
+    if (Math.abs(d) === 45) {
+      c = Math.sqrt(0.5);
+      s = m.copysign(c, r);
+    } else if (Math.abs(d) === 30) {
+      c = Math.sqrt(0.75);
+      s = m.copysign(0.5, r)
+    }
     switch (q & 3) {
     case 0:  sinx =  s; cosx =  c; break;
     case 1:  sinx =  c; cosx = -s; break;
@@ -325,22 +332,29 @@ geodesic.Accumulator = {};
   m.sincosde = function(x, t) {
     // In order to minimize round-off errors, this function exactly reduces
     // the argument to the range [-45, 45] before converting it to radians.
-    var r, q, s, c, sinx, cosx;
-    r = x % 360;
-    q = Math.round(r / 90); // If r is NaN this returns NaN
-    r = m.AngRound((r - 90 * q) + t);
-    // now abs(r) <= 45
-    r *= this.degree;
+    var d, r, q, s, c, sinx, cosx;
+    d = x % 360;
+    q = Math.round(d / 90);     // If d is NaN this returns NaN
+    d = m.AngRound((d - 90 * q) + t);
+    // now abs(d) <= 45
+    r = d * this.degree;
     // Possibly could call the gnu extension sincos
     s = Math.sin(r); c = Math.cos(r);
+    if (Math.abs(d) === 45) {
+      c = Math.sqrt(0.5);
+      s = m.copysign(c, r);
+    } else if (Math.abs(d) === 30) {
+      c = Math.sqrt(0.75);
+      s = m.copysign(0.5, r)
+    }
     switch (q & 3) {
     case 0:  sinx =  s; cosx =  c; break;
     case 1:  sinx =  c; cosx = -s; break;
     case 2:  sinx = -s; cosx = -c; break;
     default: sinx = -c; cosx =  s; break; // case 3
     }
     cosx += 0;
-    if (sinx === 0) sinx = m.copysign(sinx, x);
+    if (sinx === 0) sinx = m.copysign(sinx, x+t);
     return {s: sinx, c: cosx};
   };
 
diff --git a/geodesic/test/geodesictest.js b/geodesic/test/geodesictest.js
@@ -651,15 +651,14 @@ describe("geodesic", function() {
 
     it("GeodSolve99", function() {
       // Test case https://github.com/geographiclib/geographiclib-js/issues/3
+      // Problem was that output of sincosd(+/-45) was inconsistent because of
+      // directed rounding by Math.round.  Fix by porting special treatment of
+      // angles +/-45 and +/-30 from C++ library.
       var geod = g.WGS84,
           inv = geod.Inverse(45, 0, -45, 179.572719);
-//        inv = geod.Inverse(45, 0, -45, 179);
       assert.approx(inv.azi1, 90.00000028, 1e-8  );
       assert.approx(inv.azi2, 90.00000028, 1e-8  );
       assert.approx(inv.s12, 19987083.007, 0.5e-3);
-//    assert.approx(inv.azi1, 90.20248882, 1e-8  );
-//    assert.approx(inv.azi2, 90.20248882, 1e-8  );
-//    assert.approx(inv.s12, 19941926.020, 0.5e-3);
     });
 
   });
