Add GeographicLib v2.4 parity for trig accuracy, Hypot3, and Jacobi amplitude

GeographicLib/EllipticFunction.Priv.cs

@@ -481,6 +481,79 @@ public double DeltaH(double sn, double cn, double dn)
                 return H(sn, cn, dn) * (PI / 2) / H() - Atan2(sn, cn);
             }
 
+            /// <summary>
+            /// The Jacobi amplitude function.
+            /// </summary>
+            /// <param name="x">the argument.</param>
+            /// <returns>am(<i>x</i>, <i>k</i>).</returns>
+            public double Am(double x)
+            {
+                // This implements DLMF Sec 22.20(ii).
+                // See also Sala (1989), https://doi.org/10.1137/0520100, Sec 5.
+                var tol = Pow(DBL_EPSILON, 0.75);
+                var k2 = _k2;
+                var kp2 = _kp2;
+                if (_k2 == 0)
+                    return x;
+                if (_kp2 == 0)
+                    return Atan(Sinh(x));
+                if (_k2 < 0)
+                {
+                    k2 = -_k2 / _kp2;
+                    kp2 = 1 / _kp2;
+                    x *= Sqrt(_kp2);
+                }
+
+                Span<double> a = stackalloc double[num_], c = stackalloc double[num_];
+                a[0] = 1;
+                c[0] = Sqrt(k2);
+                var b = Sqrt(kp2);
+                var l = 1;
+                for (; l < num_ || GEOGRAPHICLIB_PANIC; ++l)
+                {
+                    a[l] = (a[l - 1] + b) / 2;
+                    c[l] = (a[l - 1] - b) / 2;
+                    b = Sqrt(a[l - 1] * b);
+                    if (!(c[l] > tol * a[l]))
+                        break;
+                }
+
+                var phi = a[l] * x * (1 << l);
+                var phi1 = 0d;
+                for (; l > 0; --l)
+                {
+                    phi1 = phi;
+                    phi = (phi + Asin(c[l] * Sin(phi) / a[l])) / 2;
+                }
+
+                return _k2 < 0 ? phi1 - phi : phi;
+            }
+
+            /// <summary>
+            /// The Jacobi amplitude function and associated elliptic functions.
+            /// </summary>
+            /// <param name="x">the argument.</param>
+            /// <param name="sn">sn(<i>x</i>, <i>k</i>).</param>
+            /// <param name="cn">cn(<i>x</i>, <i>k</i>).</param>
+            /// <param name="dn">dn(<i>x</i>, <i>k</i>).</param>
+            /// <returns>am(<i>x</i>, <i>k</i>).</returns>
+            public double Am(double x, out double sn, out double cn, out double dn)
+            {
+                var phi = Am(x);
+                if (_kp2 == 0)
+                {
+                    sn = Tanh(x);
+                    cn = dn = 1 / Cosh(x);
+                }
+                else
+                {
+                    sn = Sin(phi);
+                    cn = Cos(phi);
+                    dn = Delta(sn, cn);
+                }
+                return phi;
+            }
+
             /// <summary>
             /// The Jacobi elliptic functions.
             /// </summary>
@@ -996,4 +1069,3 @@ public void Reset(double k2, double alpha2, double kp2, double alphap2)
         }
     }
 }
-

GeographicLib/EllipticFunction.cs

@@ -341,6 +341,23 @@ public partial class EllipticFunction
         /// <returns>the periodic function π <i>H</i>(φ, <i>k</i>) / (2 <i>H</i>(<i>k</i>)) -  φ.</returns>
         public double DeltaH(double sn, double cn, double dn) => _priv.DeltaH(sn, cn, dn);
 
+        /// <summary>
+        /// The Jacobi amplitude function.
+        /// </summary>
+        /// <param name="x">the argument.</param>
+        /// <returns>am(<i>x</i>, <i>k</i>).</returns>
+        public double Am(double x) => _priv.Am(x);
+
+        /// <summary>
+        /// The Jacobi amplitude function and associated elliptic functions.
+        /// </summary>
+        /// <param name="x">the argument.</param>
+        /// <param name="sn">sn(<i>x</i>, <i>k</i>).</param>
+        /// <param name="cn">cn(<i>x</i>, <i>k</i>).</param>
+        /// <param name="dn">dn(<i>x</i>, <i>k</i>).</param>
+        /// <returns>am(<i>x</i>, <i>k</i>).</returns>
+        public double Am(double x, out double sn, out double cn, out double dn) => _priv.Am(x, out sn, out cn, out dn);
+
         /// <summary>
         /// The Jacobi elliptic functions.
         /// </summary>

GeographicLib/MathEx.cs

@@ -365,6 +365,15 @@ public static void Norm(ref double x, ref double y)
         /// <returns>Hypotenuse of a right-angled triangle computed as √(1^2+x^2).</returns>
         public static double Hypot(double x) => Hypot(1, x);
 
+        /// <summary>
+        /// Computes the square root of the sum of the squares of x, y and z.
+        /// </summary>
+        /// <param name="x"></param>
+        /// <param name="y"></param>
+        /// <param name="z"></param>
+        /// <returns>Hypotenuse of a right-angled triangle computed as √(x^2+y^2+z^2).</returns>
+        public static double Hypot3(double x, double y, double z) => Hypot(Hypot(x, y), z);
+
         /// <summary>
         /// Square a number.
         /// </summary>
@@ -373,6 +382,26 @@ public static void Norm(ref double x, ref double y)
         [MethodImpl(MethodImplOptions.AggressiveInlining)]
         public static double Sq(double x) => x * x;
 
+        [MethodImpl(MethodImplOptions.AggressiveInlining)]
+        private static bool TryReducedSinCosSpecial(double d, double r, out double s, out double c)
+        {
+            var ad = Abs(d);
+            if (2 * ad == QD)
+            {
+                c = Sqrt(0.5);
+                s = CopySign(c, r);
+                return true;
+            }
+            if (3 * ad == QD)
+            {
+                c = Sqrt(3) / 2;
+                s = CopySign(0.5, r);
+                return true;
+            }
+            s = c = 0;
+            return false;
+        }
+
         /// <summary>
         /// Evaluate the sine and cosine function with the argument in degrees.
         /// </summary>
@@ -391,10 +420,15 @@ public static void SinCosd(double x, out double sinx, out double cosx)
         {
             // 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 = Remquo(x, QD, out var q);
-            r *= Degree;
+            var d = Remquo(x, QD, out var q);
+            var r = d * Degree;
 
             var (s, c) = SinCos(r);
+            if (TryReducedSinCosSpecial(d, r, out var ss, out var cc))
+            {
+                s = ss;
+                c = cc;
+            }
 
             switch ((uint)q & 3U)
             {
@@ -436,11 +470,16 @@ public static void SinCosde(double x, double t, out double sinx, out double cosx
             // the argument to the range [-45, 45] before converting it to radians.
             // This implementation allows x outside [-180, 180], but implementations in
             // other languages may not.
-            var r = AngRound(Remquo(x, QD, out var q) + t); // now abs(r) <= 45
-            r *= Degree;
+            var d = AngRound(Remquo(x, QD, out var q) + t); // now abs(r) <= 45
+            var r = d * Degree;
 
             // g++ -O turns these two function calls into a call to sincos
             var (s, c) = SinCos(r);
+            if (TryReducedSinCosSpecial(d, r, out var ss, out var cc))
+            {
+                s = ss;
+                c = cc;
+            }
 
             switch ((uint)q & 3U)
             {
@@ -467,13 +506,15 @@ public static void SinCosde(double x, double t, out double sinx, out double cosx
         public static double Sind(double x)
         {
             // See sincosd
-            double r;
-            r = Remquo(x, QD, out var q); // now abs(r) <= 45
-            r *= Degree;
+            var d = Remquo(x, QD, out var q); // now abs(r) <= 45
+            var r = d * Degree;
 
             var p = (uint)q;
 
-            r = (p & 1U) != 0 ? Cos(r) : Sin(r);
+            if (TryReducedSinCosSpecial(d, r, out var s, out var c))
+                r = (p & 1U) != 0 ? c : s;
+            else
+                r = (p & 1U) != 0 ? Cos(r) : Sin(r);
             if ((p & 2U) != 0) r = -r;
 
             if (r == 0) r = CopySign(r, x);
@@ -492,12 +533,14 @@ public static double Sind(double x)
         public static double Cosd(double x)
         {
             // See sincosd
-            double r;
-            r = Remquo(x, QD, out var q); // now abs(r) <= 45
-            r *= Degree;
+            var d = Remquo(x, QD, out var q); // now abs(r) <= 45
+            var r = d * Degree;
 
             var p = (uint)(q + 1);
-            r = (p & 1U) != 0 ? Cos(r) : Sin(r);
+            if (TryReducedSinCosSpecial(d, r, out var s, out var c))
+                r = (p & 1U) != 0 ? c : s;
+            else
+                r = (p & 1U) != 0 ? Cos(r) : Sin(r);
             if ((p & 2U) != 0) r = -r;
             return 0d + r;
         }

GeographicLibTests/EllipticFunctionTests.cs

@@ -0,0 +1,35 @@
+using Microsoft.VisualStudio.TestTools.UnitTesting;
+
+namespace GeographicLib.Tests
+{
+    [TestClass]
+    public class EllipticFunctionTests
+    {
+        [DataTestMethod]
+        [DataRow(0d, 1.25d)]
+        [DataRow(1d, 1.25d)]
+        public void TestAm_SpecialModulusCases(double k2, double x)
+        {
+            var ell = new EllipticFunction(k2);
+            var expected = k2 == 0 ? x : System.Math.Atan(System.Math.Sinh(x));
+
+            Assert.AreEqual(expected, ell.Am(x), 1e-15);
+        }
+
+        [TestMethod]
+        public void TestAm_MatchesSncndnAndInverseF()
+        {
+            var ell = new EllipticFunction(0.5);
+            var x = 1.25;
+
+            var phi = ell.Am(x, out var sn, out var cn, out var dn);
+            ell.Sncndn(x, out var sn1, out var cn1, out var dn1);
+
+            Assert.AreEqual(phi, System.Math.Atan2(sn, cn), 1e-15);
+            Assert.AreEqual(x, ell.F(phi), 1e-14);
+            Assert.AreEqual(sn1, sn, 1e-15);
+            Assert.AreEqual(cn1, cn, 1e-15);
+            Assert.AreEqual(dn1, dn, 1e-15);
+        }
+    }
+}

GeographicLibTests/MathExTests.cs

@@ -94,6 +94,30 @@ public void TestSinCosd_Accuracy()
             Assert.AreEqual(dsin9, -dsin123456789);
         }
 
+        [DataTestMethod]
+        [DataRow(30d, 0.5, 0.8660254037844386)]
+        [DataRow(45d, 0.7071067811865476, 0.7071067811865476)]
+        [DataRow(60d, 0.8660254037844386, 0.5)]
+        [DataRow(-30d, -0.5, 0.8660254037844386)]
+        [DataRow(-45d, -0.7071067811865476, 0.7071067811865476)]
+        [DataRow(-60d, -0.8660254037844386, 0.5)]
+        [DataRow(390d, 0.5, 0.8660254037844386)]
+        [DataRow(405d, 0.7071067811865476, 0.7071067811865476)]
+        [DataRow(420d, 0.8660254037844386, 0.5)]
+        public void TestTrigd_ExactSpecialAngles(double angle, double expectedSin, double expectedCos)
+        {
+            Assert.That.EqualsExactly(expectedSin, MathEx.Sind(angle));
+            Assert.That.EqualsExactly(expectedCos, MathEx.Cosd(angle));
+
+            MathEx.SinCosd(angle, out var sinx, out var cosx);
+            Assert.That.EqualsExactly(expectedSin, sinx);
+            Assert.That.EqualsExactly(expectedCos, cosx);
+
+            MathEx.SinCosde(angle, 0, out sinx, out cosx);
+            Assert.That.EqualsExactly(expectedSin, sinx);
+            Assert.That.EqualsExactly(expectedCos, cosx);
+        }
+
         [TestMethod]
         public void TestAtan2d_Accuracy()
         {
@@ -153,6 +177,15 @@ public void TestHypot_GH75651()
             Assert.IsTrue(x > 0);
         }
 
+        [DataTestMethod]
+        [DataRow(3d, 4d, 12d, 13d)]
+        [DataRow(0d, 0d, 0d, 0d)]
+        [DataRow(0d, 0d, double.PositiveInfinity, double.PositiveInfinity)]
+        public void TestHypot3(double x, double y, double z, double expected)
+        {
+            Assert.AreEqual(expected, MathEx.Hypot3(x, y, z));
+        }
+
         [DataTestMethod]
         [DynamicData("Remquo", typeof(MathExTestData))]
         public void TestRemquo(long x, long y, long d, int q)
@@ -284,4 +317,4 @@ public void TestAsinh(long x, long expected)
         }
 #endif
     }
-}
+}

GeographicLibTests/RhumbSolveTests.cs

@@ -54,7 +54,7 @@ public void HiglyProlateOblateSphere(
         /// </summary>
         [DataTestMethod]
         [DataRow(0, 0, 0, 10001965, 89.99999347, 0.0, 0, false, DisplayName = "RhumbSovle8")]
-        [DataRow(0, 0, 60, 20003930, 89.99999347, 1654.69811, 1123576495779545, true, DisplayName = "RhumbSovle9")]
+        [DataRow(0, 0, 60, 20003930, 89.99999347, 1654.69811, 1123576495931784.2, true, DisplayName = "RhumbSolve9")]
         [DataRow(0, 0, 0, 10001966, 89.99999758, double.NaN, double.NaN, false, DisplayName = "RhumbSovle10")]
         [DataRow(0, 0, 60, 20003932, 89.99999758, double.NaN, double.NaN, true, DisplayName = "RhumbSovle11")]
         public void RhumbLinesWithDestinationsCloseToPole(

README.md

@@ -165,6 +165,13 @@ For changes adopted from GeographicLib, please refer the its change log [here](h
 - **FIX**
   - [Strong name validation fails when loading GeographicLib.dll from .NET Framework application](https://github.com/noelex/GeographicLib.NET/issues/33)
 
+### Version 2.4.0 (unreleased)
+- **NEW**
+  - Add `MathEx.Hypot3`.
+  - Add `EllipticFunction.Am` overloads for the Jacobi amplitude function.
+- **FIX**
+  - Improve `MathEx.Sind`, `MathEx.Cosd`, `MathEx.SinCosd` and `MathEx.SinCosde` accuracy for reduced angles at exact multiples of 30°, 45° and 60°.
+
 ### Version 2.3.1 (released 2024/03/20)
 - **NEW**
   - Improve performance of `Geodesic.Direct`, `Rhumb.Direct`, `GeodesicExact.Direct`, `GeodesicExact.Inverse` and their overloads/variants. These are now heap allocation free.
@@ -227,4 +234,4 @@ For changes adopted from GeographicLib, please refer the its change log [here](h
   - Fix stack overflow bug for `Forward(double lon0, double lat, double lon)` and `Reverse(double lon0, double x, double y)` in `TransverseMercatorExact`.
 
 ### Version 1.51.0 (released 2021/03/14)
-Initial stable release.
+Initial stable release.