GPU coverage fold: compute the bearing without the driver's trig

CHANGELOG.md

@@ -65,6 +65,11 @@ had changed in them - which is the gap this file exists to close.
   opened on a library of builds that looked unpublished. A row whose half of
   the catalogue is unreachable now offers a retry. `firmware.library` rows
   carry `checking` and `catalogue_down`.
+- **GPU coverage no longer depends on the driver's trig.** The fold kernel's
+  bearing cancelled two near-equal numbers, so a driver whose `sin` and `cos`
+  are only as accurate as Vulkan requires, as on an Intel UHD 620, put a
+  directional antenna's gain 0.9 dB out. The bearing, distance and elevation
+  now use the kernel's own series.
 
 ## [0.3.0] - 2026-10-24
 

internal/rf/gpu/bearingprecision_test.go

@@ -0,0 +1,166 @@
+package gpu
+
+import (
+	"math"
+	"testing"
+
+	"github.com/MeshBench/meshbench/internal/rf/geo"
+)
+
+// The fold kernel's bearing, carried out here in float32 with the worst trig
+// Vulkan allows - sin and cos each out by 2^-11 - holds to the CPU twin.
+//
+// No CI runner has a GPU, so the kernel itself is only ever compared on
+// whatever device someone owns. This is the arithmetic of coverage.wgsl's
+// bearing_deg, psin, pcos, psin_lat and atan2p written out in Go, so the
+// property can be checked everywhere: on an Intel UHD 620 the textbook form
+// put a beam's bearing nearly nine degrees out, because its x term cancels
+// two numbers near 0.46 down to about a thousandth and the trig error
+// survives whole.
+func TestTheFoldBearingSurvivesTheLeastAccurateTrig(t *testing.T) {
+	const stLat, stLon = 56.6, -4.0 // the directional station in the fold comparison
+	worstNew, worstOld := 0.0, 0.0
+	for _, sinErr := range []float32{-trigErr, trigErr} {
+		for _, cosErr := range []float32{-trigErr, trigErr} {
+			sinD := func(x float32) float32 { return float32(math.Sin(float64(x))) + sinErr }
+			cosD := func(x float32) float32 { return float32(math.Cos(float64(x))) + cosErr }
+			for r := range 96 {
+				for c := range 96 {
+					lat := 56.5 + (float64(r)+0.5)*0.5/96
+					lon := -4.2 + (float64(c)+0.5)*1.0/96
+					if math.Abs(lat-stLat) < 1e-9 && math.Abs(lon-stLon) < 1e-9 {
+						continue
+					}
+					want := geo.BearingDeg(stLat, stLon, lat, lon)
+					worstNew = math.Max(worstNew, angleDiff(want,
+						float64(bearingNew32(stLat, stLon, float32(lat), float32(lon), sinD, cosD))))
+					worstOld = math.Max(worstOld, angleDiff(want,
+						float64(bearingOld32(stLat, stLon, float32(lat), float32(lon), sinD, cosD))))
+				}
+			}
+		}
+	}
+	t.Logf("worst bearing error with trig out by 2^-11: textbook form %.3f deg, kernel's form %.4f deg",
+		worstOld, worstNew)
+	// What is left is float32 coordinates, not trig: a cell a quarter of a
+	// kilometre from the mast is a few float32 steps of latitude away, and
+	// that alone is 0.03 degrees, falling to under 0.01 by 2 km. On the fold
+	// test's beam, 12 dB over 180 degrees, 0.05 degrees is 0.005 dB against a
+	// tolerance of 0.05.
+	if worstNew > 0.05 {
+		t.Errorf("the kernel's bearing is %.4f degrees out with the least accurate trig Vulkan permits", worstNew)
+	}
+}
+
+// trigErr is the absolute error Vulkan permits a driver's sin and cos.
+const trigErr = 1.0 / 2048
+
+func angleDiff(a, b float64) float64 {
+	d := math.Mod(math.Abs(a-b), 360)
+	return math.Min(d, 360-d)
+}
+
+// series is coverage.wgsl's psin, pcos and psin_lat, over a driver's sin.
+type series struct{ sinD func(float32) float32 }
+
+func (k series) psin(x float32) float32 {
+	if x > math.Pi/4 || x < -math.Pi/4 {
+		return k.sinD(x)
+	}
+	x2 := x * x
+	return x * (1 - x2/6*(1-x2/20*(1-x2/42*(1-x2/72*(1-x2/110)))))
+}
+
+func (k series) pcos(x float32) float32 {
+	s := k.psin(x * 0.5)
+	return 1 - 2*s*s
+}
+
+func (k series) psinLat(x float32) float32 { return 2 * k.psin(x*0.5) * k.pcos(x*0.5) }
+
+// bearingNew32 is coverage.wgsl's bearing_deg.
+func bearingNew32(lat1, lon1, lat2, lon2 float32, sinD, _ func(float32) float32) float32 {
+	const rad = float32(math.Pi / 180)
+	k := series{sinD}
+	dLat := (lat2 - lat1) * rad
+	dLon := (lon2 - lon1) * rad
+	c2 := k.pcos(lat2 * rad)
+	h := k.psin(dLon * 0.5)
+	y := k.psin(dLon) * c2
+	x := k.psin(dLat) + 2*k.psinLat(lat1*rad)*c2*h*h
+	bg := atan2p32(y, x) / rad
+	if bg < 0 {
+		bg += 360
+	}
+	return bg
+}
+
+// haversine32 is coverage.wgsl's haversine_km.
+func haversine32(lat1, lon1, lat2, lon2 float32, sinD func(float32) float32) float32 {
+	const rad = float32(math.Pi / 180)
+	k := series{sinD}
+	sLat := k.psin((lat2 - lat1) * rad / 2)
+	sLon := k.psin((lon2 - lon1) * rad / 2)
+	a := sLat*sLat + k.pcos(lat1*rad)*k.pcos(lat2*rad)*sLon*sLon
+	return 2 * 6371 * float32(math.Asin(math.Min(1, math.Sqrt(float64(a)))))
+}
+
+// The kernel's distance holds too. Half a cell's latitude difference is a few
+// ten-thousandths of a radian, the size of the error Vulkan permits sin, and
+// with the driver's trig the old form put a cell 50 km away 7 km out.
+func TestTheFoldDistanceSurvivesTheLeastAccurateTrig(t *testing.T) {
+	const stLat, stLon = 56.6, -4.0
+	worst := 0.0
+	for _, sinErr := range []float32{-trigErr, trigErr} {
+		sinD := func(x float32) float32 { return float32(math.Sin(float64(x))) + sinErr }
+		for r := range 96 {
+			for c := range 96 {
+				lat := 56.5 + (float64(r)+0.5)*0.5/96
+				lon := -4.2 + (float64(c)+0.5)*1.0/96
+				want := geo.DistanceKm(stLat, stLon, lat, lon)
+				got := float64(haversine32(stLat, stLon, float32(lat), float32(lon), sinD))
+				worst = math.Max(worst, math.Abs(got-want))
+			}
+		}
+	}
+	t.Logf("worst distance error with trig out by 2^-11: %.2f m", worst*1000)
+	// Float32 coordinates again: a few metres near the mast, a few hundredths
+	// of a dB of path loss at worst.
+	if worst > 0.01 {
+		t.Errorf("the kernel's distance is %.1f m out with the least accurate trig Vulkan permits", worst*1000)
+	}
+}
+
+// bearingOld32 is the textbook form the kernel used before.
+func bearingOld32(lat1, lon1, lat2, lon2 float32, sinD, cosD func(float32) float32) float32 {
+	const rad = float32(math.Pi / 180)
+	y := sinD((lon2-lon1)*rad) * cosD(lat2*rad)
+	x := cosD(lat1*rad)*sinD(lat2*rad) - sinD(lat1*rad)*cosD(lat2*rad)*cosD((lon2-lon1)*rad)
+	bg := float32(math.Atan2(float64(y), float64(x))) / rad
+	if bg < 0 {
+		bg += 360
+	}
+	return bg
+}
+
+// atan2p32 is coverage.wgsl's atan2p.
+func atan2p32(y, x float32) float32 {
+	ax, ay := float32(math.Abs(float64(x))), float32(math.Abs(float64(y)))
+	mx := max(ax, ay)
+	if mx == 0 {
+		return 0
+	}
+	z := min(ax, ay) / mx
+	z2 := z * z
+	a := z * (0.9998660 + z2*(-0.3302995+z2*(0.1801410+z2*(-0.0851330+z2*0.0208351))))
+	if ay > ax {
+		a = math.Pi/2 - a
+	}
+	if x < 0 {
+		a = math.Pi - a
+	}
+	if y < 0 {
+		a = -a
+	}
+	return a
+}

internal/rf/gpu/coverage.wgsl

@@ -61,12 +61,16 @@ fn height_at(lat: f32, lon: f32) -> vec2<f32> {
   return vec2<f32>(top * (1.0 - ty) + bot * ty, 1.0);
 }
 
+// haversine_km uses psin and pcos (below) for the reason bearing_deg does:
+// half a cell's latitude difference is a few ten-thousandths of a radian, no
+// bigger than the error Vulkan lets a driver's sin have.
 fn haversine_km(lat1: f32, lon1: f32, lat2: f32, lon2: f32) -> f32 {
   let rad = PI / 180.0;
   let d_lat = (lat2 - lat1) * rad;
   let d_lon = (lon2 - lon1) * rad;
-  let a = sin(d_lat / 2.0) * sin(d_lat / 2.0) +
-    cos(lat1 * rad) * cos(lat2 * rad) * sin(d_lon / 2.0) * sin(d_lon / 2.0);
+  let s_lat = psin(d_lat / 2.0);
+  let s_lon = psin(d_lon / 2.0);
+  let a = s_lat * s_lat + pcos(lat1 * rad) * pcos(lat2 * rad) * s_lon * s_lon;
   return 2.0 * 6371.0 * asin(min(1.0, sqrt(a)));
 }
 
@@ -219,12 +223,76 @@ fn gain_db(bearing: f32, elev: f32) -> f32 {
   return g0 * (1.0 - te) + g1 * te;
 }
 
+// psin is sin from its series, for |x| up to pi/4, and pcos and psin_lat
+// build on it for any latitude. The bearing uses no other trig.
+//
+// Vulkan asks a driver's sin and cos for an absolute error of 2^-11 and no
+// better, and some meet exactly that. On an angle of a few thousandths of a
+// radian that is most of the answer, and on a latitude it is still a tenth of
+// a percent - several hundredths of a degree of bearing. The series is good to
+// float32 across its range on every device.
+fn psin(x: f32) -> f32 {
+  if (abs(x) > PI / 4.0) {
+    return sin(x);
+  }
+  let x2 = x * x;
+  return x * (1.0 - x2 / 6.0 * (1.0 - x2 / 20.0 * (1.0 - x2 / 42.0 * (1.0 - x2 / 72.0 * (1.0 - x2 / 110.0)))));
+}
+
+// pcos is cos for |x| up to pi/2, by the half angle.
+fn pcos(x: f32) -> f32 {
+  let s = psin(x * 0.5);
+  return 1.0 - 2.0 * s * s;
+}
+
+// psin_lat is sin for |x| up to pi/2, by the half angle.
+fn psin_lat(x: f32) -> f32 {
+  return 2.0 * psin(x * 0.5) * pcos(x * 0.5);
+}
+
+// atan2p is atan2 to about 1e-5 rad on every device (Abramowitz and Stegun
+// 4.4.49 on [0, 1], then the octants), for the same reason: the driver's own
+// is only promised to a few thousand ulp.
+fn atan2p(y: f32, x: f32) -> f32 {
+  let ax = abs(x);
+  let ay = abs(y);
+  let mx = max(ax, ay);
+  if (mx == 0.0) {
+    return 0.0;
+  }
+  let z = min(ax, ay) / mx;
+  let z2 = z * z;
+  var a = z * (0.9998660 + z2 * (-0.3302995 + z2 * (0.1801410 + z2 * (-0.0851330 + z2 * 0.0208351))));
+  if (ay > ax) {
+    a = PI / 2.0 - a;
+  }
+  if (x < 0.0) {
+    a = PI - a;
+  }
+  if (y < 0.0) {
+    a = -a;
+  }
+  return a;
+}
+
+// bearing_deg is the initial great-circle bearing, as geo.BearingDeg.
+//
+// Not in the textbook form. Its x term, cos(lat1)sin(lat2) -
+// sin(lat1)cos(lat2)cos(dlon), is two numbers near 0.46 that differ by about
+// a thousandth for a cell a few kilometres away, so the driver's trig error
+// lands on the difference whole. On an Intel UHD 620 that put a beam's
+// bearing nearly nine degrees out and its gain 0.9 dB wrong. The same
+// quantity is sin(dlat) + 2 sin(lat1) cos(lat2) sin^2(dlon/2): nothing
+// cancels, and every sine is the series above.
 fn bearing_deg(lat1: f32, lon1: f32, lat2: f32, lon2: f32) -> f32 {
   let rad = PI / 180.0;
-  let y = sin((lon2 - lon1) * rad) * cos(lat2 * rad);
-  let x = cos(lat1 * rad) * sin(lat2 * rad) -
-    sin(lat1 * rad) * cos(lat2 * rad) * cos((lon2 - lon1) * rad);
-  var bg = atan2(y, x) / rad;
+  let d_lat = (lat2 - lat1) * rad;
+  let d_lon = (lon2 - lon1) * rad;
+  let c2 = pcos(lat2 * rad);
+  let h = psin(d_lon * 0.5);
+  let y = psin(d_lon) * c2;
+  let x = psin(d_lat) + 2.0 * psin_lat(lat1 * rad) * c2 * h * h;
+  var bg = atan2p(y, x) / rad;
   if (bg < 0.0) {
     bg += 360.0;
   }
@@ -271,7 +339,7 @@ fn fold(@builtin(global_invocation_id) gid: vec3<u32>) {
     m_in = b.rem_tx + b.rem_gain - self_loss + self_gain - b.st_sens;
   } else {
     let bearing = bearing_deg(p.st_lat, p.st_lon, ll.x, ll.y);
-    let elev = atan2(rx_alt - p.st_alt, dist_km * 1000.0) * 180.0 / PI;
+    let elev = atan2p(rx_alt - p.st_alt, dist_km * 1000.0) * 180.0 / PI;
     let g = gain_db(bearing, elev);
     m_out = b.st_tx + g - loss + b.rem_gain - b.rem_sens;
     m_in = b.rem_tx + b.rem_gain - loss + g - b.st_sens;

internal/rf/gpu/coveragefold_equiv_test.go

@@ -75,6 +75,7 @@ func compareCoverageFoldCPU(t *testing.T) {
 			t.Fatalf("%s: GPU returned %d slots, CPU %d", name, len(got), len(want))
 		}
 		wrong := 0
+		worst := 0.0
 		for i := range want {
 			for _, v := range []float32{got[i].MinDB, got[i].OutDB, got[i].InDB} {
 				if math.IsNaN(float64(v)) || math.IsInf(float64(v), 0) {
@@ -89,12 +90,16 @@ func compareCoverageFoldCPU(t *testing.T) {
 				}
 				continue
 			}
+			worst = math.Max(worst, math.Abs(float64(got[i].MinDB-want[i].MinDB)))
 			if math.Abs(float64(got[i].MinDB-want[i].MinDB)) > 0.05 ||
 				math.Abs(float64(got[i].OutDB-want[i].OutDB)) > 0.05 ||
 				math.Abs(float64(got[i].InDB-want[i].InDB)) > 0.05 {
 				wrong++
 			}
 		}
+		// Said on a pass too: a device that agrees with 0.049 dB to spare of
+		// a 0.05 dB tolerance is one driver update from failing.
+		t.Logf("%s: worst margin disagreement %.4f dB over %d cells", name, worst, len(want))
 		if wrong > 0 {
 			t.Fatalf("%s: %d of %d cells disagree with the CPU twin", name, wrong, len(want))
 		}