// Shared simulation math; see docs/determinism.md. Each routine uses only
// correctly rounded +, -, * and integer bit operations in a fixed order, so
// it returns the same bits on every device. `precise` keeps the compiler
// from fusing or reordering them; it only covers the function it is written
// in, so every intermediate here is a `precise` local.
// src/runtime/sim_math.rs is the bit-exact Rust reference; change both.
// ponytail: these run in every DeterminismMode, Off included; switch them to
// native built-ins through a specialization constant if profiling shows the
// cost.
#ifndef SIM_MATH_GLSL
#define SIM_MATH_GLSL
// Newton steps approach from below and stop an ULP short of round values
// such as 1 / 1. Each routine ends by keeping whichever of y - 1 ULP, y,
// and y + 1 ULP leaves the smallest residual.
float sim_neighbour(float y, int offset) {
return uintBitsToFloat(uint(int(floatBitsToUint(y)) + offset));
}
// 1 / x for finite, normal, non-zero x, within 1 ULP.
float sim_recip(float x) {
float magnitude = abs(x);
// Bit-trick seed, then three Newton steps y = y * (2 - x * y).
precise float y = uintBitsToFloat(0x7EF311C3u - floatBitsToUint(magnitude));
for (int step = 0; step < 3; step++) {
precise float error = 2.0 - magnitude * y;
y = y * error;
}
float best = y;
precise float best_residual = abs(1.0 - magnitude * y);
for (int offset = -1; offset <= 1; offset += 2) {
float candidate = sim_neighbour(y, offset);
precise float residual = abs(1.0 - magnitude * candidate);
if (residual < best_residual) {
best = candidate;
best_residual = residual;
}
}
return x < 0.0 ? -best : best;
}
float sim_div(float a, float b) {
precise float result = a * sim_recip(b);
return result;
}
vec3 sim_div(vec3 v, float s) {
precise vec3 result = v * sim_recip(s);
return result;
}
// 1 / sqrt(x) for finite, normal, positive x, within 2 ULP.
float sim_rsqrt(float x) {
precise float half_x = 0.5 * x;
precise float y = uintBitsToFloat(0x5F375A86u - (floatBitsToUint(x) >> 1));
for (int step = 0; step < 3; step++) {
precise float square = half_x * y;
square = square * y;
precise float error = 1.5 - square;
y = y * error;
}
float best = y;
precise float best_residual = abs(1.0 - x * y * y);
for (int offset = -1; offset <= 1; offset += 2) {
float candidate = sim_neighbour(y, offset);
precise float residual = abs(1.0 - x * candidate * candidate);
if (residual < best_residual) {
best = candidate;
best_residual = residual;
}
}
return best;
}
// sqrt(x); 0 for x <= 0.
float sim_sqrt(float x) {
if (x <= 0.0) return 0.0;
precise float root = x * sim_rsqrt(x);
float best = root;
precise float best_residual = abs(root * root - x);
for (int offset = -1; offset <= 1; offset += 2) {
float candidate = sim_neighbour(root, offset);
precise float residual = abs(candidate * candidate - x);
if (residual < best_residual) {
best = candidate;
best_residual = residual;
}
}
return best;
}
float sim_dot(vec3 a, vec3 b) {
precise float sum = a.x * b.x;
sum = sum + a.y * b.y;
sum = sum + a.z * b.z;
return sum;
}
float sim_length(vec3 v) {
return sim_sqrt(sim_dot(v, v));
}
// v scaled to length 1; the zero vector stays zero.
vec3 sim_normalize(vec3 v) {
float square = sim_dot(v, v);
if (square <= 0.0) return vec3(0.0);
precise vec3 result = v * sim_rsqrt(square);
return result;
}
// m * v with each row summed in x, y, z order.
vec3 sim_mul(mat3 m, vec3 v) {
return vec3(
sim_dot(vec3(m[0].x, m[1].x, m[2].x), v),
sim_dot(vec3(m[0].y, m[1].y, m[2].y), v),
sim_dot(vec3(m[0].z, m[1].z, m[2].z), v)
);
}
// transpose(m) * v: each column dotted with v.
vec3 sim_mul_transposed(mat3 m, vec3 v) {
return vec3(sim_dot(m[0], v), sim_dot(m[1], v), sim_dot(m[2], v));
}
// float to int the way Rust's `as i32` does it: truncate toward zero,
// saturate out-of-range values, NaN to 0. GLSL's int() leaves those cases
// undefined.
int sim_to_int(float x) {
if (isnan(x)) return 0;
if (x >= 2147483648.0) return 2147483647;
if (x <= -2147483648.0) return int(0x80000000u);
return int(x);
}
ivec3 sim_to_int(vec3 v) {
return ivec3(sim_to_int(v.x), sim_to_int(v.y), sim_to_int(v.z));
}
#endif