1/**
2 * Shared WGSL fragments. Shaders are generated as strings with all sizes
3 * baked in as compile-time constants (the WGSL analog of what SHTNS does
4 * with NVRTC on CUDA: cf. init_cuda_program() in sht_gpu.cu).
5 *
6 * fp32 extended-range constants: same values SHTNS injects for a
7 * single-precision recurrence (sht_gpu.cu):
8 * SHT_ACCURACY = 1e-15
9 * SHT_SCALE_FACTOR = 2^56 = 7.2057594037927936e16
10 * A per-thread integer exponent `ny` counts how many times the running
11 * Legendre value has been multiplied by SCALE to stay in fp32 range;
12 * contributions are only accumulated once ny == 0 (value back in normal
13 * range and significant).
14 */
16export const RESCALE_WGSL = /* wgsl */ `
17const SCALE: f32 = 7.2057594e16; // rounds to exactly 2^56 in f32
18const INV_SCALE: f32 = 1.0 / 7.2057594e16;
19const ACCURACY: f32 = 1e-15;
20const RESCALE_THR: f32 = ACCURACY * SCALE + 1.0; // ~73: value became significant again
22struct Seed { y0: f32, ny: i32 }
24// Seed of the recurrence: y0 ~ sin(theta)^m by binary exponentiation with
25// rescaling (ports the HI_LLIM path of SHT/cuda_legendre.gen.cu, ~651-691).
26// The caller multiplies by amm afterwards (|amm| is O(1)).
27fn sinpow_rescaled(st: f32, m: u32) -> Seed {
28 var y0: f32 = 1.0;
29 var ny: i32 = 0;
30 if (m > 0u) {
31 var s: f32 = st;
32 var lb: u32 = m;
33 if ((lb & 1u) != 0u) { y0 = s; }
34 var nsint: i32 = 0;
35 lb = lb >> 1u;
36 while (lb > 0u) {
37 s = s * s;
38 nsint = nsint + nsint;
39 if (s < INV_SCALE) {
40 nsint = nsint - 1;
41 s = s * SCALE;
42 }
43 if ((lb & 1u) != 0u) {
44 y0 = y0 * s;
45 ny = ny + nsint;
46 if (y0 < (ACCURACY + INV_SCALE)) {
47 y0 = y0 * SCALE;
48 ny = ny - 1;
49 }
50 }
51 lb = lb >> 1u;
52 }
53 }
54 return Seed(y0, ny);
55}
56`;