mosprites at compile time

mosprite.hpp

// mosprite.hpp — compile-time sprite synthesis and n-dimensional rasterisation.
//
// A C++20 transcription of the MoSprites NG pipeline (mosprites-ng.md) in which
// every stage parameter is a type. That is not decoration. Three things follow
// from moving the policy into the type system, and only the first is about speed:
//
//   1. The inner loop contains no branch on operator, order, or colour space.
//      `if constexpr` deletes the untaken paths outright.
//
//   2. Dimension is a template parameter. ng §9 argues that the whole apparatus
//      lifts to n dimensions because coverage is separable, a = product of f_i.
//      Here that argument *is* the implementation: Walk<N,Axis> unrolls the loop
//      nest, and the 2-D sprite and the 3-D voxel model are one function.
//
//   3. Scale is a compile-time rational, so q = denom(p) is a constant, and the
//      symmetry results derived in aamazing.md §9 and corrected in ng §0 become
//      predicates the compiler can check. A configuration that claims a mirror
//      it does not have fails to compile. See `require_mirror`.
//
// Requires C++20. Header-only, no dependencies.

#pragma once
#include <array>
#include <cstdint>
#include <vector>
#include <string_view>

namespace mosprite {

// ─────────────────────────────────────────────────────── constexpr arithmetic
// Written out rather than pulled from <cmath> so every path is usable in a
// constant expression on any conforming C++20 implementation.

constexpr int  cgcd(int a, int b) { return b ? cgcd(b, a % b) : (a < 0 ? -a : a); }
constexpr int  cfloor(double x)   { int i = int(x); return (x < 0 && double(i) != x) ? i - 1 : i; }
constexpr int  cceil (double x)   { int i = int(x); return (x > 0 && double(i) != x) ? i + 1 : i; }
constexpr double cabs(double x)   { return x < 0 ? -x : x; }
constexpr double cmin(double a, double b) { return a < b ? a : b; }
constexpr int    imin(int a, int b) { return a < b ? a : b; }
constexpr int    imax(int a, int b) { return a > b ? a : b; }
constexpr double cmax(double a, double b) { return a > b ? a : b; }

// overlap of [a0,a1) with [b0,b1) — the 1-D factor of a separable coverage
constexpr double span(double a0, double a1, double b0, double b1) {
    return cmax(0.0, cmin(a1, b1) - cmax(a0, b0));
}

// Exact form. Scale is a rational, so if cell and pixel bounds are carried in
// units of 1/den the overlap is an integer and coverage is the exact rational
// overlap/den. Going through `double p = num/den` instead loses this: 21/5 is
// not representable, 15 * 4.2 is not 63, and a sprite that should tile the
// device grid exactly acquires a phantom column. The rational is in the type;
// the arithmetic should use it.
constexpr long ispan(long a0, long a1, long b0, long b1) {
    long lo = a0 > b0 ? a0 : b0, hi = a1 < b1 ? a1 : b1;
    return hi > lo ? hi - lo : 0;
}
constexpr int idiv_floor(long a, long b) { long q = a / b; return int((a % b && (a < 0) != (b < 0)) ? q - 1 : q); }
constexpr int idiv_ceil (long a, long b) { return -idiv_floor(-a, b); }

constexpr double cpow(double b, int e) { double r = 1; for (int i = 0; i < e; ++i) r *= b; return r; }
// enough for the sRGB transfer curve; Newton on x^(1/2.4) and its inverse
constexpr double cpow_f(double b, double e) {
    if (b <= 0) return 0;
    // ln via atanh series, exp via Taylor — both converge fast on [0,1]
    double z = (b - 1) / (b + 1), z2 = z * z, s = 0, t = z;
    for (int k = 0; k < 40; ++k) { s += t / (2 * k + 1); t *= z2; }
    double ln = 2 * s, x = e * ln, r = 1, term = 1;
    for (int k = 1; k < 40; ++k) { term *= x / k; r += term; }
    return r;
}

// ───────────────────────────────────────────────────── scale as a rational
// The seam lattice period q is the reduced denominator (aamazing.md §5). Making
// it a compile-time constant is what lets the symmetry theorems be checkable.

template <int Num, int Den>
struct Scale {
    static_assert(Den > 0 && Num > 0, "scale must be positive");
    static constexpr int    num   = Num / cgcd(Num, Den);
    static constexpr int    den   = Den / cgcd(Num, Den);
    static constexpr int    q     = den;                    // seam lattice period
    static constexpr double value = double(Num) / double(Den);
};

template <int N> using Int = Scale<N, 1>;

// ───────────────────────────────────────────────────── composite operators
// Each carries its closed form (ng §6) and, crucially, whether it commutes —
// which is what decides whether traversal order is a light direction and
// whether a mirror survives fractional scale.

struct Px { double r = 0, g = 0, b = 0, a = 0; };            // premultiplied, working space

struct Over {
    static constexpr std::string_view name = "over";
    static constexpr bool commutative = false;
    static constexpr Px apply(Px d, double sr, double sg, double sb, double a) {
        return { sr * a + d.r * (1 - a), sg * a + d.g * (1 - a),
                 sb * a + d.b * (1 - a), a + d.a * (1 - a) };
    }
};
struct Add {
    static constexpr std::string_view name = "add";
    static constexpr bool commutative = true;
    static constexpr Px apply(Px d, double sr, double sg, double sb, double a) {
        return { d.r + sr * a, d.g + sg * a, d.b + sb * a, cmin(1.0, d.a + a) };
    }
};
struct Mul {
    static constexpr std::string_view name = "mul";
    static constexpr bool commutative = false;
    // the /255 is the normalisation that keeps a product of two 8-bit
    // quantities an 8-bit quantity; omitting it returns ~40000 (ng §6)
    static constexpr Px apply(Px d, double sr, double sg, double sb, double a) {
        return { d.r * (1 - a) + (d.r * sr / 255.0) * a,
                 d.g * (1 - a) + (d.g * sg / 255.0) * a,
                 d.b * (1 - a) + (d.b * sb / 255.0) * a, a + d.a * (1 - a) };
    }
};

// ───────────────────────────────────────────────────────── colour spaces
struct sRGB {
    static constexpr std::string_view name = "srgb";
    static constexpr bool linear = false;
    static constexpr double encode(double c) { return c; }
    static constexpr double decode(double c) { return c; }
};
struct Linear {
    static constexpr std::string_view name = "linear";
    static constexpr bool linear = true;
    static constexpr double encode(double c) {                  // sRGB → light
        double u = c / 255.0;
        return (u <= 0.04045 ? u / 12.92 : cpow_f((u + 0.055) / 1.055, 2.4)) * 255.0;
    }
    static constexpr double decode(double c) {                  // light → sRGB
        double u = c / 255.0;
        return (u <= 0.0031308 ? u * 12.92 : 1.055 * cpow_f(u, 1.0 / 2.4) - 0.055) * 255.0;
    }
};

// ─────────────────────────────────────────────────────── traversal orders
// A total order on cells. Only meaningful under a non-commutative operator —
// see `order_is_observable` below, which the compiler can tell you about.

struct RowMajor {
    static constexpr std::string_view name = "rowMajor";
    template <int N>
    static constexpr long key(const std::array<int, N>& c, const std::array<int, N>& dim) {
        long k = 0; for (int i = N - 1; i >= 0; --i) k = k * dim[i] + c[i]; return k;
    }
};
struct ColMajor {
    static constexpr std::string_view name = "colMajor";
    template <int N>
    static constexpr long key(const std::array<int, N>& c, const std::array<int, N>& dim) {
        long k = 0; for (int i = 0; i < N; ++i) k = k * dim[i] + c[i]; return k;
    }
};
struct Serpentine {
    static constexpr std::string_view name = "serpentine";
    template <int N>
    static constexpr long key(const std::array<int, N>& c, const std::array<int, N>& dim) {
        long row = 0; for (int i = N - 1; i >= 1; --i) row = row * dim[i] + c[i];
        long x = (row % 2) ? (dim[0] - 1 - c[0]) : c[0];
        return row * dim[0] + x;
    }
};
struct Radial {
    static constexpr std::string_view name = "radial";
    template <int N>
    static constexpr long key(const std::array<int, N>& c, const std::array<int, N>& dim) {
        long d2 = 0;
        for (int i = 0; i < N; ++i) { long t = 2 * c[i] - (dim[i] - 1); d2 += t * t; }
        return d2;
    }
};
struct Morton {
    static constexpr std::string_view name = "morton";
    template <int N>
    static constexpr long key(const std::array<int, N>& c, const std::array<int, N>&) {
        long k = 0;
        for (int bit = 0; bit < 16; ++bit)
            for (int i = 0; i < N; ++i)
                k |= long((c[i] >> bit) & 1) << (bit * N + i);
        return k;
    }
};
// paint in order of field value: the apparent light is generated by the same
// arithmetic as the form (ng §5). Needs Φ's function, not Φ's output — which is
// why the field is a template parameter of the sprite, not a runtime lookup.
struct ByField { static constexpr std::string_view name = "byField"; };

// ──────────────────────────────────────────────────────────────── fields
// n-generic by construction: every field is defined on a coordinate array, so
// the same set serves 2-D sprites and 3-D voxel models (ng §9).

struct Circle  { static constexpr std::string_view name = "circle";
    template <int N> static constexpr long f(const std::array<long, N>& u) {
        long s = 0; for (int i = 0; i < N; ++i) s += u[i] * u[i]; return s; } };
struct Diamond { static constexpr std::string_view name = "diamond";
    template <int N> static constexpr long f(const std::array<long, N>& u) {
        long s = 0; for (int i = 0; i < N; ++i) s += (u[i] < 0 ? -u[i] : u[i]); return s; } };
struct Cheby   { static constexpr std::string_view name = "cheby";
    template <int N> static constexpr long f(const std::array<long, N>& u) {
        long m = 0; for (int i = 0; i < N; ++i) { long a = u[i] < 0 ? -u[i] : u[i]; if (a > m) m = a; } return m; } };
struct Xor     { static constexpr std::string_view name = "xor";
    template <int N> static constexpr long f(const std::array<long, N>& u) {
        long s = 0; for (int i = 0; i < N; ++i) s ^= (u[i] < 0 ? -u[i] : u[i]); return s; } };
struct Product { static constexpr std::string_view name = "product";
    template <int N> static constexpr long f(const std::array<long, N>& u) {
        long s = 1; for (int i = 0; i < N; ++i) s *= u[i]; return s; } };
struct Trefoil { static constexpr std::string_view name = "trefoil";   // x³ − 3xy² lifted
    template <int N> static constexpr long f(const std::array<long, N>& u) {
        long s = 0; for (int i = 0; i < N; ++i) { long t = u[(i + 1) % N];
            s += u[i] * u[i] * u[i] - 3 * u[i] * t * t; } return s; } };
struct Hyper   { static constexpr std::string_view name = "hyper";     // alternating squares
    template <int N> static constexpr long f(const std::array<long, N>& u) {
        long s = 0; for (int i = 0; i < N; ++i) s += (i % 2 ? -1 : 1) * u[i] * u[i]; return s; } };

// ────────────────────────────────────────────── hyperoctahedral group B_n
// The signed permutation group: |B_n| = 2^n · n!. B_2 is the 8-element dihedral
// symmetry menu of the original tool; B_3 is the 48 elements of the voxel work;
// B_4 has 384. Generated at compile time so ng §9's "Γ generalises from the
// dihedral groups to the B_n family" is a fact about the code, not a plan.

template <int N> constexpr int factorial() { int r = 1; for (int i = 2; i <= N; ++i) r *= i; return r; }
template <int N> inline constexpr int B_order = (1 << N) * factorial<N>();

template <int N>
struct Element {                                   // x'_i = sign_i * x_{perm_i}
    std::array<int, N> perm{};
    std::array<int, N> sign{};
    constexpr std::array<long, N> operator()(const std::array<long, N>& x) const {
        std::array<long, N> y{};
        for (int i = 0; i < N; ++i) y[i] = sign[i] * x[perm[i]];
        return y;
    }
};

template <int N>
constexpr std::array<Element<N>, B_order<N>> make_B() {
    std::array<Element<N>, B_order<N>> out{};
    std::array<int, N> perm{};
    for (int i = 0; i < N; ++i) perm[i] = i;
    int idx = 0;
    // lexicographic permutations × all sign patterns
    for (int pcount = 0; pcount < factorial<N>(); ++pcount) {
        for (int mask = 0; mask < (1 << N); ++mask) {
            Element<N> e{};
            e.perm = perm;
            for (int i = 0; i < N; ++i) e.sign[i] = (mask >> i) & 1 ? -1 : 1;
            out[idx++] = e;
        }
        // next permutation, in place
        int i = N - 2; while (i >= 0 && perm[i] >= perm[i + 1]) --i;
        if (i < 0) break;
        int j = N - 1; while (perm[j] <= perm[i]) --j;
        { int t = perm[i]; perm[i] = perm[j]; perm[j] = t; }
        for (int a = i + 1, b = N - 1; a < b; ++a, --b) { int t = perm[a]; perm[a] = perm[b]; perm[b] = t; }
    }
    return out;
}
template <int N> inline constexpr auto B = make_B<N>();

// A symmetry is a subgroup, given as a predicate over B_n elements. The sprite
// is a section of the quotient: each cell takes the value of the lexicographic
// least point in its orbit, so the grid is symmetric by construction rather
// than by a fold pass that can be got wrong.
struct Asymmetric {                                    // trivial subgroup: no symmetry
    template <int N> static constexpr bool member(const Element<N>& e) {
        for (int i = 0; i < N; ++i) { if (e.perm[i] != i || e.sign[i] != 1) return false; }
        return true; } };
struct MirrorX {                                       // reflect axis 0 only
    template <int N> static constexpr bool member(const Element<N>& e) {
        for (int i = 0; i < N; ++i) { if (e.perm[i] != i) return false; }
        for (int i = 1; i < N; ++i) { if (e.sign[i] != 1) return false; }
        return true; } };
struct MirrorAll {                                     // reflect every axis: 2^N elements
    template <int N> static constexpr bool member(const Element<N>& e) {
        for (int i = 0; i < N; ++i) { if (e.perm[i] != i) return false; }
        return true; } };
struct Dihedral   { template <int N> static constexpr bool member(const Element<N>&) { return true; } };

// ────────────────────────────────────────── theorems, as compile-time predicates
//
// aamazing.md §9 proves the *phase field* is mirror-symmetric iff q ≤ 2 and
// q | (w−1). ng §0 corrects the conclusion drawn from it: that is a statement
// about the seam lattice, not about pixels. Under a non-commutative operator
// every mirrored pair of seams is resolved in the opposite sense, so the
// rendered image is symmetric only when there are no seams at all.

template <class S, int W> inline constexpr bool phase_field_symmetric_v =
    (S::q <= 2) && ((W - 1) % S::q == 0);

template <class Op, class S, int W> inline constexpr bool mirror_survives_v =
    Op::commutative ? ((long(W) * S::num) % S::den == 0)   // q | W
                    : (S::q == 1);                          // no seams at all

// Traversal order is observable only under a non-commutative operator (ng §5).
template <class Op, class Ord> inline constexpr bool order_is_observable_v =
    !Op::commutative && !std::is_same_v<Ord, RowMajor>;

// Fraction of device pixels that lie on the seam sublattice, per axis.
template <class S> inline constexpr double seam_share_v =
    S::q == 1 ? 0.0 : double(S::q - 1) / (double(S::q) * S::value);

// Opt-in guard. Instantiate it and a configuration that claims a mirror it
// does not have becomes a compile error rather than a surprising PNG.
template <class Op, class S, int W>
struct require_mirror {
    static_assert(mirror_survives_v<Op, S, W>,
        "This (operator, scale, width) does not render a mirror. Under a "
        "non-commutative operator the mirror survives only at integer scale; "
        "under a commutative one it survives iff q divides the width. "
        "See mosprites-ng.md §0 and §6.");
};

// ──────────────────────────────────────────────────────── sprite generator Φ
// mulberry32, verbatim from the reference implementation, so a seed means the
// same thing here as it does in the browser.

struct Rng {
    std::uint32_t s;
    constexpr explicit Rng(std::uint32_t seed) : s(seed) {}
    constexpr double operator()() {
        s += 0x6D2B79F5u;
        std::uint32_t t = s;
        t = std::uint32_t((t ^ (t >> 15)) * (t | 1u));
        t ^= t + std::uint32_t((t ^ (t >> 7)) * (t | 61u));
        return double((t ^ (t >> 14)) >> 8) / double(1u << 24);
    }
};

struct RGB { std::uint8_t r = 0, g = 0, b = 0; };

// bit-plane palette: bpc bits per channel, so colours land on a regular lattice
constexpr RGB palette_entry(std::uint32_t seed, int i, int bpc) {
    Rng rng(seed + std::uint32_t(i) * 2654435761u);
    int levels = (1 << bpc) - 1;
    auto ch = [&] { return std::uint8_t(double(int(rng() * (levels + 1)) % (levels + 1)) / levels * 255.0); };
    return { ch(), ch(), ch() };
}

template <int N> struct Extent { std::array<int, N> n{}; 
    constexpr int total() const { int t = 1; for (int i = 0; i < N; ++i) t *= n[i]; return t; } };

// Φ's configuration is a type. Note `Field` is carried in the type, which is
// how order:ByField gets the *function* rather than a reconstruction of it
// (ng §2.2 — a reconstructed function is a wrong function with no error).
template <class FieldT, class SymT, int Modulus, int Ncol = 4, int Bpc = 2,
          int Stride = 1, int CoveragePct = 50>
struct SpriteCfg {
    using Field = FieldT;
    using Sym   = SymT;
    static constexpr int modulus = Modulus, ncol = Ncol, bpc = Bpc,
                         stride = Stride, coverage_pct = CoveragePct;
    static_assert(Modulus > 1, "modulus must exceed 1");
    static_assert(Ncol >= 1 && Ncol <= 16, "palette size out of range");
};

template <class Cfg, int N, int... Dims>
struct Sprite {
    static constexpr int    rank = N;
    static constexpr Extent<N> dim{ { Dims... } };
    static constexpr int    cells = (Dims * ...);
    static_assert(sizeof...(Dims) == N, "one extent per axis");

    std::array<std::uint8_t, cells> grid{};        // 0 = empty, else palette index+1
    std::array<RGB, Cfg::ncol>      colors{};

    // centred, doubled coordinates: integral, so the group acts exactly
    static constexpr std::array<long, N> centred(const std::array<int, N>& c) {
        std::array<long, N> u{};
        for (int i = 0; i < N; ++i) u[i] = 2L * c[i] - (dim.n[i] - 1);
        return u;
    }
    // the field, published as a function of cell coordinates (ng §2.2)
    static constexpr long field_at(const std::array<int, N>& c) {
        auto u = centred(c);
        for (int i = 0; i < N; ++i) u[i] *= Cfg::stride;
        return Cfg::Field::template f<N>(u);
    }
    static constexpr int index(const std::array<int, N>& c) {
        int k = 0; for (int i = N - 1; i >= 0; --i) k = k * dim.n[i] + c[i]; return k;
    }
    static constexpr std::array<int, N> coord(int k) {
        std::array<int, N> c{};
        for (int i = 0; i < N; ++i) { c[i] = k % dim.n[i]; k /= dim.n[i]; }
        return c;
    }
    constexpr std::uint8_t at(const std::array<int, N>& c) const { return grid[index(c)]; }
};

template <class Cfg, int N, int... Dims>
constexpr Sprite<Cfg, N, Dims...> generate(std::uint32_t seed) {
    using Sp = Sprite<Cfg, N, Dims...>;
    Sp s{};
    Rng rng(seed);
    for (int i = 0; i < Cfg::ncol; ++i) s.colors[i] = palette_entry(seed ^ 0x9E3779B9u, i, Cfg::bpc);
    // luminance sort keeps palette index monotone in brightness
    for (int i = 1; i < Cfg::ncol; ++i)
        for (int j = i; j > 0; --j) {
            auto L = [](RGB c) { return 0.2126 * c.r + 0.7152 * c.g + 0.0722 * c.b; };
            if (L(s.colors[j]) < L(s.colors[j - 1])) { auto t = s.colors[j]; s.colors[j] = s.colors[j - 1]; s.colors[j - 1] = t; }
        }
    const long phase_off = long(rng() * 17) - 8;
    for (int k = 0; k < Sp::cells; ++k) {
        auto c = Sp::coord(k);
        auto u = Sp::centred(c);
        // section of the quotient: every orbit resolves to one representative,
        // so the grid carries the symmetry exactly, with no fold pass
        std::array<long, N> best = u;
        for (const auto& e : B<N>) {
            if (!Cfg::Sym::template member<N>(e)) continue;
            auto v = e(u);
            for (int i = 0; i < N; ++i) { if (v[i] < best[i]) { best = v; break; } if (v[i] > best[i]) break; }
        }
        for (int i = 0; i < N; ++i) best[i] = best[i] * Cfg::stride + phase_off;
        long n = Cfg::Field::template f<N>(best);
        long vm = ((n % Cfg::modulus) + Cfg::modulus) % Cfg::modulus;
        double t = double(vm) / double(Cfg::modulus);
        if (t >= double(Cfg::coverage_pct) / 100.0) { s.grid[k] = 0; continue; }
        int band = int(t / (double(Cfg::coverage_pct) / 100.0) * Cfg::ncol);
        s.grid[k] = std::uint8_t(1 + (band < Cfg::ncol ? band : Cfg::ncol - 1));
    }
    return s;
}

// ───────────────────────────────────────────────── the rasteriser, rank-N
// Coverage is separable: a = product over axes of the 1-D overlap. Walk<N,Axis>
// unrolls that product into a compile-time loop nest, which is the entirety of
// what ng §9 means by "the n-dimensional generalisation replaces the two nested
// loops with a recursive walk accumulating a = prod f_i. Nothing else changes."

// Splat<N,Axis> composites one cell into the buffer. The recursion is over
// axes, not pixels: at each level it multiplies in that axis's 1-D overlap and
// descends, so the terminal case receives a = prod f_i already formed. In two
// dimensions this is the familiar doubly-nested loop; in three it is the same
// function with one more instantiation. That equivalence is ng §9's argument.
//
// The buffer is threaded as a pointer rather than captured in a lambda: GCC 13
// ICEs on constexpr evaluation of a mutating by-reference capture over an array
// member, and depending on a compiler bug is not a design.
template <int N, int Axis>
struct Splat {
    // lo and num are in units of 1/den, so every overlap below is an integer
    template <class Op, class Buf>
    static constexpr void go(const std::array<long, N>& lo, long num, long den,
                             const std::array<int, N>& dev, std::array<int, N>& idx,
                             double acc, Buf& buf, double sr, double sg, double sb) {
        const int a = imax(0, idiv_floor(lo[Axis], den));
        const int b = imin(dev[Axis], idiv_ceil(lo[Axis] + num, den));
        for (int i = a; i < b; ++i) {
            const long ov = ispan(lo[Axis], lo[Axis] + num, long(i) * den, long(i + 1) * den);
            if (ov <= 0) continue;
            idx[Axis] = i;
            Splat<N, Axis + 1>::template go<Op>(lo, num, den, dev, idx,
                                                acc * (double(ov) / double(den)), buf, sr, sg, sb);
        }
    }
};
template <int N>
struct Splat<N, N> {
    template <class Op, class Buf>
    static constexpr void go(const std::array<long, N>&, long, long,
                             const std::array<int, N>& dev, std::array<int, N>& idx,
                             double acc, Buf& buf, double sr, double sg, double sb) {
        std::size_t k = 0;
        for (int i = N - 1; i >= 0; --i) k = k * std::size_t(dev[i]) + std::size_t(idx[i]);
        buf[k] = Op::apply(buf[k], sr, sg, sb, acc);
    }
};

// Phase applies to where a sprite is PLACED, never to the cells inside it.
// Snapping each cell would round away the fractional offsets that are the seam
// field, which is the entire effect — the sprite would render as if at integer
// scale. Local snaps the sprite's origin so it rasterises identically wherever
// it sits; Global leaves it, and the texture flows across a sheet. ng §3.
enum class Phase { Global, Local };

template <class ScaleT, class OpT, class OrderT, class SpaceT,
          Phase Ph = Phase::Local, int MatteR = -1, int MatteG = 0, int MatteB = 0>
struct Policy {
    using S = ScaleT; using Op = OpT; using Order = OrderT; using Space = SpaceT;
    static constexpr Phase  phase  = Ph;
    static constexpr bool   matted = MatteR >= 0;
    static constexpr RGB    matte  { std::uint8_t(MatteR < 0 ? 0 : MatteR),
                                     std::uint8_t(MatteG), std::uint8_t(MatteB) };
    static constexpr double p = ScaleT::value;
    static constexpr int    q = ScaleT::q;
};

template <int N> struct Raster {
    std::array<int, N>  dev{};
    std::vector<Px>     px;                        // premultiplied, working space
    constexpr int total() const { int t = 1; for (int i = 0; i < N; ++i) t *= dev[i]; return t; }
};

// The output contract (ng §4): the kernel accumulates premultiplied colour in
// the working space and converts nothing. resolve() is the single boundary
// step, and it unpremultiplies *before* applying the transfer curve because
// both operations are defined only on unpremultiplied colour.
template <class Space>
constexpr std::array<std::uint8_t, 4> resolve_px(Px v) {
    if (v.a <= 1e-9) return { 0, 0, 0, 0 };
    double r = v.r / v.a, g = v.g / v.a, b = v.b / v.a;          // 1. unpremultiply
    if constexpr (Space::linear) {                                // 2. transfer curve
        r = Space::decode(r); g = Space::decode(g); b = Space::decode(b);
    }
    auto q8 = [](double c) { return std::uint8_t(c < 0 ? 0 : (c > 255 ? 255 : int(c + 0.5))); };
    return { q8(r), q8(g), q8(b), q8(v.a * 255.0) };
}

template <class Pol, class Sp>
Raster<Sp::rank> rasterize(const Sp& s, std::array<double, Sp::rank> origin = {}) {
    constexpr int N = Sp::rank;
    using Op = typename Pol::Op; using Space = typename Pol::Space;

    constexpr long num = Pol::S::num, den = Pol::S::den;
    std::array<long, N> org{};                       // in units of 1/den
    for (int i = 0; i < N; ++i) org[i] = long(origin[i] * den + (origin[i] < 0 ? -0.5 : 0.5));
    if constexpr (Pol::phase == Phase::Local)
        for (int i = 0; i < N; ++i) org[i] = long(idiv_floor(org[i] + den / 2, den)) * den;

    Raster<N> out;
    for (int i = 0; i < N; ++i) out.dev[i] = idiv_ceil(org[i] + long(Sp::dim.n[i]) * num, den);
    out.px.assign(std::size_t(out.total()), Px{});

    if constexpr (Pol::matted) {
        Px m{ Space::encode(Pol::matte.r), Space::encode(Pol::matte.g),
              Space::encode(Pol::matte.b), 1.0 };
        for (auto& v : out.px) v = m;
    }

    // traversal order: a permutation of cell indices, chosen at compile time
    std::vector<int> seq(std::size_t(Sp::cells));
    for (int i = 0; i < Sp::cells; ++i) seq[std::size_t(i)] = i;
    {
        std::vector<long> key(std::size_t(Sp::cells));
        for (int i = 0; i < Sp::cells; ++i) {
            auto c = Sp::coord(i);
            if constexpr (std::is_same_v<typename Pol::Order, ByField>)
                key[std::size_t(i)] = Sp::field_at(c);           // Φ's function, not a copy of it
            else
                key[std::size_t(i)] = Pol::Order::template key<N>(c, Sp::dim.n);
        }
        // stable insertion sort: order must be deterministic, and ties must not
        // depend on the sort's internals (the CA work showed scan-order ties
        // silently destroy equivariance)
        for (std::size_t i = 1; i < seq.size(); ++i)
            for (std::size_t j = i; j > 0 && key[std::size_t(seq[j])] < key[std::size_t(seq[j - 1])]; --j) {
                int t = seq[j]; seq[j] = seq[j - 1]; seq[j - 1] = t;
            }
    }

    for (int ci : seq) {
        const std::uint8_t v = s.grid[std::size_t(ci)];
        if (!v) continue;
        const RGB col = s.colors[std::size_t(v - 1)];
        const double sr = Space::encode(col.r), sg = Space::encode(col.g), sb = Space::encode(col.b);

        auto c = Sp::coord(ci);
        std::array<long, N> lo{};
        for (int i = 0; i < N; ++i) lo[i] = org[i] + long(c[i]) * num;
        std::array<int, N> idx{};
        Splat<N, 0>::template go<Op>(lo, num, den, out.dev, idx, 1.0, out.px, sr, sg, sb);
    }
    return out;
}

// ───────────────────────────────────────────── constexpr path: one pixel
// The closed forms in aamazing.md §3 and §6 are claims about a single device
// pixel, so the compile-time entry point is a single device pixel. sample()
// walks the cells in traversal order and composites those whose box overlaps
// the requested pixel — the same arithmetic as rasterize(), without allocating
// a buffer, which keeps it usable in a constant expression.
//
// (It is O(cells) per pixel. That is the right trade for an oracle: the test
// suite checks a handful of pixels, and the whole-image path is rasterize().)

template <class Pol, class Sp>
constexpr Px sample(const Sp& s, const std::array<int, Sp::rank>& q,
                    std::array<double, Sp::rank> origin = {}) {
    constexpr int N = Sp::rank;
    using Op = typename Pol::Op;
    using Space = typename Pol::Space;

    constexpr long num = Pol::S::num, den = Pol::S::den;
    std::array<long, N> org{};
    for (int i = 0; i < N; ++i) org[i] = long(origin[i] * den + (origin[i] < 0 ? -0.5 : 0.5));
    if constexpr (Pol::phase == Phase::Local)
        for (int i = 0; i < N; ++i) org[i] = long(idiv_floor(org[i] + den / 2, den)) * den;

    Px acc{};
    if constexpr (Pol::matted)
        acc = Px{ Space::encode(Pol::matte.r), Space::encode(Pol::matte.g),
                  Space::encode(Pol::matte.b), 1.0 };

    // traversal order: visit cells by ascending key, ties by index. Selection
    // rather than a sort, so no scratch storage is needed.
    long last_key = 0; int last_idx = -1;
    for (int drawn = 0; drawn < Sp::cells; ++drawn) {
        long best_key = 0; int best = -1;
        for (int ci = 0; ci < Sp::cells; ++ci) {
            auto c = Sp::coord(ci);
            long k = 0;
            if constexpr (std::is_same_v<typename Pol::Order, ByField>) k = Sp::field_at(c);
            else k = Pol::Order::template key<N>(c, Sp::dim.n);
            if (last_idx >= 0 && (k < last_key || (k == last_key && ci <= last_idx))) continue;
            if (best < 0 || k < best_key) { best_key = k; best = ci; }
        }
        if (best < 0) break;
        last_key = best_key; last_idx = best;

        const std::uint8_t v = s.grid[std::size_t(best)];
        if (!v) continue;
        const RGB col = s.colors[std::size_t(v - 1)];
        auto c = Sp::coord(best);

        double a = 1.0;
        for (int i = 0; i < N && a > 0; ++i) {
            const long lo = org[i] + long(c[i]) * num;
            a *= double(ispan(lo, lo + num, long(q[i]) * den, long(q[i] + 1) * den)) / double(den);
        }
        if (a <= 0) continue;
        acc = Op::apply(acc, Space::encode(col.r), Space::encode(col.g), Space::encode(col.b), a);
    }
    return acc;
}

} // namespace mosprite