Cubic Shader
Placeholder summary — update once the shader is written.
Source
image.glsl
// ============================================================
// CUBIC ALGEBRAIC NUMBERS IN THE UPPER HALF PLANE
//
// A cubic with leading coefficient a and root z factors as
//
// a * (x^2 - 2sx + n)(x - r)
//
// where s = Re(z), n = |z|^2. Fixing integer a, the lower
// three coefficients (b,c,d) trace a line parameterized by
// the real root r:
//
// L_a(r) = a(-2s, n, 0) + r * a(-1, 2s, -n)
//
// We loop over a = 1, ..., uAMax and raymarch each line in
// both directions. Only positive a: negating a negates all
// coefficients without changing the roots.
//
// Set uAMax = 1 for monic cubics only.
// ============================================================
// --- Search parameters ---
// uBallRadius, uRMax, uHitThreshold, uAMax and uIrreducibleOnly are
// uniforms from config.json. MAX_STEPS stays a compile-time constant:
// changing it is a recompile, not a control.
const int MAX_STEPS = 2048;
// --- View ---
//const vec2 VIEW_CENTER = vec2(-0.662359, 0.56228);
//const vec2 VIEW_CENTER = vec2(-0.419643, 0.606291);
const vec2 VIEW_CENTER = vec2(0., .75);
const float VIEW_HEIGHT = VIEW_CENTER.y*2.;
// ============================================================
// Lattice SDF
// ============================================================
float sdfLatticeBalls(vec3 p, float radius) {
return length(p - round(p)) - radius;
}
// ============================================================
// Irreducibility and primitivity
//
// A hit is skipped if the cubic is reducible over Q, or if it
// not primitive.
//
// Primitive: if b, c, d are all divisible by a divisor q > 1
// of a, the cubic is q times one already found at a/q.
//
// Reducible: the trace parameter r IS the real root of the
// cubic. A cubic over Q is reducible iff it has a rational
// root, and by the rational root theorem any such root has
// the form p/q with q | a. So for each divisor q of a, we
// round q*r to the nearest integer and evaluate the
// polynomial there. If it vanishes, the cubic is reducible.
// Uses float arithmetic to avoid int overflow; threshold 0.5
// is safe because any nonzero value of an integer polynomial
// at a rational point is at least 1/q^3 in absolute value.
// ============================================================
bool shouldSkip(int a, float r, ivec3 bcd) {
int b = bcd.x, c = bcd.y, d = bcd.z;
int absA = abs(a);
for (int q = 1; q <= absA; q++) {
if (absA % q != 0) continue;
if (uPrimitiveOnly && q > 1
&& b % q == 0 && c % q == 0 && d % q == 0) return true;
if (!uIrreducibleOnly) continue;
float fq = float(q);
float fn = round(fq * r);
float fval = float(a) * fn * fn * fn
+ float(b) * fn * fn * fq
+ float(c) * fn * fq * fq
+ float(d) * fq * fq * fq;
if (abs(fval) < 0.5) return true;
}
return false;
}
// ============================================================
// Line tracing
//
// March along base + t*dir for t >= 0.
// The trace parameter t = |r|; the sign argument recovers
// the actual real root r = sign * t.
// ============================================================
float traceHalf(vec3 base, vec3 dir, float tMax,
float radius, int a, float sign) {
float speed = length(dir);
float t = 0.0;
for (int i = 0; i < MAX_STEPS; i++) {
vec3 p = base + t * dir;
float d = sdfLatticeBalls(p, radius);
if (d < uHitThreshold) {
if (uIrreducibleOnly || uPrimitiveOnly) {
float r = sign * t;
ivec3 bcd = ivec3(round(p));
if (shouldSkip(a, r, bcd)) {
t += (radius + uHitThreshold)
/ speed;
continue;
}
}
return t;
}
t += d / speed;
if (t > tMax) { return 1000000.; }
}
return 1000000.;
}
vec2 traceLine(vec3 base, vec3 dir, float tMax,
float radius, int a) {
float tp = traceHalf(base, dir, tMax, radius, a, 1.0);
float tn = traceHalf(base, -dir, tMax, radius, a, -1.0);
return vec2(tp,tn);
}
// ============================================================
// Coordinate map
// ============================================================
vec2 pixelToHalfPlane(vec2 fragCoord, vec2 resolution, vec2 cen, float ht) {
vec2 uv = fragCoord / resolution;
float aspect = resolution.x / resolution.y;
return cen
+ (uv - 0.5)
* vec2(ht * aspect, ht);
}
// ============================================================
// Main
// ============================================================
void mainImage(out vec4 fragColor, in vec2 fragCoord) {
vec2 cen = VIEW_CENTER;
float ht = VIEW_HEIGHT/1.5;
vec2 z = pixelToHalfPlane(fragCoord, iResolution.xy, cen, ht);
float ball = uBallRadius*pow(ht,.33333)*pow(z.y,.3333333);
float maxR = uRMax/pow(ht,2.);
fragColor = vec4(vec3(1.), 1.0);
if (z.y <= 0.001) {
return;
}
float s = z.x;
float n = dot(z, z);
vec3 baseUnit = vec3(-2.0 * s, n, 0.0);
vec3 dirUnit = vec3(-1.0, 2.0 * s, -n);
// X is the strength of red hits (positive real root) at this pixel,
// Y of blue hits (negative real root). Each layer's weight is clamped
// to [0,1] (with |r| < 1 the raw 1/r exceeds 1), optionally scaled by
// 1/a^uLayerFalloff. Layers combine by max: a later layer replaces
// what came before only where its hit has a smaller real root, so
// each channel shows the single best cubic at that pixel.
float X = 0.;
float Y = 0.;
for (int a = 1; a <= uAMax; a++) {
// uAccumulate off restores the original behaviour: each a layer
// starts afresh, so only a == uAMax reaches the screen.
if (!uAccumulate) { X = 0.; Y = 0.; }
vec3 base = float(a) * baseUnit;
vec3 dir = float(a) * dirUnit;
vec2 r = traceLine(base, dir, maxR, ball, a);
// A miss returns a huge t from traceHalf; treat it as exactly zero
// weight, otherwise the tiny 1/t counts as a hit and blocks later
// layers.
float layer = pow(float(a), -uLayerFalloff);
float wx = r.x > maxR ? 0. : layer * clamp(1./(pow(ht,.333)*r.x), 0., 1.);
float wy = r.y > maxR ? 0. : layer * clamp(1./(pow(ht,.333)*r.y), 0., 1.);
X = max(X, wx);
Y = max(Y, wy);
}
// Pure red (X=1,Y=0) -> (1,0,0), pure blue (X=0,Y=1) -> (0,0,1), and
// both together -> purple (0.5,0,1). Green is removed by the strongest
// hit only and never accumulates. A red hit protects the blue channel
// fully; a blue hit lets red fall only as far as half.
fragColor = vec4(1. - Y * (1. - 0.5 * X),
1. - max(X, Y),
1. - X * (1. - Y),
1.);
}