concept-collection / turing-surface-cache
turing-surface-cache / src / mgpu / noise.ts
51 lines · 1.5 KBBlameHistoryRaw
1/**
2 * The seeded perturbation a model's `init` starts from.
3 *
4 * Host-side rather than in the .m, so a run is reproducible from an integer
5 * seed and the same field can be handed to any model.
6 */
8/**
9 * Seeded uniform deviates in [0, 1): mulberry32.
10 *
11 * Integer arithmetic and one division by 2^32, so any faithful port of it
12 * produces bit-identical values — which is what lets the native benchmark under
13 * bench/shtns/ seed the same run.
14 */
15export function makeRand(seed: number): () => number {
16 let s = seed >>> 0;
17 return (): number => {
18 s = (s + 0x6d2b79f5) >>> 0;
19 let t = s;
20 t = Math.imul(t ^ (t >>> 15), t | 1);
21 t ^= t + Math.imul(t ^ (t >>> 7), t | 61);
22 return ((t ^ (t >>> 14)) >>> 0) / 4294967296;
23 };
26/** Seeded normal deviates: mulberry32 + Box-Muller. */
27export function makeRandn(seed: number): () => number {
28 const rand = makeRand(seed);
29 let spare: number | null = null;
30 return () => {
31 if (spare !== null) {
32 const v = spare;
33 spare = null;
34 return v;
35 }
36 let u = 0;
37 while (u === 0) u = rand();
38 const r = Math.sqrt(-2 * Math.log(u));
39 const th = 2 * Math.PI * rand();
40 spare = r * Math.sin(th);
41 return r * Math.cos(th);
42 };
45/** `amp`-scaled normal deviates, one per grid point, in index order. */
46export function seededNoise(npts: number, amp: number, seed: number): Float32Array {
47 const randn = makeRandn(seed);
48 const out = new Float32Array(npts);
49 for (let i = 0; i < npts; i++) out[i] = amp * randn();
50 return out;