March 29th, 2016 NOTE: Since I’ve been dragging my heels at completing this I’m tossing it out as a WIP. In the process of making my guess to the origin of the sin based hashing function found in many bruteforce shaders I became curious if I could adapt one of my existing 2D hashing functions to the limitations of a WebGL based site like Shadertoy. The problem statement horseshoes and hand-grenades Hashing is probably a poor choice of terminology. Really the goal is a cheap approximation of sampling white noise for some input domain (integer lattice). This hashing function was built with the following design goals: Runtime evaluation (computationally cheap). The output should visually appears as white-noise and likewise for its power spectra approximation. Partially separable so the cost of hashing a neighborhood of points is cheaper than individual evaluations. Approximately uniform distribution. All other statistically properties are don’t care. The original version The 32-bit integer version I came up with looks like this: uint32_t hash(uint32_t x, uint32_t y) { x *= W0; // x' = Fx(x) y *= W1; // y' = Fy(y) x ^= y; // combine x *= M; // MLCG constant return x; } First a structural comment: Other than the XOR, all the mixing is from modular integer products so it suffers from the same problems as MLCGs 1. The short story is to use the top bits which will happen automatically if simply converting to floating point. There are logically three parts to the function: Take the input coordinate P=(x,y) and map it to new point P’=(x’,y’) such that neighbors of P generally do not map to points close to P’. Combine these with XOR which plays an important role. First it cheaply combines to two partial results while maintaining uniformity but more importantly it is using a different algebraic structure: vector addition in $\mathbb{F}_2$. Mixing groups is a common construct in modern hashing and PRNG methods. Perform a final mix with a well chosen MLCG constant. To perform part (1) I chose using independent Weyl sequences 2. Either values that together perform well at low-discrepancy such as from 3, which are specified as real values (0.5545497, 0.308517), they are scaled by 232 and rounded to odd. Or simply grabbing some square-free 4 and relatively prime (coprime) 5 values. To perform part (3) I chose a constant from Pierre L’Ecuyer’s classic paper 1. A possible set of constants are (with prime factors in the comments): #define W0 0x3504f333 // 3*2309*128413 #define W1 0xf1bbcdcb // 7*349*1660097 #define M 741103597 // 13*83*686843 What I didn’t do Since semi-random choices of constants was working well enough for me I didn’t bother to Do the Right Thing which would have been to perform a search for good values. A partial list of things we know: All the constants are odd to not reduce the domain (bijection). Given the symmetry of Weyl constants 2, only need to consider values with top-bit clear. Swapping the Weyl constants doesn’t change the properties (reflection about x=y). Weyl constants in a close neighborhood have similar properties. All the constants should be square-free and relatively prime. Property 5 is only of interest for finalizing a choice of constants. The probability that two randomly chosen odd integers are coprime is $ \frac{8}{\pi^2} \thickapprox .8106 $ 5. Coupled with property 4 if two share a prime factor then one can be tweaked by to a nearby odd value. Another thing is my originally stated design goal for Fx(x) and Fy(y) is overly strictly. We are not producing a 2D low-discrepancy sequence (LDS), we have 2 1D functions that need constants with some measure of linear independence. It is perfectly possible that one could be a small odd constant (like 1). The original reasoning for that contraint is as follows: If Fx & Fy behave as independent 1D LDSs, then when combined (by XOR) the 2D result should be LDS horizontally and vertical. Additionally if as a pair the constants produce a reason LDS then we should likewise be good along diagonals. Along the same lines for ISAs with a fair range of fast integer operations one of these mapping functions (Fx,Fy) could be replaced by something other than a Weyl sequence. Notably byteswap 6 or bit-rotation 7. Translating to shaders Using the high level language GLSL as example code. Then the original translates to: uint u32_weyl(in uvec2 c) { c = uvec2(W0,W1)*c; return M*(c.x^c.y); } which is fine for modern shaders and ISAs. STOP: Below here is only of potential interest if you cannot use XOR. (Update: since WebGL 2 is out and supported this is of approaching zero interest) To move toward the limited features of older shader models the first thing we have to change is to lose the XOR and the only obvious choice is multiplication. That brings up the following issues: No more mixed algebraic structure. The final mixing product (M) is now useless. Any value of M can be folded back into W0 and W1. Using the product will zero out results along the coordinate axes. The first two of these issues requires beefing up the initial 1D mappings. If they are made to be sufficiently independent and issue 3 is handled in some manner, then the product of two 1D white noise results yields a respectable 2D white noise. We can make the 1D mappings more pseudo-random like by moving to a nested Weyl sequence 2 and eliminate the zeroing out issue by adding in the opposite coordinate as follows: int i32_weyl(in ivec2 c) { c = ivec2(W0,W1)*c*c + ivec2(c.y,c.x); return c.x*c.y; } Which works reasonably well assuming the ISA has proper 32-bit integer support. If we knew the number of bits supported by the ISA, then we could build other bit-width versions. However web-based sites like Shadertoy do not AFAIK provide such information. In any case we can translate the previous into floating point like this: float f32_nweyl(vec2 c) { c = fract(c*fract(c*vec2(W0,W1)+vec2(c.y,c.x))); return fract(2048.0*c.x*c.y); } where the extra scaling factor (2048) is to throw away the top 12-bit. We want a uniform and not product distribution. By introducing the requirement that the input values never be zero (integer values offset by say 0.5), then we can remove the addition. In floating point we have input domain issues (unless we add extra operations) as the magnitude of the input increases. I made that even worse by tossing the final fract of (Fx,Fy) which also allows the removal of the extra scaling on the mix operation. Moving back to a scalar style, my final result looked like this: #define W0 0.5545497 #define W1 0.308517 float hash(in vec2 c) { float x = c.x*fract(c.x * W0); float y = c.y*fract(c.y * W1); return fract(x*y); } I have some usage examples with WebGL limitation (both integer and FP) on ShaderToy. If all of this portion sounds rather sloppy that’s because it was. I was more interest in the proof-of-concept than attempting to create a reasonable hashing function without XOR. Notice in my final version I don’t even make an attempt at creating reasonable constants, they are simply the unmodifed values pulled from 3. I “played around” a bit to amuse myself but didn’t do any real thinking here. There’s evidence that I’m doing this portion at least partially wrong, if you look at the integer triangle noise on ShaderToy the left hand side is a standard Weyl and not nested. It has some obvious directional defects but otherwise produces almost reasonable results. References and Footnotes “Tables of linear congruential generators of different sizes and good lattice structure”, Pierre L’Ecuyer, 1999 (PDF) ↩ ↩2 Weyl sequence overview post ↩ ↩2 ↩3 Mollwollfumble’s Science blog post: Subrandom Numbers ↩ ↩2 Square-free: an integer whose prime factorization has no repeats ↩ Relatively prime on Mathworld ↩ ↩2 32-bit single state xorshift with extra mixing ↩ “PCG: A Family of Simple Fast Space-Efficient Statistically Good Algorithms for Random Number Generation”, Melissa E. O’Neill, (paper) ↩ Comments math (32) Weyl (4) , low-discrepancy (3) , sampling (7)