Freeze SmoothLife model contract and fixtures
This commit is contained in:
627
docs/model-contract.md
Normal file
627
docs/model-contract.md
Normal file
@@ -0,0 +1,627 @@
|
||||
# SmoothLife model contract
|
||||
|
||||
This document is the implementation authority for the mathematical model. It is intentionally independent of the retired application and is sufficient to implement CPU and GPU backends. `docs/historical-options.md` narrows the three historical choices; no other compatibility behavior exists.
|
||||
|
||||
## 1. Normative conventions
|
||||
|
||||
- Unless a section explicitly says `f32`, the semantic oracle and golden values use IEEE-754 binary64 (`f64`), round-to-nearest, ties-to-even. A production backend may use `f32`, but it must compare to the semantic fixtures under the tolerances below; reassociation or contraction is permitted only when the resulting values remain within the applicable tolerance.
|
||||
- `TAU = 2π`. The old shader identifier named `pi` represented `TAU`, not π.
|
||||
- `clamp01(x) = min(1, max(0, x))`. Configuration and input state must be finite; NaN and infinity are validation errors rather than values to be ordered by `min`/`max`.
|
||||
- `mix(x,y,q) = (1-q)x + qy` in that written operation order.
|
||||
- `wrap(i,N) = i mod N` with a result in `[0,N)`, including for negative `i`.
|
||||
- A sum is accumulated in the canonical enumeration order stated below. CPU reference sums use `f64` without parallel reassociation.
|
||||
- Equality signs and strict inequalities in this document are normative.
|
||||
|
||||
Contract fixtures encode authoritative f64 inputs/outputs as JSON numbers using enough decimal digits for binary64 round-trip. Branch/equality classification, indices, clamps, and values produced only by basic operations must be bit-exact in the scalar f64 CPU oracle. Values involving `sin`, `cos`, `atan`, `acos`, `exp`, or long reductions use `|actual-reference| ≤ atol+rtol*|reference|`: `atol=rtol=2^-48` for scalar curve/geometry fixtures and `2^-42` for normalized fields, integration fields, and sphere area/constant-field checks.
|
||||
|
||||
For a production f32 CPU/GPU backend compared directly with the f64 semantic fixtures, use `atol=rtol=2^-20` for scalar curves/rules and kernel weights, `2^-16` for one-step normalized planar or multiscale fields, and `2^-15` for one-step sphere or delayed-time fields. These are acceptance limits, not permission to alter branch selection, indexing, support membership, clamp points, or history-layer choice. NaN or infinity always fails. The retained packed-FFT stage-specific f32 tolerances are in section 6.2.
|
||||
|
||||
The common state invariant is:
|
||||
|
||||
> State is a scalar `A ∈ [0,1]` after every committed update.
|
||||
|
||||
Intermediate rule targets and derivatives need not be in `[0,1]`. They are clamped only at the explicit points below.
|
||||
|
||||
## 2. Domains, storage, and indexing
|
||||
|
||||
Base domains are a periodic one-dimensional circle, two-dimensional torus, or three-dimensional torus. Every active extent is positive. Storage is x-fast row-major:
|
||||
|
||||
```text
|
||||
index1(x) = x
|
||||
index2(x,y) = x + Nx*y
|
||||
index3(x,y,z) = x + Nx*(y + Ny*z)
|
||||
```
|
||||
|
||||
Coordinates are wrapped before indexing. Canonical field enumeration is `z`, then `y`, then `x`, with `x` innermost; absent dimensions are omitted. The sample count is `V=Nx`, `Nx*Ny`, or `Nx*Ny*Nz`.
|
||||
|
||||
For an even extent `N`, the signed periodic offset represented by an index `i∈[0,N)` is exactly
|
||||
|
||||
```text
|
||||
offset(i,N) = i when i < N/2
|
||||
i - N otherwise
|
||||
```
|
||||
|
||||
Thus the Nyquist sample `i=N/2` has offset `-N/2`. For odd `N`, use `i` when `i≤floor(N/2)` and `i-N` otherwise. Euclidean lattice distance is the square root of the sum of squared signed offsets.
|
||||
|
||||
A sphere is one closed spherical surface represented by six active face arrays, not six independent boundaries. Delayed time has a two-dimensional periodic spatial torus and a circular temporal history.
|
||||
|
||||
## 3. Sampled kernels and convolution
|
||||
|
||||
### 3.1 Kernel functions
|
||||
|
||||
For the base model, validate `ra>0`, `rr>0`, and `rb>0`, then define
|
||||
|
||||
```text
|
||||
ri = ra / rr
|
||||
w = ra / rb
|
||||
```
|
||||
|
||||
For `w>0`, the sampled ramp is
|
||||
|
||||
```text
|
||||
L(r;a,w) = 0 when r < a-w/2
|
||||
1 when r > a+w/2
|
||||
(r-a)/w + 1/2 otherwise
|
||||
```
|
||||
|
||||
At `r=a-w/2` the third branch evaluates to exactly zero; at `r=a+w/2` it evaluates to exactly one. The kernels are
|
||||
|
||||
```text
|
||||
KD(r) = 1 - L(r;ri,w)
|
||||
KR(r) = L(r;ri,w) * (1 - L(r;ra,w))
|
||||
```
|
||||
|
||||
Weights are sampled at lattice points. A support entry exists iff its evaluated weight is strictly greater than zero. Generate the complete nonzero support by examining every unique periodic offset in the domain; never stop at the historical connected-component cutoff. Entries are stored in canonical field-index order.
|
||||
|
||||
For a radial kernel `K`, circular convolution is
|
||||
|
||||
```text
|
||||
(A * K)(x) = Σ_o A(wrap(x-o)) K(o)
|
||||
```
|
||||
|
||||
where wrapping is component-wise. Neighborhood fields are
|
||||
|
||||
```text
|
||||
M = (A * KD) / Σ_o KD(o)
|
||||
N = (A * KR) / Σ_o KR(o)
|
||||
```
|
||||
|
||||
The two sums must be finite and strictly positive or validation fails. Division occurs after the complete numerator sum. A known two-dimensional check, using the canonical signed-offset grid and `ra=10`, `rr=3`, `rb=10`, is
|
||||
|
||||
```text
|
||||
sum(KR) ≈ 279.216312
|
||||
sum(KD) ≈ 35.524035
|
||||
```
|
||||
|
||||
The check is diagnostic rather than a decimal truncation target.
|
||||
|
||||
### 3.2 Kernel suitability diagnostics
|
||||
|
||||
A configuration remains deterministic when support wraps, but must emit a non-fatal `support_touches_nyquist` warning if any nonzero support sample has `|o_j|=N_j/2` on an even axis. Equivalently for an isotropic base ring, the warning begins when `ra+w/2` is strictly greater than half an active extent; equality has zero outer-boundary weight and does not warn by itself. Emit `support_covers_domain` if every unique periodic offset has a nonzero disk or ring weight. Zero or non-finite normalization is an error, not a warning. These diagnostics never truncate support or alter weights.
|
||||
|
||||
## 4. Rule curves
|
||||
|
||||
All smooth widths `e`, including `sn` and `sm`, must be finite and strictly positive. Let
|
||||
|
||||
```text
|
||||
u = (x-a+e/2)/e
|
||||
logistic(x;a,e) = 1 / (1 + exp(-4(x-a)/e))
|
||||
```
|
||||
|
||||
The rising curve `P_t(x,a,e)` is selected by `t∈0..7`:
|
||||
|
||||
```text
|
||||
0 hard: x >= a ? 1 : 0
|
||||
1 linear: 0 if x < a-e/2
|
||||
1 if x > a+e/2
|
||||
u otherwise
|
||||
2 Hermite: 0 if x < a-e/2
|
||||
1 if x > a+e/2
|
||||
u²(3-2u) otherwise
|
||||
3 sine: 0 if x < a-e/2
|
||||
1 if x > a+e/2
|
||||
0.5*sin((TAU/2)*(x-a)/e)+0.5 otherwise
|
||||
4 logistic: logistic(x;a,e)
|
||||
5 atan: atan((x-a)*(TAU/2)/e)/(TAU/2)+0.5
|
||||
6 atan/cos: (1.1*atan((x-a)/e)/(TAU/4)*cos(1.4*(x-a))+1)/2
|
||||
7 overshoot: (logistic(x;a,e)-0.5)
|
||||
*(1+exp(-(x-a)²/e²))+0.5
|
||||
```
|
||||
|
||||
At `x=a±e/2`, compact curves 1–3 use their formula branch (not the outer constant branch); that branch selection is fixed for fixtures. Curves 4–7 are not compactly clamped.
|
||||
|
||||
For window types `t∈0..7`,
|
||||
|
||||
```text
|
||||
W_t(n;a,b,e) = P_t(n,a,e) * (1-P_t(n,b,e))
|
||||
```
|
||||
|
||||
A hard window is exactly `[a,b)`: it is one at `n=a` and zero at `n=b`. Types 8 and 9 are complete windows and do not invoke `P_8` or `P_9`:
|
||||
|
||||
```text
|
||||
base = logistic(n;a,e) * (1-logistic(n;b,e))
|
||||
mid = (a+b)/2
|
||||
g = exp(-(20*(n-mid))²)
|
||||
W_8 = base*(1-0.2*g)
|
||||
W_9 = base*(1+0.2*g)
|
||||
```
|
||||
|
||||
The mixer curve has its own selector `m∈0..7` only:
|
||||
|
||||
```text
|
||||
Q = P_m(M, 0.5, sm)
|
||||
```
|
||||
|
||||
Curves 6 and 7, mixer values, interpolated thresholds, windows, and resulting targets can overshoot. Do not clamp them during rule evaluation.
|
||||
|
||||
With `B=W_t(N;b1,b2,sn)` and `D=W_t(N;d1,d2,sn)`, the four rule constructions are exactly
|
||||
|
||||
```text
|
||||
1: S = mix(B,D,M)
|
||||
2: S = mix(B,D,Q)
|
||||
3: S = W_t(N, mix(b1,d1,M), mix(b2,d2,M), sn)
|
||||
4: S = W_t(N, mix(b1,d1,Q), mix(b2,d2,Q), sn)
|
||||
```
|
||||
|
||||
`M` and `N` are sampled normalized neighborhoods. No implicit clamp occurs between any expression above.
|
||||
|
||||
## 5. Dynamics and integration
|
||||
|
||||
For state `Y`, let `S(Y)` mean recomputing its neighborhoods and target.
|
||||
|
||||
```text
|
||||
Discrete: next = clamp01(S(A))
|
||||
Growth: f(Y) = 2*S(Y)-1
|
||||
Relaxation: f(Y) = S(Y)-Y
|
||||
```
|
||||
|
||||
Discrete dynamics ignore `dt` and the selected integrator. Experimental numeric modes 3 and 4 are unsupported.
|
||||
|
||||
### 5.1 Euler and Adams–Bashforth
|
||||
|
||||
Euler commits
|
||||
|
||||
```text
|
||||
A_next = clamp01(A + dt*k_n), k_n=f(A)
|
||||
```
|
||||
|
||||
AB3 has deterministic startup and clamps only the committed result:
|
||||
|
||||
```text
|
||||
no prior derivative: Δ = k_n
|
||||
one prior derivative: Δ = (3*k_n-k_{n-1})/2
|
||||
at least two derivatives: Δ = (23*k_n-16*k_{n-1}+5*k_{n-2})/12
|
||||
A_next = clamp01(A + dt*Δ)
|
||||
```
|
||||
|
||||
After a successful commit, shift in the derivative evaluated at the pre-step state. Undefined or invalid history is never partially reused.
|
||||
|
||||
### 5.2 RK4 and relaxation reference
|
||||
|
||||
Let `A0` be the state at step start. Every intermediate state and the final state is clamped element-wise:
|
||||
|
||||
```text
|
||||
k1 = f1(A0)
|
||||
Y2 = clamp01(A0 + dt*k1/2)
|
||||
k2 = f2(Y2)
|
||||
Y3 = clamp01(A0 + dt*k2/2)
|
||||
k3 = f3(Y3)
|
||||
Y4 = clamp01(A0 + dt*k3)
|
||||
k4 = f4(Y4)
|
||||
A_next = clamp01(A0 + dt*(k1+2*k2+2*k3+k4)/6)
|
||||
```
|
||||
|
||||
For Growth, every `fj(Y)=2*S(Y)-1`. For Relaxation the local option is:
|
||||
|
||||
```text
|
||||
StageState (default): f1(A0)=S(A0)-A0; fj(Y)=S(Y)-Y
|
||||
StepOrigin: f1(A0)=S(A0)-A0; fj(Y)=S(Y)-A0, j=2..4
|
||||
```
|
||||
|
||||
`Rk4RelaxationReference` exists structurally only for Relaxation+RK4. It is absent—not stored and ignored—for Discrete, Growth, Euler, and AB3.
|
||||
|
||||
### 5.3 AB history invalidation
|
||||
|
||||
Preserve AB derivative history across presentation-only changes, inspection changes, pause/resume, and scheduler-rate changes. Invalidate it completely on:
|
||||
|
||||
- reset;
|
||||
- ordinary state load;
|
||||
- initializer or seed change/use;
|
||||
- rule change;
|
||||
- dynamics change;
|
||||
- timestep change;
|
||||
- integrator change;
|
||||
- kernel geometry change;
|
||||
- any historical-option change;
|
||||
- topology or shape change;
|
||||
- variant change;
|
||||
- backend change.
|
||||
|
||||
Switching away from AB3 and later back never revives derivatives. An exact continuation checkpoint may restore history only after its full validated model and backend-independent arithmetic descriptor matches: schema/version, variant, topology, shape, scalar representation, rule, dynamics, `dt`, integrator, kernels, applicable historical options, and arithmetic contract. A backend identifier may differ only where that backend is certified to implement the same arithmetic descriptor. Any missing, malformed, non-finite, wrong-shaped, or mismatched derivative discards all AB history. This list is authoritative.
|
||||
|
||||
## 6. FFT algorithms
|
||||
|
||||
Planar and multiscale convolution have one local algorithm selection.
|
||||
|
||||
### 6.1 Standard
|
||||
|
||||
For sample count `V`, Standard uses the conventional DFT
|
||||
|
||||
```text
|
||||
F[k] = Σ_x A[x] exp(-i*TAU*k·x/N)
|
||||
A[x] = (1/V) Σ_k F[k] exp(+i*TAU*k·x/N)
|
||||
```
|
||||
|
||||
applied separably. Multiply state and sampled-kernel spectra, inverse transform, then divide by the corresponding sampled `f64` kernel sum. Forward is unscaled and inverse is scaled by `1/V`.
|
||||
|
||||
### 6.2 LegacyPackedUnitary arithmetic
|
||||
|
||||
`LegacyPackedUnitary` is a complete GPU algorithm, not a compatibility mode. Every active extent must be a power of two, `Nx≥2`, and x must be even. One-dimensional, two-dimensional, and three-dimensional arrays use adjacent-real x packing and a half-width spectrum.
|
||||
|
||||
For `m=Nx/2`, pack
|
||||
|
||||
```text
|
||||
p[j] = complex(A[2j], A[2j+1]), j=0..m-1
|
||||
```
|
||||
|
||||
and compute a unitary radix-2 FFT `P=U_m(p)`. For `k=0..m`, with indices modulo `m`,
|
||||
|
||||
```text
|
||||
Pk = P[k mod m]
|
||||
Qk = conjugate(P[(m-k) mod m])
|
||||
E = (Pk+Qk)/2
|
||||
O = (Pk-Qk)/(2i)
|
||||
X[k] = (E + exp(-i*TAU*k/Nx)*O)/sqrt(2)
|
||||
```
|
||||
|
||||
This is the unitary real transform. Transform the full y axis, then the full z axis, for every stored x frequency; absent axes are skipped. Each is the same unitary complex transform. The stored spectrum shape is `(Nx/2+1,Ny,Nz)` with omitted dimensions understood.
|
||||
|
||||
Inverse y and z stages run in reverse axis order. Recover packed x data for `k=0..m-1` by
|
||||
|
||||
```text
|
||||
E = (X[k] + conjugate(X[m-k]))/sqrt(2)
|
||||
O = conjugate(exp(-i*TAU*k/Nx))
|
||||
*(X[k] - conjugate(X[m-k]))/sqrt(2)
|
||||
P[k] = E + i*O
|
||||
```
|
||||
|
||||
then apply the unitary inverse `U_m⁻¹` and unpack real and imaginary components to adjacent real samples.
|
||||
|
||||
A unitary inverse of a spectral product produces circular convolution divided by `sqrt(V)`. Therefore each spectral product is corrected, before inverse stages, by exactly
|
||||
|
||||
```text
|
||||
correction = sqrt(V) / kernel_sum
|
||||
H[k] = correction * (F_A[k] * F_K[k])
|
||||
```
|
||||
|
||||
where `kernel_sum` is the complete sampled-kernel sum. This yields the normalized neighborhood directly.
|
||||
|
||||
The historical plan is an output-indexed radix-2 plan. For an ordinary complex axis of length `n=2^B`, execute stages `e=1..B`. At output index `q`, let `ℓ=2^e` and `j=q mod ℓ`. Select inputs
|
||||
|
||||
```text
|
||||
j < ℓ/2: ia=q, ib=q+ℓ/2
|
||||
otherwise: ia=q-ℓ/2, ib=q
|
||||
```
|
||||
|
||||
At `e=1` only, replace both selected indices by their `B`-bit reversals. Later stages use them directly. The plan twiddle and output are
|
||||
|
||||
```text
|
||||
w = exp(sign*i*TAU*j/ℓ) # sign=-1 forward, +1 inverse
|
||||
out[q] = (input[ia] + w*input[ib]) / sqrt(2)
|
||||
```
|
||||
|
||||
For `j≥ℓ/2`, the twiddle already contains the minus sign of the upper butterfly; do not rewrite that output as a separately ordered subtraction. Packed x first runs this plan at length `m=Nx/2` through stages `1..log2(m)`.
|
||||
|
||||
The x real/complex conversion is a distinct tangle stage and fixes the operation order behind the mathematical formulas above. For each stored `k=0..m`, set `a=P[k mod m]`, `b=conjugate(P[(m-k) mod m])`, and
|
||||
|
||||
```text
|
||||
forward w = exp(-i*TAU*(k/Nx+1/4))
|
||||
X[k] = (a+b + (a-b)*w) * (0.5/sqrt(2))
|
||||
```
|
||||
|
||||
For inverse conversion only `k=0..m-1` is produced. Set `a=X[k]`, `b=conjugate(X[m-k])`, and
|
||||
|
||||
```text
|
||||
inverse w = exp(+i*TAU*(k/Nx+1/4))
|
||||
P[k] = (a+b + (a-b)*w) * (0.5*sqrt(2))
|
||||
```
|
||||
|
||||
Then run ordinary packed-x stages with inverse sign. Forward axis order is packed x stages, forward tangle, y, z. Inverse order is z, y, inverse tangle, packed x stages. Stages ping-pong; no pass reads and writes the same resource.
|
||||
|
||||
All arithmetic in this subsection after input conversion is operation-ordered IEEE binary32. Let `fl` mean one round-to-nearest-ties-even `f32` operation. Subnormals are preserved and fused multiply-add is forbidden. An ordinary plan output is exactly
|
||||
|
||||
```text
|
||||
p0=fl(wr*br); p1=fl(wi*bi)
|
||||
r = fl(fl(ar+p0)-p1)
|
||||
p2=fl(wr*bi); p3=fl(wi*br)
|
||||
i = fl(fl(ai+p2)+p3)
|
||||
out.r=fl(r*INV_SQRT_2_F32)
|
||||
out.i=fl(i*INV_SQRT_2_F32)
|
||||
```
|
||||
|
||||
A tangle first forms `dr=fl(ar-br)`, `di=fl(ai-bi)`, evaluates `(dr+i*di)*w` with the same four-product order, forms `sr=fl(ar+br)` and `si=fl(ai+bi)`, adds product to sum component-wise, then multiplies by its rounded tangle scale. Conjugating `b` is an exact sign-bit change before those operations.
|
||||
|
||||
`INV_SQRT_2_F32`, both tangle scales, twiddle components, and `correction` are each rounded once to `f32` before use; twiddles are generated from the corresponding binary64 `sin`/`cos`. For spectral multiplication, first compute `br=fl(F_K.r*correction)` and `bi=fl(F_K.i*correction)`, then multiply `F_A` by that scaled complex kernel using the same four-product sequence. Packing and unpacking are exact component copies. This operation order defines the CPU stage oracle.
|
||||
|
||||
Integer plans and bit-reversal indices must match exactly. Contract fixtures encode packed values, every bit-reversal and butterfly stage, real/complex conversion, spectra, correction, and normalized convolution as hexadecimal `f32` bits. The scalar CPU oracle must be bit-exact to them. GPU fixture comparison uses
|
||||
|
||||
```text
|
||||
|actual-reference| <= atol + rtol*|reference|
|
||||
```
|
||||
|
||||
with finite values required:
|
||||
|
||||
| Check | `atol` | `rtol` |
|
||||
| --- | ---: | ---: |
|
||||
| packing and one-butterfly microfixtures | `0` (bit-exact) | `0` |
|
||||
| each complete FFT stage and real/complex conversion | `2^-20` | `2^-20` |
|
||||
| final half spectrum | `2^-18` | `2^-18` |
|
||||
| normalized convolution | `2^-16` | `2^-16` |
|
||||
|
||||
The operation-ordered f32 exception is local to this algorithm. Rule evaluation, integration, and their semantic fixtures remain defined by the f64 oracle, while a production GPU may keep convolution results in f32 and validate the subsequent common pipeline under the f32 backend tolerances above; no readback or widening is required during ordinary execution. An explicit request on CPU or unsupported GPU hardware fails clearly and never substitutes Standard. Complete kernels and safe resources are used by both algorithms.
|
||||
|
||||
## 7. Multiscale variant
|
||||
|
||||
There are exactly three scales, indexed in order `0,1,2`. Each has its own sampled disk/ring kernels, rule, and timestep `dt_i`. It supports Growth and corrected Relaxation only; Discrete is rejected because no target-aggregation policy is defined.
|
||||
|
||||
Two neighborhood interpretations are crossed with three composition policies, giving exactly six combinations:
|
||||
|
||||
| Inputs | Sequential | Ordered clamped sum | Mean increment |
|
||||
| --- | --- | --- | --- |
|
||||
| Independent | yes | yes | yes |
|
||||
| Chained | yes | yes | yes |
|
||||
|
||||
At any evaluated state `Y`, Independent supplies
|
||||
|
||||
```text
|
||||
scale i: (N_i,M_i) = (ring_i(Y), disk_i(Y))
|
||||
```
|
||||
|
||||
Chained supplies
|
||||
|
||||
```text
|
||||
scale 0: (N_0,M_0) = (ring_0(Y), ring_1(Y))
|
||||
scale 1: (N_1,M_1) = (ring_1(Y), ring_2(Y))
|
||||
scale 2: (N_2,M_2) = (ring_2(Y), disk_2(Y))
|
||||
```
|
||||
|
||||
Tuple order is always `(N,M)`. A needed field is recomputed from the specified state even if another scale has the same nominal field.
|
||||
|
||||
For a scale target `S_i(Y)`, define
|
||||
|
||||
```text
|
||||
Growth response: g_i(Y;Aref) = 2*S_i(Y)-1
|
||||
Relaxation response: g_i(Y;Aref) = S_i(Y)-Aref
|
||||
increment: I_i = dt_i*g_i
|
||||
```
|
||||
|
||||
Composition is:
|
||||
|
||||
1. **Sequential stage relaxation.** Set `Y_0=A0`. For `i=0,1,2`, recompute every required field from `Y_i`; use `Aref=Y_i` for Relaxation; then `Y_{i+1}=clamp01(Y_i+dt_i*g_i(Y_i;Y_i))`. Commit `Y_3`.
|
||||
2. **Ordered clamped sum.** Compute all three targets from one unchanged snapshot `A0`; Relaxation uses `Aref=A0`. Set `Y_0=A0`, then in scale order `Y_{i+1}=clamp01(Y_i+dt_i*g_i(A0;A0))`. Commit `Y_3`. Clamping the accumulator does not change any already computed response.
|
||||
3. **Mean increment.** Compute all targets from snapshot `A0`, with Relaxation reference `A0`, then commit `clamp01(A0+(I_0+I_1+I_2)/3)` with one final clamp.
|
||||
|
||||
This is per-scale Euler composition; the base Euler/AB3/RK4 selector does not apply. Shared-snapshot methods never read an intermediate accumulator. Additive discrete responses, stale chained fields, feedback reads, and undefined relaxation sources are rejected.
|
||||
|
||||
## 8. Sphere variant
|
||||
|
||||
### 8.1 Geometry and area
|
||||
|
||||
`K` must be positive and even; the default is `K=128`. There are six active `K×K` face arrays, x-fast within each face, in this fixed order and frame:
|
||||
|
||||
| face | normal `n` | `u_axis` | `v_axis` |
|
||||
| --- | --- | --- | --- |
|
||||
| 0 `+X` | `(1,0,0)` | `(0,1,0)` | `(0,0,1)` |
|
||||
| 1 `+Y` | `(0,1,0)` | `(-1,0,0)` | `(0,0,1)` |
|
||||
| 2 `-X` | `(-1,0,0)` | `(0,-1,0)` | `(0,0,1)` |
|
||||
| 3 `-Y` | `(0,-1,0)` | `(1,0,0)` | `(0,0,1)` |
|
||||
| 4 `+Z` | `(0,0,1)` | `(1,0,0)` | `(0,1,0)` |
|
||||
| 5 `-Z` | `(0,0,-1)` | `(1,0,0)` | `(0,-1,0)` |
|
||||
|
||||
`u_axis×v_axis=normal`. Let `R=K/2`. A face sample `(x,y)` is at
|
||||
|
||||
```text
|
||||
u = 2*(x+0.5)/K - 1
|
||||
v = 2*(y+0.5)/K - 1
|
||||
d = normalize(n + tan(u*π/4)*u_axis + tan(v*π/4)*v_axis)
|
||||
```
|
||||
|
||||
Cell area is the spherical area of its four mapped corner rays. For unit rays `a,b,c`, define
|
||||
|
||||
```text
|
||||
tri(a,b,c) = 2*atan2(abs(a·(b×c)), 1+a·b+b·c+c·a)
|
||||
```
|
||||
|
||||
With corners in `(u,v)` order `d00,d10,d11,d01`,
|
||||
|
||||
```text
|
||||
cell_area = R²*(tri(d00,d10,d11)+tri(d00,d11,d01))
|
||||
```
|
||||
|
||||
Canonical global enumeration is face `0..5`, then y, then x. Area totals must sum to `4πR²` within f64 fixture tolerance.
|
||||
|
||||
Geodesic distance is
|
||||
|
||||
```text
|
||||
distance(a,b) = R*acos(clamp(a·b,-1,1))
|
||||
```
|
||||
|
||||
The inner radius is first defined in planar units and both radii are converted once:
|
||||
|
||||
```text
|
||||
ri_planar = ra_planar/3
|
||||
ri = R*acos(1-ri_planar²/(2R²))
|
||||
ra = R*acos(1-ra_planar²/(2R²))
|
||||
```
|
||||
|
||||
Thus `ri=ra/3` is the planar-radius relation; nonlinear geodesic conversion is applied separately to each radius. Validate `0<ra_planar≤2R`; do not silently clamp either conversion argument. Both kernel transition widths are exactly `1`:
|
||||
|
||||
```text
|
||||
KD(r)=1-L(r;ri,1)
|
||||
KR(r)=L(r;ri,1)*(1-L(r;ra,1))
|
||||
```
|
||||
|
||||
### 8.2 Corrected sphere
|
||||
|
||||
`SphereModel::Corrected` is storage-independent. For every center sample, enumerate all samples in all six active arrays, evaluate geodesic distance and every nonzero kernel weight, and accumulate in canonical global order:
|
||||
|
||||
```text
|
||||
M(c) = Σ_j area(j)*KD(distance(c,j))*A(j)
|
||||
/ Σ_j area(j)*KD(distance(c,j))
|
||||
N(c) = Σ_j area(j)*KR(distance(c,j))*A(j)
|
||||
/ Σ_j area(j)*KR(distance(c,j))
|
||||
```
|
||||
|
||||
The denominators are the actual per-center sampled sums, not analytic cap areas or a shared center-independent approximation. This guarantees a constant field maps to the same constant, including at face edges and cube corners. Complete geodesic support includes every positive weight on any face; atlas adjacency and storage gutters must not affect the result. Nonpositive/non-finite denominators are errors. Emit a warning if the outer nonzero support crosses a hemisphere (`ra+0.5>πR/2`), but still globally enumerate it.
|
||||
|
||||
### 8.3 Legacy sphere
|
||||
|
||||
The precise retained Legacy geometry, deterministic side-gutter masking, normalization, and constraints are specified in `docs/historical-options.md`. It uses the same face frames and active-cell areas, but a face-local stencil rather than corrected global enumeration. Both sphere models clamp dot products, initialize all storage, and use ping-pong commits.
|
||||
|
||||
Sphere update modes are
|
||||
|
||||
```text
|
||||
Direct: A_next = clamp01(S)
|
||||
Smooth: A_next = clamp01(A + 0.1*(2*S-1))
|
||||
```
|
||||
|
||||
`0.1` is fixed and is not a configurable timestep.
|
||||
|
||||
## 9. Delayed-time variant
|
||||
|
||||
The delayed-time field is a periodic `Nx×Ny` torus with depth `D=16`. Its sampled spatial kernels use `ri=ra/3` and transition width `1` for both disk and ring. Let `h=ceil(ra+0.5)` and enumerate the direct stencil in `dy=-h..h`, then `dx=-h..h` order. Evaluate both weights at every offset in that square; zero-weight offsets may be skipped by an optimized implementation but do not affect either sum. Spatial coordinates wrap only when sampling history, so distinct stencil offsets remain distinct even if a small torus maps them to the same cell. Emit a non-fatal geometry warning when `h` exceeds half either spatial extent.
|
||||
|
||||
History is `H[0..15]`. `head` is the next layer to overwrite and
|
||||
|
||||
```text
|
||||
latest = wrap(head-1,16)
|
||||
delay(o) = floor(distance(o)+0.5)
|
||||
layer(o) = wrap(latest-delay(o),16)
|
||||
```
|
||||
|
||||
The `+0.5` equality rounds upward: a distance exactly `q+0.5` selects delay `q+1`. Delay zero reads the latest committed state. The delayed neighborhoods are
|
||||
|
||||
```text
|
||||
M(x) = Σ_o KD(distance(o))*H[layer(o)][wrap(x-o)] / Σ_o KD(distance(o))
|
||||
N(x) = Σ_o KR(distance(o))*H[layer(o)][wrap(x-o)] / Σ_o KR(distance(o))
|
||||
```
|
||||
|
||||
Both sums use the direct `dy`, then `dx` stencil order above. Temporal layer choice does not change normalization.
|
||||
|
||||
A step reads the entire old history and computes from `A=H[latest]`:
|
||||
|
||||
```text
|
||||
Direct: next = clamp01(S)
|
||||
Smooth: next = clamp01(H[latest] + 0.1*(2*S-1))
|
||||
```
|
||||
|
||||
It then writes `next` to `H[head]` and only afterward sets `head=wrap(head+1,16)`. No pass reads the layer it is currently writing. Initialization replicates one field into every layer and sets `head=0`, so the initial latest layer is 15. The old non-causal head anomaly is not configurable.
|
||||
|
||||
## 10. Deterministic PRNG
|
||||
|
||||
All seeded initializers use ChaCha12 with the original 16-word layout:
|
||||
|
||||
```text
|
||||
state[0..4] = LE words of ASCII "expand 32-byte k"
|
||||
state[4..12] = eight little-endian key words
|
||||
state[12..14] = little-endian 64-bit block counter
|
||||
state[14..16] = little-endian 64-bit stream value
|
||||
```
|
||||
|
||||
The 32-byte key is the seed encoded as one little-endian `u64` followed by 24 zero bytes. Counter and stream both start at zero. A block applies six ChaCha double rounds (12 rounds total), adds the original state word-by-word with wrapping `u32` addition, emits words `0..15` sequentially in little-endian order, and increments only the 64-bit counter modulo `2^64`.
|
||||
|
||||
A quarter round on `(a,b,c,d)` is exactly
|
||||
|
||||
```text
|
||||
a += b; d ^= a; d = rotl(d,16)
|
||||
c += d; b ^= c; b = rotl(b,12)
|
||||
a += b; d ^= a; d = rotl(d, 8)
|
||||
c += d; b ^= c; b = rotl(b, 7)
|
||||
```
|
||||
|
||||
Each double round applies columns `(0,4,8,12)`, `(1,5,9,13)`, `(2,6,10,14)`, `(3,7,11,15)`, then diagonals `(0,5,10,15)`, `(1,6,11,12)`, `(2,7,8,13)`, `(3,4,9,14)`.
|
||||
|
||||
`next_u32` consumes sequential output words. `next_u64` consumes two words and places the first in the low half:
|
||||
|
||||
```text
|
||||
lo = next_u32(); hi = next_u32()
|
||||
value = u64(lo) | (u64(hi)<<32)
|
||||
```
|
||||
|
||||
A unit binary64 draw is
|
||||
|
||||
```text
|
||||
unit = f64(next_u64() >> 11) * 2^-53 # [0,1)
|
||||
```
|
||||
|
||||
For an unbiased integer in `[0,b)`, require `0<b≤2^64`, set `limit=2^64-(2^64 mod b)`, repeatedly draw `x=next_u64()` while `x≥limit`, then return `x mod b`. Rejected words are consumed. For `b=2^64`, every word is accepted directly. An integer range offsets this result; a real `[lo,hi)` draw is `lo+(hi-lo)*unit` in written order. There are no hidden warm-up draws, per-thread streams, or implementation-library distribution calls.
|
||||
|
||||
## 11. Initializers
|
||||
|
||||
All initializers are newly authored, topology-native algorithms. They do not reproduce, capture, import, or migrate legacy preset state. Painting loops use canonical storage enumeration and consume no random values.
|
||||
|
||||
### 11.1 Periodic planar splats
|
||||
|
||||
Clear the field to zero. For dimension `d`, calculate once
|
||||
|
||||
```text
|
||||
count = floor(product(axis_extent)
|
||||
/ product(min(2*ra,axis_extent))) + 1
|
||||
```
|
||||
|
||||
using f64 products in axis order x, y, z. For each splat, draw continuous center coordinates in axis order (`x`, then `y`, then `z` as applicable), each uniform in `[0,extent)`, then draw radius uniform in `[0.5*ra,ra)`. Thus each splat consumes `d+1` unit draws.
|
||||
|
||||
For every lattice sample, compute component-wise shortest periodic distance from its integer coordinate to the continuous center, then Euclidean distance. Set the sample to one iff `distance < radius`; equality is not painted. Splats overwrite with one and overlaps consume no draws. This applies to the circle, 2-D torus, and 3-D torus.
|
||||
|
||||
### 11.2 Geodesic sphere overlays
|
||||
|
||||
Clear all six active face arrays to zero. Perform exactly 1,000 overlays. For each overlay, in this order:
|
||||
|
||||
1. `face = bounded_u64(6)`;
|
||||
2. `cx = bounded_u64(K)` and `cy = bounded_u64(K)`;
|
||||
3. `radius = 2 + bounded_u64(6)`, giving integers 2 through 7 inclusive;
|
||||
4. `value = bounded_u64(2)` converted to `0.0` or `1.0`.
|
||||
|
||||
Use the selected active sample `(face,cx,cy)` as the center. Enumerate every sample on every face and assign `value` iff its clamped-dot geodesic distance from the center is strictly `< radius`. Equality is not painted. Later overlays win. Global geodesic painting deliberately crosses every face seam and cube corner and is used for both sphere models; no atlas gutter is initialized by randomness.
|
||||
|
||||
### 11.3 Periodic delayed-time boxes
|
||||
|
||||
Require positive extents. Clear one `Nx×Ny` field to zero and perform exactly 1,000 boxes. For each box draw, in order,
|
||||
|
||||
```text
|
||||
x0 = bounded_u64(Nx)
|
||||
y0 = bounded_u64(Ny)
|
||||
w = 10 + bounded_u64(10) # 10..19 inclusive
|
||||
h = 10 + bounded_u64(10) # 10..19 inclusive
|
||||
```
|
||||
|
||||
Set to one all wrapped coordinates `(x0+dx,y0+dy)` with `dx∈[0,w)` and `dy∈[0,h)`. Boxes are half-open on their high edges and later boxes also write one. Replicate the completed field bit-for-bit into all 16 history layers and set `head=0`. There is no gradual fill and no undefined layer.
|
||||
|
||||
## 12. Newly authored preset policy
|
||||
|
||||
Bundled presets are authored directly in the versioned product schema. There is no legacy row capture, importer, conversion, deduplication, migration, count target, grouping target, or value-correspondence requirement. Each preset has a stable ID, name, description, tags, authoring provenance, explicit variant/rules/dynamics/shape/backend, explicit seed and initializer, and presentation recommendations. It stores only historical-option fields applicable to that run.
|
||||
|
||||
Unless a newly authored preset explicitly overrides them, authoring baselines are:
|
||||
|
||||
- Base: shape `1024`, `512×512`, or `64×64×64` by dimension; Euler; Standard FFT; seed 1; periodic planar splats; palette 2; 3-D volume style 2. Euler stores no RK4-reference field. A newly applicable Relaxation+RK4 field defaults to StageState.
|
||||
- Multiscale: `512×512`; Sequential; Independent; Growth or Relaxation only; Standard FFT; seed 1; periodic planar splats; palette 7.
|
||||
- Sphere: `K=128`, `R=64`; Corrected; seed 1; geodesic overlays; palette 1. A Legacy preset must be separately named and explicitly select only that local option.
|
||||
- Delayed time: `512×512×16` history; causal delay; seed 1; one periodic box field replicated to every layer; palette 7.
|
||||
|
||||
These are modern product-authoring choices, not inferred migration defaults. Overrides are explicit and remain subject to this contract.
|
||||
|
||||
## 13. Clamp-point summary
|
||||
|
||||
No rule curve, mixer, threshold interpolation, kernel convolution, or target is implicitly clamped. Clamp exactly at:
|
||||
|
||||
- Discrete direct commit;
|
||||
- Euler and every AB committed state;
|
||||
- every RK4 intermediate state and final state;
|
||||
- each Sequential multiscale stage;
|
||||
- each Ordered-clamped-sum accumulator stage;
|
||||
- the one Mean-increment multiscale commit;
|
||||
- sphere Direct/Smooth commits;
|
||||
- delayed-time Direct/Smooth commits;
|
||||
- dot products immediately before every `acos`, to `[-1,1]`.
|
||||
|
||||
Initializers produce only zero or one. This list is exhaustive.
|
||||
Reference in New Issue
Block a user