# Unified CRT Torus Braid DAG: Graded Sidon Energy with Directional Hierarchical Pruning

## Provenance (Clean Room)

This document synthesizes mathematical patterns from two external works, as
recorded in `CITATION.cff` references [0] and [1]. All implementation is
original to SilverSight — the external works inspired the *structural
analogies* and *design patterns*, not the code or theorems.

| External work | Pattern borrowed | Our adaptation |
|---------------|-----------------|----------------|
| ppf-contact-solver cubic barrier | Curvature linear in gap: ψ''(g) = 4(1-g/ĝ) | SidonEnergy(gap) = 4(1-gap/M) |
| ppf-contact-solver elasticity-inclusive stiffness | Coupling constraint stiffness into material stiffness | Coupled axis-swap × adjustment DAG traversal |
| ppf-contact-solver eigen-filtering | max(λ, 0) SPD projection | DAG branch pruning by Sidon-energy sign |
| ppf-contact-solver two-pass allocation | Dry pass (sparsity discovery) + fill-in (value computation) | DAG topology BFS + Sidon value fill |
| NAADF AADF 6-direction encoding | 6 directional distances instead of 1 scalar SDF | Chiral (identity, reflection) pairing per strand |
| NAADF 3-level nested hierarchy | Voxel → Block → Chunk | Crossing → Strand → DAG |
| NAADF max safe step formula | step_d = (1+bound-offset)/\|rayDir\| | Capacity = spacing/2 per direction |
| NAADF iterative distance transform | 3-iteration AADF propagation | Sidon potential bound propagation |
| NAADF hash-deduplicated blocks | Content-addressable 64-voxel store | Content-addressable strand/χ-state store |

---

## 1. Graded Sidon Energy (Cubic Barrier Analog)

### Current problem

The DAG currently uses `is_sidon(A)` as a binary predicate. A set is either
Sidon or it isn't. This gives no gradient signal when traversing — crossings
either succeed (reach Sidon) or fail (don't), with no intermediate information
about *which crossings bring us closer*.

### Reformulation

Define a **Sidon energy** that is graded by the gap distance relative to M:

```
Let gap(A) = min_{a < a' < a'' < a'''} |(a' + a'') - (a''' + a)|   (collision gap)

If gap(A) > M:
    SidonEnergy = 0                     (Sidon achieved — no residual energy)
Else:
    ℰ = 4 * (1 - gap(A) / M)           (graded residual, 0 < ℰ ≤ 4)
```

The cubic barrier ψ(g) from ppf-contact-solver uses curvature = 4(1 - g/ĝ),
which is *linear in the gap*. Our SidonEnergy uses the same form, where:
- g → gap(A): how close the closest collision is to the wrapping modulus M
- ĝ → M: the wrapping modulus (the threshold at which collisions are guaranteed
  to wrap differently)

Properties:
- When gap(A) = 0 (collision at zero difference): ℰ = 4 (maximum energy)
- When gap(A) = M/2: ℰ = 2 (half energy)
- When gap(A) > M (Sidon): ℰ = 0 (converged)
- The derivative dℰ/dgap = -4/M is constant: energy decreases linearly as the
  gap opens up

### DAG integration

Replace the binary `is_sidon` gate with:

```
def sidon_residual(A: List[int], M: int) -> float:
    """Return 0 if Sidon, else ℰ ∈ (0, 4]."""
    gap = min_collision_gap(A)
    return 0 if gap > M else 4.0 * (1.0 - gap / M)

def crossing_accepted(parent_energy: float, child_energy: float) -> bool:
    """A crossing is accepted iff it does not increase Sidon energy."""
    return child_energy <= parent_energy + EPSILON
```

This is the clean-room analog of eigen-filtering (max(λ, 0)): crossings that
*would increase the residual* are pruned, guaranteeing monotonic DAG descent.

---

## 2. Directional Chiral Decomposition (AADF Analog)

### Current model

Each strand has a chiral pair (L_id, L_ref). A crossing toggles between
over (+) and under (-), adjusting these two values. The adjustment is
isotropic: +2 on the active value, -1 on the passive value.

### Reformulation

The NAADF insight is that directional bounds are independent. A ray moving
+x doesn't care about the distance in +y. We apply the same independence
to crossing directions:

For strand i with pair (L_id, L_ref), define **4 directional capacities**:

```
cap⁺_id(i)  =  # of over-crossings possible before L_id exceeds band
cap⁻_id(i)  =  # of under-crossings possible before L_id drops below band
cap⁺_ref(i) =  # of over-crossings possible before L_ref exceeds band
cap⁻_ref(i) =  # of under-crossings possible before L_ref drops below band
```

Each directional capacity is the number of steps before the adjusted value
reaches the nearest other strand's band. This mirrors NAADF's 6-directional
AADF (2 bits per direction at block level, 5 bits at chunk level).

### Capacity encoding

Pack the 4 directional capacities into a single 8-bit word per strand:

```
Bits 0-1: cap⁺_id(i)  (2 bits, range 0-3 at strand level)
Bits 2-3: cap⁻_id(i)  (2 bits, range 0-3)
Bits 4-5: cap⁺_ref(i) (2 bits, range 0-3)
Bits 6-7: cap⁻_ref(i) (2 bits, range 0-3)
```

At the DAG level (aggregated across all strands), promote to 4 bits per
direction (range 0-15). The promotion rule mirrors NAADF's bound promotion:
`dag_cap = min(strand_cap × 4 + intra_strand_offset, 15)`.

---

## 3. Three-Level Nested DAG Hierarchy

### NAADF hierarchy (inspiration)

| Level | Grid | Voxels per | Bits per element | Purpose |
|-------|------|------------|-----------------|---------|
| Voxel | 4³ | 1 | 16 (2 bits × 6 dir + 1 occupancy) | Per-voxel state |
| Block | 4³ blocks | 64 | 32 (2 bits × 6 dir + 2 state) | Small-group aggregate |
| Chunk | N/16 chunks | 4096 | 32 (5 bits × 6 dir + 2 state) | Large-region empty skip |

### Our hierarchy

| Level | Contains | States per | Bits per element | Encoding | Analog to |
|-------|----------|-----------|-----------------|----------|-----------|
| Crossing | Single σ⁺/σ⁻ | 1 crossing | 8 (4 dir caps) | Per-crossing directional capacity | Voxel-level AADF |
| Strand | 2 crossings (id+ref) | ~15 crossings | 16 (4 dir caps × 4 bits) | Aggregated strand capacity | Block-level AADF |
| DAG | 2N moduli (N strands) | All crossings | 32 (4 dir caps × 8 bits) | Global Sidon potential | Chunk-level AADF |

### Hierarchy rules

**Propagation (bottom-up):** After each crossing on strand i, recompute the
strand-level capacities by scanning the current crossing-level capacities.
Then aggregate to DAG-level:

```
strand_cap[d] = min(crossing_cap[d] for crossing in strand)    # worst case
dag_cap[d]    = min(strand_cap[d] for strand in dag)            # global worst case
```

**Skip (top-down):** If the DAG-level capacity in a direction is 0, no strand
has remaining capacity in that direction — the entire DAG branch can be pruned.
This mirrors NAADF's chunk-level empty skip: if the chunk bound says 31 empty
voxels ahead, skip all 31 without descending.

```
def skip_direction(dag, direction: str) -> bool:
    """Return True if no strand can absorb another crossing in this direction."""
    return dag_cap[direction] == 0
```

---

## 4. Coupled Axis-Swap × Adjustment DAG Traversal

### Current problem

Axis-swap and adjustment are separate models. Axis-swap permutes reflection
moduli (FA-invariant, YB ✓), adjustment changes modulus values (FA-changing,
YB ✗). The DAG tries both independently, but they don't interact.

### Reformulation

The ppf-contact-solver's **elasticity-inclusive dynamic stiffness** couples
the constraint stiffness (contact) into the material stiffness (elasticity).
We do the same: couple the topology stiffness (axis-swap) into the resource
stiffness (adjustment).

**How it works:**

```
def coupled_crossing(pairs, s, direction):
    """Single coupled crossing: axis-swap then adjust."""

    # Step 1: Axis-swap topology (changes modulus ordering only)
    pairs = axis_swap(pairs, s)

    # Step 2: Adjustment with direction-dependent intensity
    # The adjustment step size is scaled by the topology permutation:
    #   - If strand s just received a new modulus via swap, the adjustment
    #     is larger (the topology "stiffens" the crossing)
    #   - If strand s kept its modulus (no swap effect), adjustment is
    #     the standard ±2/∓1
    stiffness = topology_stiffness(pairs, s)   # ∈ [1, 2]
    pairs = adjust(pairs, s, direction * stiffness)

    return pairs if pairwise_coprime(pairs) else None
```

Where `topology_stiffness` is 2 if the swap changed the reflection modulus
ordering, 1 otherwise. This is the clean-room analog of the elasticity-
inclusive stiffness term (barrier.cu:48-85).

### Three-phase traversal

```
Phase 1 (Graded Sidon search):
  Use coupled crossings with small bands (gap ~ 60, capacity ~ 15 per strand).
  Track SidonEnergy residual. Prune branches that increase ℰ.

Phase 2 (DAG topology expansion):
  Once ℰ < threshold, expand to wide bands (gap ~ 500, capacity ~ 125 per strand).
  The DAG-level capacity-0 check prunes exhausted regions.

Phase 3 (Content-addressable dedup):
  Deduplicate identical modulus configurations via hash map (NAADF analog of
  chunkCalc.fx:57-115). If two DAG nodes have identical moduli and FA values,
  they share the same strand-state record.
```

---

## 5. Algorithm: Unified DAG Build

```python
def build_unified_dag(A0, S, n_strands, max_steps):
    """Build the CRT torus DAG with graded Sidon energy + hierarchical pruning."""

    # Initialize: chiral pairs with per-directional capacity
    pairs = chiral_pairs(n_strands, max_crossings=15)
    caps = compute_directional_capacities(pairs)

    # Root node: initial Sidon energy
    root = DAGNode(pairs, A0, caps)
    root.energy = sidon_energy(A0, modulus_product(pairs))

    queue = [root]
    visited = {}

    while queue:
        node = queue.pop(0)

        # Skip if DAG-level capacity is 0 in all directions
        for d in ['id_over', 'id_under', 'ref_over', 'ref_under']:
            if dag_capacity(node, d) <= 0:
                continue  # can't cross in this direction

        # Coupled crossing on each strand
        for s in range(n_strands):
            for direction in ['over', 'under']:
                child = coupled_crossing(node, s, direction)
                if child is None:
                    continue  # coprimality failed

                # Compute Sidon energy and prune if it increases
                child.energy = sidon_energy(child.A, child.M)
                if child.energy > node.energy + EPSILON:
                    continue  # monotonicity violated: prune

                # Dedup against visited states
                h = hash(child.moduli)
                if h in visited and visited[h].A == child.A:
                    continue  # content-addressable dedup

                visited[h] = child
                queue.append(child)
```

---

## 6. Implementation Plan

| Component | File | Status | Notes |
|-----------|------|--------|-------|
| `sidon_energy` | `scripts/sidon_energy.py` | New | Graded energy from gap/M using cubic form |
| `directional_capacity` | `scripts/full_chiral_dag.py` | Extend | Add 4-directional capacity tracking |
| `coupled_crossing` | `scripts/full_chiral_dag.py` | Extend | Axis-swap then adjust with topology stiffness |
| `hierarchical_skip` | `scripts/full_chiral_dag.py` | Extend | DAG-level capacity-0 pruning |
| `content_addressable_store` | `scripts/full_chiral_dag.py` | New | Hash-dedup modulus configs |
| `unified_dag_build` | `scripts/run_8strand_search.py` | Extend | Replace BFS with unified algorithm |
| `SidonEnergy` theorem | `formal/CoreFormalism/SidonEnergy.lean` | New | Lemma: ℰ monotonic under accepted crossings |
| `CoupledCrossing` lemma | `formal/CoreFormalism/CoupledCrossing.lean` | New | YB preservation under coupled model |

---

## 7. SidonEnergy Gradient via Asymmetric Scoring Identity

### Mixedbread's scoring identity

The key mathematical insight from mixedbread's asymmetric quantization
(CITATION.cff [2], blog 2026-06-29):

```
q · b = 2 * Σ_{b_i=+1} q_i - Σ q_i
```

A binary × int8 dot product needs only:
1. Precompute Σ q_i (once per query)
2. Sum query dimensions where document bit = +1 (the "selected" sum)
3. Apply the identity: 2×selected − total

No multiplication per dimension needed. Just a conditional add and a shift.

### Our analog: SidonEnergy gradient

For a crossing on strand i, define a **crossing sign** s_i ∈ {+1, −1} encoding
over/under, and a **crossing contribution** c_i (the change in collision gap
attributable to strand i). The SidonEnergy before and after a crossing relates
as:

```
ℰ_after = ℰ_before − (2 * s_i * c_i) / M
```

Derivation:
- ℰ = max(0, 4(1 − gap/M)) for non-Sidon states
- Δgap = s_i * c_i (the gap change from strand i's crossing, signed)
- Δℰ = 4/M * (−Δgap) = −4 * s_i * c_i / M
- So ℰ_after = ℰ_before − (2 * 2 * s_i * c_i / M)

The factor of **2** appears for the same reason as in mixedbread's identity:
the active crossing direction contributes with double weight (+2 step) while
the passive direction contributes with single weight (−1 step). This is an
**asymmetric scoring kernel** embedded in the CRT arithmetic.

### Practical benefit

Replace:
```python
child_energy = sidon_energy(child.A, child.M)   # full recompute
if child_energy > node.energy: continue          # O(N²) collision check
```

With:
```python
# Gradient update: O(1) per crossing
delta = -4 * crossing_sign * crossing_contribution / total_modulus
child_energy = node.energy + delta
if child_energy > node.energy: continue          # same monotonicity gate
```

This is the clean-room analog of the NEON SDOT kernel: precompute the
"query sum" (ℰ_before, once per DAG node), then for each outgoing crossing,
compute only the "selected" part (the contribution of the crossed strand)
and apply the identity.

### Asymmetric precision in the DAG

Mixedbread's crucial storage insight applies directly:

| Side | Mixedbread | Our DAG |
|------|-----------|---------|
| Query / Topology | int8 (high precision, short-lived) | Braid word (1 bit per crossing, cheap to store) |
| Document / Sidon | binary (low precision, dominates storage) | SidonEnergy and FA values (full precision, recomputed on-demand) |

Just as mixedbread stores document vectors as 1-bit signs and keeps query at
int8, we store the **braid word** (which crossings happened) as a bit field
(1 bit per crossing type), while the **FA values and SidonEnergy** are
recomputed from scratch on each node access.

A complete 8-strand crossing history fits in **1 byte** (1 bit per strand for
over/under, or 2 bits per strand with directional encoding). This mirrors
mixedbread's 32× storage reduction: the braid word is the "binary document"
that dominates storage cost, while the DAG traversal is the "int8 query" that
dominates compute cost.

### Scoring kernel for the CRT residue update

Mixedbread's kernel avoids full multiply via precomputed query planes. Our
kernel for CRT residue update:

```python
def crossing_residue(residue_before: int, step: int, mod: int) -> int:
    """CRT residue update: no multiply needed.

    Analogous to mixedbread's q·b identity avoiding full dot product.
    """
    return (residue_before + step) % mod   # single add + modulo
```

A full CRT reconstruction of FA after N crossings would be O(N × num_moduli).
With the gradient identity, each crossing update is O(1): just update the
residue for the crossed strand and compute Δℰ from the signed contribution.

---

## 8. Matrix Orthogonalization of the Modulus Configuration

### Newton-Schulz for the CRT modulus matrix

The mLSTM maintains a memory matrix C ∈ ℝ^{d×d}. Each read is a matrix-vector
product. The Newton-Schulz iteration enforces orthogonality:

```
M ← (3M − M·M^T·M) / 2      (5 iterations → M^T·M ≈ I)
```

This prevents mode collapse: a few strong directions dominating the memory,
crowding out weaker memories. The +15-45% NAR accuracy gain comes from this
equalization (CITATION.cff [3], blog 2026-06-30).

**Our analog:** The CRT modulus configuration is an n×2 matrix (n strands,
2 moduli each). The "mode collapse" analog is coprimality exhaustion: a few
strands' moduli grow large while others stay small, eventually hitting the
Q16_16 bound while other strands have unused capacity.

Define the **modulus orthogonality** constraint:

```
For all i ≠ j:  gcd(L_id_i, L_id_j) = 1
                 gcd(L_id_i, L_ref_j) = 1
                 gcd(L_ref_i, L_ref_j) = 1
```

This is already satisfied by construction (individual primes with spacing).
But the **distribution** of moduli can become unbalanced after many crossings.
The Newton-Schulz analog is a **redistribution step** that normalizes the
modulus set:

```
1. Compute Frobenius norm of the modulus matrix:  ‖M‖_F = √(Σ L_i²)
2. Normalize:  L_i ← L_i / ‖M‖_F × target_norm
3. Re-discretize to nearest integer coprime with all other moduli
```

This is not a literal NS iteration (our moduli are discrete, not continuous),
but the **intent** is the same: prevent a few directions from dominating.

### Read-only orthogonalization

The critical design choice from the mLSTM experiment: **orthogonalize during
reads, don't write back**. Writing back the orthogonalized memory degraded
performance because it destroyed the information stored in the memory state.

**Our mapping:**

```
Read path (axis-swap):     orthogonalize → apply YB constraint
                              → axis-swap is FA-invariant (CRT symmetry)
                              → YB ensures the braid word is consistent
                              → safe to orthogonalize: no information loss

Write path (adjustment):   don't orthogonalize
                              → adjustment changes FA values
                              → orthogonalizing after adjustment would
                                destroy the Sidon state we just created
                              → read-only orthogonalization preserves the
                                Sidon state while keeping the topology clean
```

This is exactly the read-only pattern from the mLSTM experiment, and it
validates our dual-model decomposition: axis-swap (topology, orthogonalized)
and adjustment (resource, unconstrained).

### Capacity equalization

Muon's optimizer orthogonalizes momenta to prevent strong directions from
dominating. The result is that weaker directions get lifted. In our model,
this maps to **capacity redistribution**:

```
After each k crossings:
  1. Compute cap_remaining per direction per strand
  2. If max(cap) / min(cap) > 4:   # imbalance threshold
       Redistribute: strand with min cap gets a modulus reset
       (new modulus further from neighbors)
  3. Graft the new modulus into the existing DAG node
     (verify coprimality first)
```

This prevents the "mode collapse" where one strand exhausts its capacity
while others have slack. The threshold of 4 is arbitrary — like NS iteration
count, it needs empirical calibration.

---

## 9. Open Questions

1. **Topology stiffness function** — should the coupling factor be binary (1 or
   2) or continuous? The ppf-contact-solver uses a continuous stiffness term
   (the projected elastic Hessian), but our crossing values are discrete
   (integers). A continuous stiffness would round to integer, which may lose
   the coupling benefit.

2. **Directional capacity vs actual coprimality** — capacity is a heuristic
   (band spacing / step size). Actual coprimality can fail earlier if adjusted
   values happen to coincide with another strand's modulus. The capacity bound
   is safe (never overestimates) but may be conservative.

3. **Three-phase transition thresholds** — when does Phase 1 end and Phase 2
   begin? The SidonEnergy threshold (ℰ < 0.1?) needs empirical calibration.

4. **Hash dedup collision rate** — NAADF uses open-addressing with linear
   probing. Our modulus state space is smaller (16 ints per node), so a simple
   Python dict may suffice. Formal verification will need a hash lemma.
