leg_synth: advance the recurrence the way leg_analys does

scripts/diagnose-leg.ts on an RTX PRO 6000 (driver 590.48, Blackwell) follows the m = 0 recurrence out of the production shader term by term: l ilat 0 ilat 1 ilat 64 ilat 127 0 ok ok ok ok 1 ok ok ok ok 2 x0.00e+0 x0.00e+0 x0.00e+0 x0.00e+0 3 x-6.67e-1 x-6.70e-1 x4.45e-1 x-6.67e-1 y_2 is not drifting, it is exactly zero at every latitude, and y_3 comes to -6.67e-1 of its value — which is exactly b_3*y_1/y_3, the right coefficients applied to y_2 = 0. So on the first iteration let c0 = ab[base + (l + 2u - m)]; y0 = c0.x * ct * y1 + c0.y * y0; produces zero, while the c1 read two lines later is fine. The buffer is not at fault: leg_analys reads the same ab and is correct to 1e-7 on the same device, and m = 62 and 63 — the only orders whose loop breaks before this runs — are the only correct ones. It also explains why the near-equator latitude stays close to right while the poles reach 1e11: there ct is ~0.012, so the term that goes missing barely contributes. Nothing in that form is invalid WGSL, so this is a miscompiled load. The change makes the advance structurally identical to leg_analys', which the same driver compiles correctly: both coefficients fetched unconditionally, and the new y0 carried in a temporary rather than assigned and then read back by the y1 update. Two shaders doing the same recurrence should agree on how anyway. NOT YET CONFIRMED to fix that machine — there is no Nvidia GPU here. It does not regress Intel Xe via Dawn, where test:node passed before and passes now, and the reasoning above is what it is based on.

Jeremy Magland <jmagland@flatironinstitute.org> committed 20 days ago commit f290310ff72f parent 5e5b2f8 Browse files

1 changed file+22−4

src/sht/wgsl/leg.tsmodified+22−4View file

@@ -73,12 +73,30 @@ fn leg_synth(@builtin(global_invocation_id) gid: vec3u,

7373 y1 *= INV_SCALE;

7474 }

7575 if (l + 2u > LMAX) { break; }

76- let c0 = ab[base + (l + 2u - m)];

77- y0 = c0.x * ct * y1 + c0.y * y0;

76+ // Advance (y_l, y_{l+1}) to (y_{l+2}, y_{l+3}).

77+ //

78+ // Written in exactly the shape leg_analys uses below — both coefficients

79+ // fetched unconditionally, the new y0 carried in a temporary rather than

80+ // assigned and then read back by the y1 update. The shorter form,

81+ //

82+ // let c0 = ab[base + (l + 2u - m)];

83+ // y0 = c0.x * ct * y1 + c0.y * y0;

84+ // if (l + 3u <= LMAX) { ... y1 = c1.x * ct * y0 + c1.y * y1; }

85+ //

86+ // says the same thing and is what this was, but NVIDIA's Vulkan compiler

87+ // (driver 590.48, Blackwell) mis-compiles it: c0 reads as (0, 0) on the

88+ // first iteration, so y_{l+2} comes out exactly zero and every later term

89+ // follows a different solution of the recurrence, reaching ~1e11 by l = 63.

90+ // leg_analys, doing the same arithmetic in this shape, was correct on the

91+ // same driver. See scripts/diagnose-leg.ts, which is how that was found.

92+ let a0 = ab[base + (l + 2u - m)];

93+ var a1 = vec2f(0.0);

7894 if (l + 3u <= LMAX) {

79- let c1 = ab[base + (l + 3u - m)];

80- y1 = c1.x * ct * y0 + c1.y * y1;

95+ a1 = ab[base + (l + 3u - m)];

8196 }

97+ let t0 = a0.x * ct * y1 + a0.y * y0;

98+ y1 = a1.x * ct * t0 + a1.y * y1;

99+ y0 = t0;

82100 l += 2u;

83101 }

84102 fm[m * NLAT + ilat] = acc;