From ed515625934a6fe75b9c30cfff0918d693bfe5a7 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Mon, 30 Mar 2026 13:35:18 +0000 Subject: [PATCH 1/3] Initial plan From 8b9f1dcf475117d5f801a48f8a27bf941c3986f1 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Mon, 30 Mar 2026 13:48:58 +0000 Subject: [PATCH 2/3] feat: add GeographicLib v2.4 parity APIs Agent-Logs-Url: https://github.com/habizdev/GeographicLib.NET/sessions/c8e70dc3-399a-43d9-9694-d0cb700f65bb Co-authored-by: AgileJoshua <564800+AgileJoshua@users.noreply.github.com> --- GeographicLib/EllipticFunction.Priv.cs | 74 ++++++++++++++++++++- GeographicLib/EllipticFunction.cs | 17 +++++ GeographicLib/MathEx.cs | 59 ++++++++++++---- GeographicLibTests/EllipticFunctionTests.cs | 35 ++++++++++ GeographicLibTests/MathExTests.cs | 35 +++++++++- GeographicLibTests/RhumbSolveTests.cs | 2 +- README.md | 9 ++- 7 files changed, 215 insertions(+), 16 deletions(-) create mode 100644 GeographicLibTests/EllipticFunctionTests.cs diff --git a/GeographicLib/EllipticFunction.Priv.cs b/GeographicLib/EllipticFunction.Priv.cs index 6fb78b2..9f69722 100644 --- a/GeographicLib/EllipticFunction.Priv.cs +++ b/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); } + /// + /// The Jacobi amplitude function. + /// + /// the argument. + /// am(x, k). + 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 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; + } + + /// + /// The Jacobi amplitude function and associated elliptic functions. + /// + /// the argument. + /// sn(x, k). + /// cn(x, k). + /// dn(x, k). + /// am(x, k). + 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; + } + /// /// The Jacobi elliptic functions. /// @@ -996,4 +1069,3 @@ public void Reset(double k2, double alpha2, double kp2, double alphap2) } } } - diff --git a/GeographicLib/EllipticFunction.cs b/GeographicLib/EllipticFunction.cs index f9c4bee..79257f0 100644 --- a/GeographicLib/EllipticFunction.cs +++ b/GeographicLib/EllipticFunction.cs @@ -341,6 +341,23 @@ public partial class EllipticFunction /// the periodic function π H(φ, k) / (2 H(k)) - φ. public double DeltaH(double sn, double cn, double dn) => _priv.DeltaH(sn, cn, dn); + /// + /// The Jacobi amplitude function. + /// + /// the argument. + /// am(x, k). + public double Am(double x) => _priv.Am(x); + + /// + /// The Jacobi amplitude function and associated elliptic functions. + /// + /// the argument. + /// sn(x, k). + /// cn(x, k). + /// dn(x, k). + /// am(x, k). + public double Am(double x, out double sn, out double cn, out double dn) => _priv.Am(x, out sn, out cn, out dn); + /// /// The Jacobi elliptic functions. /// diff --git a/GeographicLib/MathEx.cs b/GeographicLib/MathEx.cs index 2699e11..5b6ec8e 100644 --- a/GeographicLib/MathEx.cs +++ b/GeographicLib/MathEx.cs @@ -365,6 +365,15 @@ public static void Norm(ref double x, ref double y) /// Hypotenuse of a right-angled triangle computed as √(1^2+x^2). public static double Hypot(double x) => Hypot(1, x); + /// + /// Computes the square root of the sum of the squares of x, y and z. + /// + /// + /// + /// + /// Hypotenuse of a right-angled triangle computed as √(x^2+y^2+z^2). + public static double Hypot3(double x, double y, double z) => Hypot(Hypot(x, y), z); + /// /// Square a number. /// @@ -391,10 +400,21 @@ 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); + var ad = Abs(d); + if (2 * ad == QD) + { + c = Sqrt(0.5); + s = CopySign(c, r); + } + else if (3 * ad == QD) + { + c = Sqrt(3) / 2; + s = CopySign(0.5, r); + } switch ((uint)q & 3U) { @@ -436,11 +456,22 @@ 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); + var ad = Abs(d); + if (2 * ad == QD) + { + c = Sqrt(0.5); + s = CopySign(c, r); + } + else if (3 * ad == QD) + { + c = Sqrt(3) / 2; + s = CopySign(0.5, r); + } switch ((uint)q & 3U) { @@ -467,13 +498,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 ad = Abs(d); var p = (uint)q; - r = (p & 1U) != 0 ? Cos(r) : Sin(r); + r = (p & 1U) != 0 + ? (2 * ad == QD ? Sqrt(0.5) : (3 * ad == QD ? Sqrt(3) / 2 : Cos(r))) + : CopySign(2 * ad == QD ? Sqrt(0.5) : (3 * ad == QD ? 0.5 : Sin(r)), r); if ((p & 2U) != 0) r = -r; if (r == 0) r = CopySign(r, x); @@ -492,12 +525,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 ad = Abs(d); var p = (uint)(q + 1); - r = (p & 1U) != 0 ? Cos(r) : Sin(r); + r = (p & 1U) != 0 + ? (2 * ad == QD ? Sqrt(0.5) : (3 * ad == QD ? Sqrt(3) / 2 : Cos(r))) + : CopySign(2 * ad == QD ? Sqrt(0.5) : (3 * ad == QD ? 0.5 : Sin(r)), r); if ((p & 2U) != 0) r = -r; return 0d + r; } diff --git a/GeographicLibTests/EllipticFunctionTests.cs b/GeographicLibTests/EllipticFunctionTests.cs new file mode 100644 index 0000000..79e3d54 --- /dev/null +++ b/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); + } + } +} diff --git a/GeographicLibTests/MathExTests.cs b/GeographicLibTests/MathExTests.cs index d6f7316..ccc1b59 100644 --- a/GeographicLibTests/MathExTests.cs +++ b/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 } -} \ No newline at end of file +} diff --git a/GeographicLibTests/RhumbSolveTests.cs b/GeographicLibTests/RhumbSolveTests.cs index 9de0346..27f9854 100644 --- a/GeographicLibTests/RhumbSolveTests.cs +++ b/GeographicLibTests/RhumbSolveTests.cs @@ -54,7 +54,7 @@ public void HiglyProlateOblateSphere( /// [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 = "RhumbSovle9")] [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( diff --git a/README.md b/README.md index f852a91..50aa4c3 100644 --- a/README.md +++ b/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. \ No newline at end of file +Initial stable release. From e3d5293ab8e5fe0c92311ede3e58970709237bd5 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Mon, 30 Mar 2026 13:52:56 +0000 Subject: [PATCH 3/3] chore: finalize v2.4 parity updates Agent-Logs-Url: https://github.com/habizdev/GeographicLib.NET/sessions/c8e70dc3-399a-43d9-9694-d0cb700f65bb Co-authored-by: AgileJoshua <564800+AgileJoshua@users.noreply.github.com> --- GeographicLib/MathEx.cs | 60 +++++++++++++++------------ GeographicLibTests/RhumbSolveTests.cs | 2 +- 2 files changed, 35 insertions(+), 27 deletions(-) diff --git a/GeographicLib/MathEx.cs b/GeographicLib/MathEx.cs index 5b6ec8e..e1bdf94 100644 --- a/GeographicLib/MathEx.cs +++ b/GeographicLib/MathEx.cs @@ -382,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; + } + /// /// Evaluate the sine and cosine function with the argument in degrees. /// @@ -404,16 +424,10 @@ public static void SinCosd(double x, out double sinx, out double cosx) var r = d * Degree; var (s, c) = SinCos(r); - var ad = Abs(d); - if (2 * ad == QD) - { - c = Sqrt(0.5); - s = CopySign(c, r); - } - else if (3 * ad == QD) + if (TryReducedSinCosSpecial(d, r, out var ss, out var cc)) { - c = Sqrt(3) / 2; - s = CopySign(0.5, r); + s = ss; + c = cc; } switch ((uint)q & 3U) @@ -461,16 +475,10 @@ public static void SinCosde(double x, double t, out double sinx, out double cosx // g++ -O turns these two function calls into a call to sincos var (s, c) = SinCos(r); - var ad = Abs(d); - if (2 * ad == QD) + if (TryReducedSinCosSpecial(d, r, out var ss, out var cc)) { - c = Sqrt(0.5); - s = CopySign(c, r); - } - else if (3 * ad == QD) - { - c = Sqrt(3) / 2; - s = CopySign(0.5, r); + s = ss; + c = cc; } switch ((uint)q & 3U) @@ -500,13 +508,13 @@ public static double Sind(double x) // See sincosd var d = Remquo(x, QD, out var q); // now abs(r) <= 45 var r = d * Degree; - var ad = Abs(d); var p = (uint)q; - r = (p & 1U) != 0 - ? (2 * ad == QD ? Sqrt(0.5) : (3 * ad == QD ? Sqrt(3) / 2 : Cos(r))) - : CopySign(2 * ad == QD ? Sqrt(0.5) : (3 * ad == QD ? 0.5 : Sin(r)), 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); @@ -527,12 +535,12 @@ public static double Cosd(double x) // See sincosd var d = Remquo(x, QD, out var q); // now abs(r) <= 45 var r = d * Degree; - var ad = Abs(d); var p = (uint)(q + 1); - r = (p & 1U) != 0 - ? (2 * ad == QD ? Sqrt(0.5) : (3 * ad == QD ? Sqrt(3) / 2 : Cos(r))) - : CopySign(2 * ad == QD ? Sqrt(0.5) : (3 * ad == QD ? 0.5 : Sin(r)), 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; } diff --git a/GeographicLibTests/RhumbSolveTests.cs b/GeographicLibTests/RhumbSolveTests.cs index 27f9854..7933f39 100644 --- a/GeographicLibTests/RhumbSolveTests.cs +++ b/GeographicLibTests/RhumbSolveTests.cs @@ -54,7 +54,7 @@ public void HiglyProlateOblateSphere( /// [DataTestMethod] [DataRow(0, 0, 0, 10001965, 89.99999347, 0.0, 0, false, DisplayName = "RhumbSovle8")] - [DataRow(0, 0, 60, 20003930, 89.99999347, 1654.69811, 1123576495931784.2, 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(