diff --git a/types/geo/point.go b/types/geo/point.go index 820582b0f..d039ea1fa 100644 --- a/types/geo/point.go +++ b/types/geo/point.go @@ -125,6 +125,9 @@ func (p Point) SphericalAngleTo(q Point) (Radians, error) { sLat, sLng := float64(qLat.Radians()), float64(qLng.Radians()) cosA := math.Sin(rLat)*math.Sin(sLat) + math.Cos(rLat)*math.Cos(sLat)*math.Cos(rLng-sLng) + // Subtle floating point imprecision can lead to cosA being outside + // the domain of arccosine [-1, 1]. Clamp the input to avoid NaN return. + cosA = min(max(-1.0, cosA), 1.0) return Radians(math.Acos(cosA)), nil } diff --git a/types/geo/point_test.go b/types/geo/point_test.go index f0d0cb3ab..32a73180a 100644 --- a/types/geo/point_test.go +++ b/types/geo/point_test.go @@ -448,65 +448,79 @@ func TestPointMarshalUint64(t *testing.T) { }) } +const earthRadius = 6371.000 // volumetric mean radius (km) +const kmToRad = 1 / earthRadius + +// Test corpus for exercising PointSphericalAngleTo. +var corpusPointSphericalAngleTo = []struct { + name string + x geo.Point + y geo.Point + want geo.Radians + wantErr string +}{ + { + name: "same-point-null-island", + x: geo.MakePoint(0, 0), + y: geo.MakePoint(0, 0), + want: 0.0 * geo.Radian, + }, + { + name: "same-point-north-pole", + x: geo.MakePoint(+90, 0), + y: geo.MakePoint(+90, +90), + want: 0.0 * geo.Radian, + }, + { + name: "same-point-south-pole", + x: geo.MakePoint(-90, 0), + y: geo.MakePoint(-90, -90), + want: 0.0 * geo.Radian, + }, + { + name: "north-pole-to-south-pole", + x: geo.MakePoint(+90, 0), + y: geo.MakePoint(-90, -90), + want: math.Pi * geo.Radian, + }, + { + name: "toronto-to-montreal", + x: geo.MakePoint(+43.6532, -79.3832), + y: geo.MakePoint(+45.5019, -73.5674), + want: 504.26 * kmToRad * geo.Radian, + }, + { + name: "sydney-to-san-francisco", + x: geo.MakePoint(-33.8727, +151.2057), + y: geo.MakePoint(+37.7749, -122.4194), + want: 11948.18 * kmToRad * geo.Radian, + }, + { + name: "new-york-to-paris", + x: geo.MakePoint(+40.7128, -74.0060), + y: geo.MakePoint(+48.8575, +2.3514), + want: 5837.15 * kmToRad * geo.Radian, + }, + { + name: "seattle-to-tokyo", + x: geo.MakePoint(+47.6061, -122.3328), + y: geo.MakePoint(+35.6764, +139.6500), + want: 7700.00 * kmToRad * geo.Radian, + }, + { + // Subtle floating point imprecision can propagate and lead to + // trigonometric functions receiving inputs outside their + // domain, thus returning NaN. + // Test one such case. + name: "floating-point-precision-test", + x: geo.MakePoint(-6.0, 0.0), + y: geo.MakePoint(-6.0, 0.0), + want: 0.0 * geo.Radian, + }, +} + func TestPointSphericalAngleTo(t *testing.T) { - const earthRadius = 6371.000 // volumetric mean radius (km) - const kmToRad = 1 / earthRadius - for _, tt := range []struct { - name string - x geo.Point - y geo.Point - want geo.Radians - wantErr string - }{ - { - name: "same-point-null-island", - x: geo.MakePoint(0, 0), - y: geo.MakePoint(0, 0), - want: 0.0 * geo.Radian, - }, - { - name: "same-point-north-pole", - x: geo.MakePoint(+90, 0), - y: geo.MakePoint(+90, +90), - want: 0.0 * geo.Radian, - }, - { - name: "same-point-south-pole", - x: geo.MakePoint(-90, 0), - y: geo.MakePoint(-90, -90), - want: 0.0 * geo.Radian, - }, - { - name: "north-pole-to-south-pole", - x: geo.MakePoint(+90, 0), - y: geo.MakePoint(-90, -90), - want: math.Pi * geo.Radian, - }, - { - name: "toronto-to-montreal", - x: geo.MakePoint(+43.6532, -79.3832), - y: geo.MakePoint(+45.5019, -73.5674), - want: 504.26 * kmToRad * geo.Radian, - }, - { - name: "sydney-to-san-francisco", - x: geo.MakePoint(-33.8727, +151.2057), - y: geo.MakePoint(+37.7749, -122.4194), - want: 11948.18 * kmToRad * geo.Radian, - }, - { - name: "new-york-to-paris", - x: geo.MakePoint(+40.7128, -74.0060), - y: geo.MakePoint(+48.8575, +2.3514), - want: 5837.15 * kmToRad * geo.Radian, - }, - { - name: "seattle-to-tokyo", - x: geo.MakePoint(+47.6061, -122.3328), - y: geo.MakePoint(+35.6764, +139.6500), - want: 7700.00 * kmToRad * geo.Radian, - }, - } { + for _, tt := range corpusPointSphericalAngleTo { t.Run(tt.name, func(t *testing.T) { got, err := tt.x.SphericalAngleTo(tt.y) if tt.wantErr == "" && err != nil { @@ -536,6 +550,23 @@ func TestPointSphericalAngleTo(t *testing.T) { } } +func FuzzPointSphericalAngleTo(f *testing.F) { + for _, tt := range corpusPointSphericalAngleTo { + xLat, xLng, _ := tt.x.LatLngFloat64() + yLat, yLng, _ := tt.y.LatLngFloat64() + f.Add(xLat, xLng, yLat, yLng) + } + + f.Fuzz(func(t *testing.T, xLat float64, xLng float64, yLat float64, yLng float64) { + x := geo.MakePoint(geo.Degrees(xLat), geo.Degrees(xLng)) + y := geo.MakePoint(geo.Degrees(yLat), geo.Degrees(yLng)) + got, _ := x.SphericalAngleTo(y) + if math.IsNaN(float64(got)) { + t.Errorf("got NaN result with xLat=%.15f xLng=%.15f yLat=%.15f yLng=%.15f", xLat, xLng, yLat, yLng) + } + }) +} + func approx[T ~float64](x, y T) bool { return math.Abs(float64(x)-float64(y)) <= 1e-5 }