diff --git a/CHANGELOG.md b/CHANGELOG.md index 0b0b0da1..a38e02b4 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/internal/rf/gpu/bearingprecision_test.go b/internal/rf/gpu/bearingprecision_test.go new file mode 100644 index 00000000..bb454874 --- /dev/null +++ b/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 +} diff --git a/internal/rf/gpu/coverage.wgsl b/internal/rf/gpu/coverage.wgsl index e9f87479..4482db3c 100644 --- a/internal/rf/gpu/coverage.wgsl +++ b/internal/rf/gpu/coverage.wgsl @@ -61,12 +61,16 @@ fn height_at(lat: f32, lon: f32) -> vec2 { return vec2(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) { 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; diff --git a/internal/rf/gpu/coveragefold_equiv_test.go b/internal/rf/gpu/coveragefold_equiv_test.go index 52a22d1f..aa781510 100644 --- a/internal/rf/gpu/coveragefold_equiv_test.go +++ b/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)) }