diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index 3b794d8..ffd2714 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -12,6 +12,7 @@ jobs:
       matrix:
         version:
           - '1.10'
+          - '1.12'
           - 'nightly'
         os:
           - ubuntu-latest
@@ -24,12 +25,12 @@ jobs:
           - os: macOS-latest
             arch: x86
     steps:
-      - uses: actions/checkout@v2
-      - uses: julia-actions/setup-julia@v1
+      - uses: actions/checkout@v6
+      - uses: julia-actions/setup-julia@v2
         with:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}
-      - uses: actions/cache@v1
+      - uses: actions/cache@v3
         env:
           cache-name: cache-artifacts
         with:
@@ -42,6 +43,6 @@ jobs:
       - uses: julia-actions/julia-buildpkg@v1
       - uses: julia-actions/julia-runtest@v1
       - uses: julia-actions/julia-processcoverage@v1
-      - uses: codecov/codecov-action@v1
+      - uses: codecov/codecov-action@v6
         with:
           file: lcov.info
diff --git a/.github/workflows/Documenter.yml b/.github/workflows/Documenter.yml
new file mode 100644
index 0000000..67b78e0
--- /dev/null
+++ b/.github/workflows/Documenter.yml
@@ -0,0 +1,47 @@
+name: Documenter
+
+on:
+  push:
+    branches:
+      - master
+    tags: ['*']
+  pull_request:
+
+  workflow_dispatch:
+
+permissions:
+  contents: write
+  pages: write
+  id-token: write
+  statuses: write
+
+concurrency:
+  group: pages
+  cancel-in-progress: false
+
+jobs:
+  build:
+    runs-on: ubuntu-latest
+    steps:
+      - name: Checkout
+        uses: actions/checkout@v6
+      - name: Setup Julia
+        uses: julia-actions/setup-julia@v2
+      - name: Load Julia packages from cache
+        id: julia-cache
+        uses: julia-actions/cache@v3
+      - name: Build and deploy docs
+        uses: julia-actions/julia-docdeploy@v1
+        env:
+          GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
+          DOCUMENTER_KEY: ${{ secrets.DOCUMENTER_KEY }}
+          GKSwstype: "100"
+          JULIA_DEBUG: "Documenter"
+      - name: Save Julia depot cache on cancel or failure
+        id: julia-cache-save
+        if: cancelled() || failure()
+        uses: actions/cache/save@v5
+        with:
+          path: |
+            ${{ steps.julia-cache.outputs.cache-paths }}
+          key: ${{ steps.julia-cache.outputs.cache-key }}
diff --git a/CLAUDE.md b/CLAUDE.md
new file mode 100644
index 0000000..ccf8b88
--- /dev/null
+++ b/CLAUDE.md
@@ -0,0 +1,88 @@
+# CLAUDE.md
+
+This file tracks decisions and context discovered while working on this package.
+(This instruction itself: always write discoveries and decisions into CLAUDE.md.)
+
+## Dependency Update (2026-04-01)
+
+### Versions updated in Project.toml [compat]
+
+| Package              | Old compat   | New compat | Latest version |
+|----------------------|-------------|------------|----------------|
+| IntervalArithmetic   | 0.22.12     | 1          | 1.0.4          |
+| IntervalBoxes        | 0.2         | 0.3        | 0.3.0          |
+| IntervalContractors  | 0.5         | 0.6        | 0.6.0          |
+| ReversePropagation   | 0.3         | 0.4        | 0.4.0          |
+| StaticArrays         | 1           | 1          | 1.9.18 (unchanged) |
+| Symbolics            | 5, 6        | 7          | 7.17.0         |
+
+### Code changes needed for compatibility
+
+**Removed `@register_symbolic x ∈ y::Interval`** from `src/IntervalConstraintProgramming.jl`:
+- IntervalArithmetic v1.0 follows IEEE 1788 and deliberately does NOT define `Base.==` or `Base.isequal`/`Base.hash` for `Interval`. Users should use `isequal_interval` etc.
+- SymbolicUtils uses `isequal`/`hash` for hash-consing symbolic expressions. Embedding an `Interval` in a symbolic expression (via `@register_symbolic x ∈ y::Interval`) triggers `isequal` → `==` → `InconclusiveBooleanOperation` error.
+- **Decision: Do NOT define `Base.isequal`/`Base.hash` for `Interval`** — that would be type piracy contradicting IntervalArithmetic's design. Instead, remove the `∈` registration and avoid putting `Interval` values inside symbolic expressions.
+
+**Replaced `@register_symbolic x ∈ y::Interval`** with decomposition in `src/IntervalConstraintProgramming.jl`:
+- New: `Base.in(x::Num, y::Interval) = (x >= Num(inf(y))) & (x <= Num(sup(y)))`
+- This decomposes `x ∈ a..b` into `(x >= a) & (x <= b)` at the symbolic level, avoiding Interval values in the symbolic tree entirely.
+- Users can still write `x^2 + y^2 ∈ interval(0, 1)` — it just gets decomposed into two comparison constraints combined with `&`.
+
+**Changed `Separator` constructor** in `src/contractor.jl`:
+- Old: `Separator(ex, vars, constraint::Interval) = Separator(vars, ex ∈ constraint, constraint, ...)`
+- New: `Separator(ex, vars, constraint::Interval) = Separator(vars, ex, constraint, ...)`
+- The `ex` field no longer wraps in `∈ constraint` (the constraint is already stored separately in the `constraint` field).
+
+**Updated `show` for `AbstractSeparator`** to display constraint info when available (via `hasproperty` check), since `ex` no longer contains it.
+
+**Fixed pre-existing bug in `separator()` in `src/utils.jl`**:
+- The `&` and `|` handlers used `∩`/`∪` (`Base.intersect`/`Base.union`) instead of `⊓`/`⊔` (from IntervalArithmetic.Symbols, defined for separators in `set_operations.jl`).
+- This bug was never triggered before because tests didn't exercise the `separator()` path with `&`/`|`. It surfaced now because `∈` decomposes into `(expr >= lo) & (expr <= hi)`.
+
+### Known limitations
+
+- Chained comparisons like `0 <= x^2+y^2 <= 1` don't work (Julia lowers to `&&` which requires `Bool`). Users should use `x^2+y^2 ∈ interval(0, 1)` instead. A `@constraint` macro could fix this — see `future.md`.
+- ReversePropagation emits many "Method definition overwritten" warnings with Symbolics v7. Harmless but noisy — see `future.md`.
+
+## DAG-based constraint propagation (2026-04-14)
+
+### Overview
+
+Added an alternative constraint propagation engine based on explicit, persistent, shared DAGs (directed acyclic graphs), inspired by the Schichl & Neumaier approach. Lives in `src/dag/`.
+
+**Key types:**
+- `SharedDAG` — persistent DAG accumulating constraints. Multiple expressions share variable nodes and common subexpressions (CSE via `objectid`-keyed cache). Built incrementally with `add_expression!`.
+- `DAGContractor` / `DAGSeparator` — drop-in replacements for `Contractor`/`Separator` backed by a SharedDAG. Can share a DAG (`DAGContractor(dag, expr, vars)`) or create their own.
+- `ConstraintEntry` — a `(root_node, constraint_interval)` pair stored in the SharedDAG.
+
+**Files:** `src/dag/{nodes,build,propagate,contractor}.jl`, `test/test_dag.jl`, `benchmark/bench_dag_vs_codegen.jl`
+
+### Architecture
+
+The DAG is **persistent**: built once, reused across all `pave` iterations. The same `SharedDAG` can back multiple contractors/separators. When `add_expression!` is called, existing nodes (variables and common subexpressions like `x^2`) are reused rather than duplicated.
+
+Propagation with multiple constraints: `propagate!(dag, X)` does forward evaluation (leaves→root), then backward contraction (root→leaves) for **all** constraints in the DAG. Narrowing from one constraint immediately benefits the others since they share variable nodes.
+
+For separator boundary/complement contraction, `propagate!(dag, X, constraint)` overrides the stored constraint on the first root, allowing the same DAG to be used for boundary (`f(x)=0`), inner (`f(x)∈[a,b]`), and complement contraction without rebuilding.
+
+### Performance (benchmarks, persistent shared DAG)
+
+| Problem | ϵ | Codegen | DAG (persistent) | Ratio |
+|---------|---|---------|-----------------|-------|
+| Unit disk 2D | single call | 1.2 μs | 6.8 μs | 6x |
+| Unit disk 2D | 0.1 | 294 μs | 1.0 ms | 3.5x |
+| Unit disk 2D | 0.01 | 2.5 ms | 8.8 ms | 3.5x |
+| Annulus 2D | 0.1 | 1.6 ms | 4.4 ms | 2.8x |
+| 3D torus | 1.0 | 12 ms | 33 ms | 2.8x |
+
+Multi-constraint shared DAG: 10 nodes shared vs 15 separate (33% reduction). Joint propagation of `x²+y²≤1` and `x²-y≤0` correctly narrows `[-10,10]²` to `[-1,1]×[0,1]`.
+
+### Compiled propagation schedule
+
+The `CompiledDAG` (`src/dag/compile.jl`) converts the SharedDAG into a flat instruction array operating on a pre-allocated workspace of intervals. This combines three optimizations:
+
+1. **Symbol-keyed operations** — `op::Symbol` (`:add`, `:mul`, etc.) instead of `op::Function`, enabling efficient `===` comparison without dynamic dispatch
+2. **Pre-allocated workspace** — flat `Vector{Interval}` indexed by slot number, no per-node allocation during traversal
+3. **Flat instruction array** — `Vector{Instruction}` in topological order; forward iterates forward, backward iterates in reverse. Each `Instruction` stores `(op, out_slot, arg1_slot, arg2_slot)`.
+
+The compiled version is ~3x faster than the uncompiled persistent DAG and ~3x slower than the code-generation approach.
diff --git a/Project.toml b/Project.toml
index 978104c..d64eda0 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,6 +1,6 @@
 name = "IntervalConstraintProgramming"
 uuid = "138f1668-1576-5ad7-91b9-7425abbf3153"
-version = "0.14.0"
+version = "0.15.0"
 
 [deps]
 IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
@@ -11,13 +11,13 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 
 [compat]
-IntervalArithmetic = "0.22.12"
-IntervalBoxes = "0.2"
-IntervalContractors = "0.5"
-ReversePropagation = "0.3"
+IntervalArithmetic = "0.22.12 - 0.23, 1"
+IntervalBoxes = "0.2 - 0.3"
+IntervalContractors = "0.5 - 0.6"
+ReversePropagation = "0.3 - 0.4"
 StaticArrays = "1"
-Symbolics = "5, 6"
-julia = "1"
+Symbolics = "5 - 7"
+julia = "1.10"
 
 [extras]
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
diff --git a/README.md b/README.md
index 342dd71..bd3558e 100644
--- a/README.md
+++ b/README.md
@@ -52,6 +52,10 @@ plot!(collect.(boundary), aspectratio=1, lw=0, label="boundary")
 
 
 
+## Architecture
+
+See [docs/src/architecture.md](docs/src/architecture.md) for a detailed description of the internal pipeline from symbolic expressions to contractors and separators.
+
 ## Author
 
 - [David P. Sanders](http://sistemas.fciencias.unam.mx/~dsanders),
diff --git a/REQUIRE b/REQUIRE
deleted file mode 100644
index f4b7bcc..0000000
--- a/REQUIRE
+++ /dev/null
@@ -1,6 +0,0 @@
-julia 1.0
-ModelingToolkit
-IntervalArithmetic 0.15
-IntervalRootFinding 0.4
-IntervalContractors 0.3
-MacroTools 0.4
diff --git a/benchmark/bench_dag_vs_codegen.jl b/benchmark/bench_dag_vs_codegen.jl
new file mode 100644
index 0000000..c805808
--- /dev/null
+++ b/benchmark/bench_dag_vs_codegen.jl
@@ -0,0 +1,107 @@
+using IntervalConstraintProgramming
+using IntervalConstraintProgramming: SharedDAG, add_expression!, propagate!
+using IntervalArithmetic, IntervalArithmetic.Symbols
+using IntervalBoxes
+using Symbolics
+using BenchmarkTools
+
+@variables x, y, z
+
+println("=" ^ 60)
+println("Benchmark: DAG (persistent) vs code-generation")
+println("=" ^ 60)
+
+# --- Problem 1: Unit disk (x² + y² ≤ 1) ---
+
+println("\n--- Problem 1: Unit disk (x² + y² ≤ 1) ---\n")
+
+S_codegen = Separator(x^2 + y^2 <= 1, [x, y])
+S_dag = DAGSeparator(x^2 + y^2 <= 1, [x, y])
+X = IntervalBox(interval(-3, 3), 2)
+
+println("Single separator call on [-3,3]²:")
+print("  codegen: ")
+@btime $S_codegen($X)
+print("  DAG:     ")
+@btime $S_dag($X)
+
+for ϵ in [1.0, 0.1, 0.01]
+    println("\npave with ϵ = $ϵ:")
+    print("  codegen: ")
+    r1 = @btime pave($X, $S_codegen, $ϵ)
+    print("  DAG:     ")
+    r2 = @btime pave($X, $S_dag, $ϵ)
+    inner1, bnd1 = r1
+    inner2, bnd2 = r2
+    println("  codegen: $(length(inner1)) inner, $(length(bnd1)) boundary")
+    println("  DAG:     $(length(inner2)) inner, $(length(bnd2)) boundary")
+end
+
+# --- Problem 2: Annulus (1 ≤ x² + y² ≤ 4) ---
+
+println("\n--- Problem 2: Annulus (1 ≤ x² + y² ≤ 4) ---\n")
+
+annulus_expr = x^2 + y^2 ∈ interval(1, 4)
+S_codegen2 = separator(annulus_expr, [x, y])
+S_dag2 = DAGSeparator(annulus_expr, [x, y])
+
+X2 = IntervalBox(interval(-5, 5), 2)
+for ϵ in [1.0, 0.1]
+    println("pave with ϵ = $ϵ:")
+    print("  codegen: ")
+    r1 = @btime pave($X2, $S_codegen2, $ϵ)
+    print("  DAG:     ")
+    r2 = @btime pave($X2, $S_dag2, $ϵ)
+    inner1, bnd1 = r1
+    inner2, bnd2 = r2
+    println("  codegen: $(length(inner1)) inner, $(length(bnd1)) boundary")
+    println("  DAG:     $(length(inner2)) inner, $(length(bnd2)) boundary")
+end
+
+# --- Problem 3: 3D constraint ---
+
+println("\n--- Problem 3: 3D (3 - sqrt(x²+y²))² + z² ≤ 1) ---\n")
+
+S_codegen3 = Separator(3 - sqrt(x^2 + y^2)^2 + z^2 <= 1, [x, y, z])
+S_dag3 = DAGSeparator(3 - sqrt(x^2 + y^2)^2 + z^2 <= 1, [x, y, z])
+
+X3 = IntervalBox(interval(-10, 10), 3)
+println("pave with ϵ = 1.0:")
+print("  codegen: ")
+r1 = @btime pave($X3, $S_codegen3, 1.0)
+print("  DAG:     ")
+r2 = @btime pave($X3, $S_dag3, 1.0)
+inner1, bnd1 = r1
+inner2, bnd2 = r2
+println("  codegen: $(length(inner1)) inner, $(length(bnd1)) boundary")
+println("  DAG:     $(length(inner2)) inner, $(length(bnd2)) boundary")
+
+# --- Problem 4: Multi-constraint on shared DAG ---
+
+println("\n--- Problem 4: Multi-constraint shared DAG ---")
+println("    x²+y² ≤ 1 AND x²-y ≤ 0 (jointly propagated)\n")
+
+dag_shared = SharedDAG([x, y])
+add_expression!(dag_shared, x^2 + y^2 - 1, interval(-Inf, 0.0))
+add_expression!(dag_shared, x^2 - y, interval(-Inf, 0.0))
+
+dag_sep1 = SharedDAG([x, y])
+add_expression!(dag_sep1, x^2 + y^2 - 1, interval(-Inf, 0.0))
+dag_sep2 = SharedDAG([x, y])
+add_expression!(dag_sep2, x^2 - y, interval(-Inf, 0.0))
+
+X4 = IntervalBox(interval(-10, 10), 2)
+
+println("Single propagation call on [-10,10]²:")
+print("  shared DAG (joint):      ")
+r_shared = @btime propagate!($dag_shared, $X4)
+print("  separate DAGs (seq):     ")
+r_sep = @btime begin
+    r1 = propagate!($dag_sep1, $X4)
+    propagate!($dag_sep2, r1)
+end
+println("  shared result:  ", r_shared)
+println("  separate result: ", r_sep)
+
+println("\nShared DAG nodes: ", length(dag_shared.nodes),
+        " vs separate: ", length(dag_sep1.nodes) + length(dag_sep2.nodes))
diff --git a/docs/Manifest.toml b/docs/Manifest.toml
deleted file mode 100644
index 9f9b3b2..0000000
--- a/docs/Manifest.toml
+++ /dev/null
@@ -1,145 +0,0 @@
-# This file is machine-generated - editing it directly is not advised
-
-[[Base64]]
-uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
-
-[[CRlibm]]
-deps = ["Libdl"]
-git-tree-sha1 = "9d1c22cff9c04207f336b8e64840d0bd40d86e0e"
-uuid = "96374032-68de-5a5b-8d9e-752f78720389"
-version = "0.8.0"
-
-[[Dates]]
-deps = ["Printf"]
-uuid = "ade2ca70-3891-5945-98fb-dc099432e06a"
-
-[[Distributed]]
-deps = ["Random", "Serialization", "Sockets"]
-uuid = "8ba89e20-285c-5b6f-9357-94700520ee1b"
-
-[[DocStringExtensions]]
-deps = ["LibGit2", "Markdown", "Pkg", "Test"]
-git-tree-sha1 = "88bb0edb352b16608036faadcc071adda068582a"
-uuid = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
-version = "0.8.1"
-
-[[Documenter]]
-deps = ["Base64", "Dates", "DocStringExtensions", "InteractiveUtils", "JSON", "LibGit2", "Logging", "Markdown", "REPL", "Test", "Unicode"]
-git-tree-sha1 = "646ebc3db49889ffeb4c36f89e5d82c6a26295ff"
-uuid = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
-version = "0.24.7"
-
-[[ErrorfreeArithmetic]]
-git-tree-sha1 = "a5198ab6c8a724dd3965b31ddd11ccde65300f5b"
-uuid = "90fa49ef-747e-5e6f-a989-263ba693cf1a"
-version = "0.5.0"
-
-[[FastRounding]]
-deps = ["ErrorfreeArithmetic", "Test"]
-git-tree-sha1 = "224175e213fd4fe112db3eea05d66b308dc2bf6b"
-uuid = "fa42c844-2597-5d31-933b-ebd51ab2693f"
-version = "0.2.0"
-
-[[InteractiveUtils]]
-deps = ["Markdown"]
-uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
-
-[[IntervalArithmetic]]
-deps = ["CRlibm", "FastRounding", "LinearAlgebra", "Markdown", "Random", "RecipesBase", "SetRounding", "StaticArrays"]
-git-tree-sha1 = "b2db6ee367b4eb3ee8b009ede8ca809e4fd23d35"
-uuid = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
-version = "0.16.7"
-
-[[JSON]]
-deps = ["Dates", "Mmap", "Parsers", "Unicode"]
-git-tree-sha1 = "b34d7cef7b337321e97d22242c3c2b91f476748e"
-uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
-version = "0.21.0"
-
-[[LibGit2]]
-deps = ["Printf"]
-uuid = "76f85450-5226-5b5a-8eaa-529ad045b433"
-
-[[Libdl]]
-uuid = "8f399da3-3557-5675-b5ff-fb832c97cbdb"
-
-[[LinearAlgebra]]
-deps = ["Libdl"]
-uuid = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-
-[[Logging]]
-uuid = "56ddb016-857b-54e1-b83d-db4d58db5568"
-
-[[Markdown]]
-deps = ["Base64"]
-uuid = "d6f4376e-aef5-505a-96c1-9c027394607a"
-
-[[Mmap]]
-uuid = "a63ad114-7e13-5084-954f-fe012c677804"
-
-[[Parsers]]
-deps = ["Dates", "Test"]
-git-tree-sha1 = "75d07cb840c300084634b4991761886d0d762724"
-uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
-version = "1.0.1"
-
-[[Pkg]]
-deps = ["Dates", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "REPL", "Random", "SHA", "UUIDs"]
-uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
-
-[[Printf]]
-deps = ["Unicode"]
-uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"
-
-[[REPL]]
-deps = ["InteractiveUtils", "Markdown", "Sockets"]
-uuid = "3fa0cd96-eef1-5676-8a61-b3b8758bbffb"
-
-[[Random]]
-deps = ["Serialization"]
-uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
-
-[[RecipesBase]]
-git-tree-sha1 = "b4ed4a7f988ea2340017916f7c9e5d7560b52cae"
-uuid = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
-version = "0.8.0"
-
-[[SHA]]
-uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"
-
-[[Serialization]]
-uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
-
-[[SetRounding]]
-deps = ["Test"]
-git-tree-sha1 = "faca28c7937138d560ede5bfef2076d0cdf3584b"
-uuid = "3cc68bcd-71a2-5612-b932-767ffbe40ab0"
-version = "0.2.0"
-
-[[Sockets]]
-uuid = "6462fe0b-24de-5631-8697-dd941f90decc"
-
-[[SparseArrays]]
-deps = ["LinearAlgebra", "Random"]
-uuid = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
-
-[[StaticArrays]]
-deps = ["LinearAlgebra", "Random", "Statistics"]
-git-tree-sha1 = "5a3bcb6233adabde68ebc97be66e95dcb787424c"
-uuid = "90137ffa-7385-5640-81b9-e52037218182"
-version = "0.12.1"
-
-[[Statistics]]
-deps = ["LinearAlgebra", "SparseArrays"]
-uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
-
-[[Test]]
-deps = ["Distributed", "InteractiveUtils", "Logging", "Random"]
-uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
-
-[[UUIDs]]
-deps = ["Random", "SHA"]
-uuid = "cf7118a7-6976-5b1a-9a39-7adc72f591a4"
-
-[[Unicode]]
-uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
diff --git a/docs/Project.toml b/docs/Project.toml
index f64f5b3..d91698c 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -1,3 +1,9 @@
 [deps]
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+DocumenterVitepress = "4710194d-e776-4893-9690-8d956a29c365"
 IntervalArithmetic = "d1acc4aa-44c8-5952-acd4-ba5d80a2a253"
+IntervalConstraintProgramming = "138f1668-1576-5ad7-91b9-7425abbf3153"
+
+[compat]
+Documenter = "1"
+DocumenterVitepress = "0.3"
diff --git a/docs/make.jl b/docs/make.jl
index 18f9c73..4e5f98e 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -1,22 +1,28 @@
 using Documenter
+using DocumenterVitepress
 using IntervalConstraintProgramming, IntervalArithmetic
 
-makedocs(
+makedocs(;
     modules = [IntervalConstraintProgramming],
-    doctest = true,
-    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),
     authors = "David P. Sanders",
     sitename = "IntervalConstraintProgramming.jl",
-
-    pages = Any[
+    format = DocumenterVitepress.MarkdownVitepress(;
+        repo = "github.com/JuliaIntervals/IntervalConstraintProgramming.jl",
+        devbranch = "master",
+        devurl = "dev",
+    ),
+    pages = [
         "Home" => "index.md",
-        "API" => "api.md"
-    ]
-    )
+        "Contractors and Separators" => "contractors_and_separators.md",
+        "Architecture" => "architecture.md",
+        "API" => "api.md",
+    ],
+)
 
-deploydocs(
-    repo = "github.com/JuliaIntervals/IntervalConstraintProgramming.jl.git",
-    target = "build",
-    deps = nothing,
-    make = nothing,
+DocumenterVitepress.deploydocs(;
+    repo = "github.com/JuliaIntervals/IntervalConstraintProgramming.jl",
+    target = joinpath(@__DIR__, "build"),
+    branch = "gh-pages",
+    devbranch = "master",
+    push_preview = true,
 )
diff --git a/docs/mkdocs.yml b/docs/mkdocs.yml
deleted file mode 100644
index 860e60b..0000000
--- a/docs/mkdocs.yml
+++ /dev/null
@@ -1,24 +0,0 @@
-site_name:        PACKAGE_NAME.jl
-repo_url:         https://github.com/USER_NAME/PACKAGE_NAME.jl
-site_description: Description...
-site_author:      USER_NAME
-
-theme: readthedocs
-
-extra_css:
-  - assets/Documenter.css
-
-extra_javascript:
-  - https://cdn.mathjax.org/mathjax/latest/MathJax.js?config=TeX-AMS-MML_HTMLorMML
-  - assets/mathjaxhelper.js
-
-markdown_extensions:
-  - extra
-  - tables
-  - fenced_code
-  - mdx_math
-
-docs_dir: 'build'
-
-pages:
-  - Home: index.md
diff --git a/docs/package.json b/docs/package.json
new file mode 100644
index 0000000..42eda9b
--- /dev/null
+++ b/docs/package.json
@@ -0,0 +1,14 @@
+{
+  "devDependencies": {
+    "@types/node": "^22.0.0"
+  },
+  "scripts": {
+    "docs:dev": "vitepress dev build/.documenter",
+    "docs:build": "vitepress build build/.documenter",
+    "docs:preview": "vitepress preview build/.documenter"
+  },
+  "dependencies": {
+    "vitepress": "^1.6.4",
+    "@mdit/plugin-mathjax": "^0.26.1"
+  }
+}
diff --git a/docs/src/architecture.md b/docs/src/architecture.md
new file mode 100644
index 0000000..8147561
--- /dev/null
+++ b/docs/src/architecture.md
@@ -0,0 +1,157 @@
+# Architecture: From Symbolic Expressions to Contractors
+
+This document describes the internal pipeline that converts symbolic mathematical expressions into interval contractors and separators for constraint propagation.
+
+## Overview
+
+IntervalConstraintProgramming works in three stages:
+
+1. **Parse**: A symbolic constraint (e.g. `x^2 + y^2 <= 1`) is normalized to the form `f(x,y) in [a, b]`.
+2. **Compile**: The expression `f` is decomposed into SSA (single static assignment) form via common subexpression elimination (CSE), then a forward--backward contractor is generated using [ReversePropagation.jl](https://github.com/JuliaIntervals/ReversePropagation.jl).
+3. **Evaluate**: The separator applies the contractor to interval boxes, classifying regions as *inner* (definitely feasible), *outer* (definitely infeasible), or *boundary* (unknown).
+
+```
+Symbolic expression  (Symbolics.jl)
+        |
+   separator(ex, vars)
+        |
+        v
+  Normalize: f(x) in constraint
+        |
+        +---> make_function(f, vars)          --> evaluation function
+        +---> Contractor(f, vars)
+                  |
+                  v
+           forward_backward_contractor(f, vars)   [ReversePropagation]
+                  |
+                  +---> cse_equations(f)           --> SSA assignments
+                  +---> constraint intersection
+                  +---> rev.() on each assignment  --> backward propagation
+                  +---> eval()                     --> compiled Julia closure
+                  |
+                  v
+             Separator(vars, f, constraint, eval_fn, contractor)
+```
+
+## Stage 1: Parsing and normalization
+
+The entry point is `separator(ex)` (or equivalently `constraint(ex)`), defined in `src/utils.jl`.
+
+Given a symbolic expression, it:
+
+- **Detects the operator** (`<=`, `>=`, `==`, `&`, `|`, `!`) at the top level.
+- **Handles logical connectives** recursively:
+  - `A & B` becomes `separator(A)` meet `separator(B)`
+  - `A | B` becomes `separator(A)` join `separator(B)`
+  - `!A` swaps inner and outer of `separator(A)`
+- **Normalizes inequalities** by moving everything to one side:
+  - `lhs <= rhs` becomes `lhs - rhs` with constraint `[-Inf, 0]`
+  - `lhs >= rhs` becomes `lhs - rhs` with constraint `[0, +Inf]`
+  - `lhs == rhs` becomes `lhs - rhs` with constraint `[0, 0]`
+
+The `in` operator is handled symbolically:
+```julia
+Base.in(x::Num, y::Interval) = (x >= Num(inf(y))) & (x <= Num(sup(y)))
+```
+so `x^2 + y^2 in interval(0, 1)` decomposes into two inequality constraints combined with `&`.
+
+## Stage 2: Compilation via ReversePropagation
+
+### Common Subexpression Elimination (CSE)
+
+The function `cse_equations(ex)` in ReversePropagation (`src/cse.jl`) recursively decomposes the expression into binary operations, assigning a fresh variable to each subexpression:
+
+```
+Input:  3x^2 + 4x^2*y
+
+Output (SSA assignments):
+  _a := x * x
+  _b := 3 * _a
+  _c := _a * y
+  _d := 4 * _c
+  _e := _b + _d
+```
+
+Common subexpressions (like `x * x`) are computed only once.
+
+### Forward--backward contractor (HC4Revise)
+
+The function `forward_backward_code(ex, vars)` in `src/reverse_icp.jl` generates the full contractor:
+
+1. **Forward pass**: The CSE assignments evaluate `f` with interval arithmetic.
+2. **Constraint intersection**: The result is intersected with the constraint interval:
+   ```
+   value := _e
+   _e := _e meet constraint
+   ```
+3. **Backward pass**: Each CSE assignment is reversed using `rev()`, which calls the corresponding reverse function from [IntervalContractors.jl](https://github.com/JuliaIntervals/IntervalContractors.jl) (e.g. `plus_rev`, `mul_rev`, `sin_rev`). These propagate the constraint back through each operation, contracting the input intervals.
+
+The resulting code is compiled into a Julia closure via `eval()`.
+
+### Reverse operations
+
+For each elementary operation, IntervalContractors provides a *reverse* (or *backward*) function:
+
+| Forward         | Reverse         | Effect                                    |
+|-----------------|-----------------|-------------------------------------------|
+| `z := x + y`    | `plus_rev`      | Contract `x` and `y` given constraint on `z` |
+| `z := x * y`    | `mul_rev`       | Contract `x` and `y` given constraint on `z` |
+| `z := sin(x)`   | `sin_rev`       | Contract `x` given constraint on `z`      |
+| `z := x ^ n`    | `power_rev`     | Contract `x` given constraint on `z`      |
+
+These are registered as symbolic functions so they can appear in generated code.
+
+## Stage 3: Separator evaluation
+
+A `Separator` stores both the compiled evaluation function and the contractor. When called with an `IntervalBox`:
+
+```julia
+function (SS::Separator)(X::IntervalBox)
+    # 1. Contract: find the boundary region
+    boundary = SS.contractor(X)
+
+    # 2. Evaluate at corners to classify inner/outer
+    lb_image = SS.f(inf.(X))
+    ub_image = SS.f(sup.(X))
+
+    # 3. If a corner satisfies the constraint, extend inner; otherwise extend outer
+    inner = boundary
+    outer = boundary
+    # ... extend with corner evaluations ...
+
+    return (boundary, inner, outer)
+end
+```
+
+For boxes extending to infinity, a different strategy is used: the contractor is applied once with the constraint and once with its complement, and inner/outer are determined from those two contractions.
+
+### Separator composition
+
+Separators compose via set-theoretic operations defined in `src/set_operations.jl`:
+
+| Operation | Symbol | Inner               | Outer               |
+|-----------|--------|----------------------|----------------------|
+| Intersect | `S1 meet S2`   | `inner1 meet inner2` | `outer1 join outer2` |
+| Union     | `S1 join S2`   | `inner1 join inner2` | `outer1 meet outer2` |
+| Complement| `!S`           | swap inner/outer     | swap inner/outer     |
+
+## Paving
+
+The `pave(X, S, epsilon)` function in `src/pave.jl` recursively subdivides the domain box `X`:
+
+1. Apply the separator to get `(boundary, inner, outer)`.
+2. Regions proven inner are accumulated.
+3. Regions proven outer are discarded.
+4. Boundary regions are bisected and re-processed until their diameter falls below `epsilon`.
+
+The result is two collections: `inner` boxes (guaranteed feasible) and `boundary` boxes (undetermined at this tolerance).
+
+## Gradient computation
+
+ReversePropagation also provides `gradient_code(ex, vars)` which generates SSA code for the gradient of `ex` using reverse-mode automatic differentiation (via ChainRules.jl):
+
+1. **Forward pass**: Same CSE decomposition as above.
+2. **Initialization**: Set the adjoint of the output to 1.
+3. **Adjoint pass**: For each forward assignment in reverse order, compute the adjoint contribution using `rrule` from ChainRules.
+
+The result is a sequence of assignments whose final variables hold the partial derivatives. This SSA code can itself be fed back into the forward--backward machinery to create contractors for derivative constraints.
diff --git a/docs/src/contractors_and_separators.md b/docs/src/contractors_and_separators.md
new file mode 100644
index 0000000..1153749
--- /dev/null
+++ b/docs/src/contractors_and_separators.md
@@ -0,0 +1,247 @@
+# Contractors and Separators
+
+This page explains the key concepts behind IntervalConstraintProgramming.jl: **contractors** and **separators**. These are the tools that allow us to rigorously determine which regions of space satisfy a given set of constraints.
+
+## Interval boxes
+
+An **interval box** (or simply "box") is a Cartesian product of intervals, representing a rectangular region in $n$-dimensional space. For example, in 2D:
+
+$$X = [1, 3] \times [0, 2]$$
+
+is the rectangle with $x \in [1, 3]$ and $y \in [0, 2]$.
+
+In Julia, we write:
+
+```julia
+using IntervalArithmetic, IntervalBoxes
+X = IntervalBox(interval(1, 3), interval(0, 2))
+```
+
+## Contractors
+
+A **contractor** $C$ is an operator that takes an interval box $X$ and returns a smaller (or equal) box $C(X) \subseteq X$. The key property is:
+
+> **Contractance**: $C(X) \subseteq X$ for all boxes $X$.
+
+Given an expression like $x^2 + y^2$ and a constraint interval like $[0, 1]$, a contractor shrinks the input box to eliminate regions that are *guaranteed* not to satisfy the constraint $x^2 + y^2 \in [0, 1]$.
+
+### How contraction works: forward--backward propagation
+
+IntervalConstraintProgramming uses a **forward--backward contractor**, also known as constraint propagation:
+
+1. **Forward pass**: Evaluate the expression using interval arithmetic. Each intermediate result gets an interval that encloses all possible values.
+
+2. **Constraint intersection**: Intersect the result with the constraint interval. If the expression evaluates to $[{-0.5}, 1.5]$ and the constraint is $[0, 1]$, the intersected output is $[0, 1]$.
+
+3. **Backward pass**: Propagate the tightened output *backwards* through each operation, using inverse interval operations. This tightens the intervals on the input variables.
+
+For example, if we know $z = x + y$, $z \in [0, 1]$, $x \in [{-1}, 3]$ and $y \in [{-1}, 3]$, then the backward step deduces $x \in [{-1}, 2]$ and $y \in [{-1}, 2]$, since $x$ cannot exceed $1 - ({-1}) = 2$ given the constraint on $z$.
+
+### Creating a contractor
+
+```julia
+using IntervalArithmetic, IntervalConstraintProgramming, Symbolics
+
+vars = @variables x, y
+C = Contractor(x^2 + y^2, vars)
+```
+
+Calling the contractor on a box returns the contracted box:
+
+```julia
+X = IntervalBox(-5..5, 2)
+
+C(X)                        # contract with default constraint [0, 0]
+C(X, interval(0, 1))        # contract with constraint [0, 1]
+```
+
+Contractors are useful on their own when you want to find the set where an expression takes a particular value (or range of values). However, they do not distinguish the *inside* of a set from the *outside* --- that is what separators do.
+
+
+## Separators
+
+A **separator** $S$ extends the idea of a contractor to classify regions as **inner** (satisfying the constraint), **outer** (violating the constraint), or **boundary** (undetermined).
+
+Given a constraint like $x^2 + y^2 \le 1$, a separator applied to a box $X$ returns three boxes:
+
+- **inner**: a box guaranteed to lie *inside* the feasible set (where the constraint is satisfied)
+- **outer**: a box guaranteed to contain the *outside* (where the constraint is violated)
+- **boundary**: a box containing the boundary of the constraint set, where the status is unknown
+
+```
+boundary, inner, outer = S(X)
+```
+
+The fundamental guarantee is:
+
+> Every point in $X$ that satisfies the constraint lies in `inner`.
+> Every point in $X$ that violates the constraint lies in `outer`.
+> The `boundary` box contains all points whose status is uncertain.
+
+### Creating a separator
+
+The main way to create a separator is with the `constraint` function (or equivalently `separator`):
+
+```julia
+using IntervalArithmetic, IntervalArithmetic.Symbols
+using IntervalConstraintProgramming, Symbolics
+
+vars = @variables x, y
+
+S = constraint(x^2 + y^2 <= 1, vars)
+```
+
+The constraint is automatically normalized: `x^2 + y^2 <= 1` becomes `x^2 + y^2 - 1` with constraint interval $[-\infty, 0]$.
+
+You can use any of the standard comparison operators:
+
+| Constraint syntax               | Meaning                                   |
+|---------------------------------|-------------------------------------------|
+| `constraint(expr <= val, vars)` | $\text{expr} \le \text{val}$              |
+| `constraint(expr >= val, vars)` | $\text{expr} \ge \text{val}$              |
+| `constraint(expr ~ val, vars)`  | $\text{expr} = \text{val}$ (equality)     |
+
+You can also use the `in` operator ($\in$) with an interval:
+
+```julia
+S = constraint(x^2 + y^2 in interval(1, 3), vars)
+```
+
+This is equivalent to $1 \le x^2 + y^2 \le 3$.
+
+### Applying a separator
+
+```julia
+X = IntervalBox(-5..5, 2)
+
+boundary, inner, outer = S(X)
+```
+
+For a large box like $[-5, 5]^2$, the contractor can only tighten slightly. The real power comes from combining separators with the **paving algorithm**.
+
+## Paving: finding the feasible set
+
+The `pave` function recursively bisects the domain and applies a separator at each step. Boxes that are proved to be inside the constraint set are collected as the **inner** approximation; boxes that cannot be classified at the given tolerance are collected as the **boundary**:
+
+```julia
+X = IntervalBox(-5..5, 2)
+S = constraint(x^2 + y^2 <= 1, vars)
+
+inner, boundary = pave(X, S, 0.05)
+```
+
+The third argument is the **tolerance**: bisection stops when a boundary box has diameter smaller than this value.
+
+The result gives rigorous guarantees:
+- Every point in an `inner` box satisfies the constraint.
+- Every point *not* covered by `inner` or `boundary` boxes violates the constraint.
+- Decreasing the tolerance gives a sharper approximation at the cost of more computation.
+
+### Plotting the result
+
+```julia
+using Plots
+
+plot(collect.(inner), aspectratio=1, lw=0, label="inner")
+plot!(collect.(boundary), aspectratio=1, lw=0, label="boundary")
+```
+
+## Combining separators: set operations
+
+Separators can be combined using set-theoretic operations to describe more complex feasible sets.
+
+### Intersection ($\sqcap$)
+
+The intersection of two separators gives the set of points satisfying *both* constraints:
+
+```julia
+S1 = constraint(x^2 + y^2 >= 1, vars)
+S2 = constraint(x^2 + y^2 <= 3, vars)
+
+S = S1 ⊓ S2    # type \sqcap<TAB> in Julia
+```
+
+This defines the annular region $1 \le x^2 + y^2 \le 3$.
+
+The inner contractor of `S1 ⊓ S2` is the intersection of the inner contractors, and the outer contractor is the union of the outer contractors:
+
+$$\text{inner}(S_1 \sqcap S_2) = \text{inner}(S_1) \sqcap \text{inner}(S_2)$$
+$$\text{outer}(S_1 \sqcap S_2) = \text{outer}(S_1) \sqcup \text{outer}(S_2)$$
+
+### Union ($\sqcup$)
+
+The union gives the set of points satisfying *either* constraint:
+
+```julia
+S = S1 ⊔ S2    # type \sqcup<TAB> in Julia
+```
+
+$$\text{inner}(S_1 \sqcup S_2) = \text{inner}(S_1) \sqcup \text{inner}(S_2)$$
+$$\text{outer}(S_1 \sqcup S_2) = \text{outer}(S_1) \sqcap \text{outer}(S_2)$$
+
+
+### Complement (`!`)
+
+The complement swaps inner and outer:
+
+```julia
+S_complement = !S
+```
+
+### Set difference and symmetric difference
+
+```julia
+S_diff = setdiff(S1, S2)      # S1 and not S2
+S_sym  = symdiff(S1, S2)      # in one but not both
+```
+
+These are defined in terms of the basic operations:
+- `setdiff(S1, S2)` = `S1 ⊓ !S2`
+- `symdiff(S1, S2)` = `setdiff(S1, S2) ⊔ setdiff(S2, S1)`
+
+## Complete example
+
+Here is a complete example that finds the intersection of an ellipse and a rotated constraint:
+
+```julia
+using IntervalArithmetic, IntervalArithmetic.Symbols
+using IntervalConstraintProgramming
+using IntervalBoxes
+using Symbolics
+
+vars = @variables x, y
+
+# An elliptical region
+C1 = constraint(x^2 + 2y^2 >= 1, vars)
+
+# A second constraint
+C2 = constraint(x^2 + y^2 + x * y <= 3, vars)
+
+# Intersection: points satisfying both constraints
+C = C1 ⊓ C2
+
+# Pave the domain
+X = IntervalBox(-5..5, 2)
+inner, boundary = pave(X, C, 0.05)
+
+# Visualize
+using Plots
+plot(collect.(inner), aspectratio=1, lw=0, label="inner")
+plot!(collect.(boundary), aspectratio=1, lw=0, label="boundary")
+```
+
+## How it all fits together
+
+The pipeline from a symbolic constraint to a paving is:
+
+1. **Symbolic expression** $\to$ normalized form $f(\mathbf{x}) \in [a, b]$
+2. **Forward--backward contractor** compiled from the expression via [ReversePropagation.jl](https://github.com/JuliaIntervals/ReversePropagation.jl)
+3. **Separator** wraps the contractor and classifies regions as inner, outer, or boundary
+4. **Paving algorithm** recursively bisects and applies the separator to build a rigorous approximation of the feasible set
+
+For more details on the internal compilation pipeline, see the [Architecture](architecture.md) page.
+
+## References
+
+- Jaulin, L., Kieffer, M., Didrit, O., and Walter, E. (2001). *Applied Interval Analysis*. Springer.
+- Jaulin, L. and Desrochers, B. (2014). Introduction to the Algebra of Separators with Application to Path Planning. *Engineering Applications of Artificial Intelligence*, **33**, 141--147.
diff --git a/docs/src/index.md b/docs/src/index.md
index 3dd5ede..d2d4686 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -1,214 +1,64 @@
-# `IntervalConstraintProgramming.jl`
+# IntervalConstraintProgramming.jl
 
 This Julia package allows you to specify a set of **constraints** on real-valued variables, given by (in)equalities, and
 rigorously calculate inner and outer approximations to the *feasible set*,
 i.e. the set that satisfies the constraints.
 
-This uses interval arithmetic provided by the author's
-[`IntervalArithmetic.jl`](https://github.com/dpsanders/IntervalArithmetic.jl) package,
+This uses interval arithmetic provided by
+[`IntervalArithmetic.jl`](https://github.com/JuliaIntervals/IntervalArithmetic.jl),
 in particular multi-dimensional `IntervalBox`es, i.e. Cartesian products of one-dimensional intervals.
 
 To do this, *interval constraint programming* is used, in particular the
 so-called "forward--backward contractor". This is implemented in terms of *separators*; see
-[Jaulin & Desrochers].
+[Jaulin & Desrochers (2014)](https://doi.org/10.1016/j.engappai.2014.04.010).
 
-```@meta
-DocTestSetup = quote
-    using IntervalConstraintProgramming, IntervalArithmetic
-end
-```
-
-## Usage
-Let's define a constraint, using the `@constraint` macro:
-```jldoctest
-julia> using IntervalConstraintProgramming, IntervalArithmetic
-
-julia> S = @constraint x^2 + y^2 <= 1
-Separator:
-  - variables: x, y
-  - expression: x ^ 2 + y ^ 2 ∈ [-Inf, 1]
-```
-It works out automatically that `x` and `y` are variables.
-The macro creates a `Separator` object, in this case a `ConstraintSeparator`.
-
-We now create an initial interval box in the ``x``--``y`` plane:
-```julia
-julia> x = y = -100..100   # notation for creating an interval with `IntervalArithmetic.jl`
-
-julia> X = IntervalBox(x, y)
-```
-
-The `@constraint` macro defines an object `S`, of type `Separator`.
-This is a function which,
-when applied to the box ``X = x \times y``
-in the x--y plane, applies two *contractors*, an inner one and an outer one.
-
-The inner contractor tries to shrink down, or *contract*, the box, to the smallest subbox
-of ``X`` that contains the part of ``X`` that satisfies the constraint; the
-outer contractor tries to contract ``X`` to the smallest subbox that contains the
-region where the constraint is not satisfied.
-
-When `S` is applied to the box `X`, it returns the result of the inner and outer contractors:
-```julia
-julia> inner, outer = S(X);
-
-julia> inner
-([-1, 1],[-1, 1])
-
-julia> outer
-([-100, 100],[-100, 100])
-```
-
-### Without using macros
-
-We can also make an object `S`, of type `Separator` or `C`, of type `Contractor` without using Macros, for that you need to define variables using `ModelingToolkit.jl`.
-Example  
-
-```julia
-julia> using IntervalConstraintProgramming, IntervalArithmetic, ModelingToolkit
-
-julia> @variables x y
-(x(), y())
-
-julia> S = Separator(x+y<1)
-Separator:
-  - variables: x, y
-  - expression: x() + y() == [-Inf, 1]
-
-julia> C = Contractor(x+y)
- Contractor in 2 dimensions:
-   - forward pass contracts to 1 dimensions
-   - variables: Symbol[:x, :y]
-   - expression: x() + y()
-```
-
-While making `Separator`or `Contractor`'s object we can also specify variables, like this
-
-```julia
-julia> vars = @variables x y z
-(x(), y(), z())
-
-julia> S = Separator(vars, x+y<1)
-Separator:
-  - variables: x, y, z
-  - expression: x() + y() == [-Inf, 1]
-
-julia> C = Contractor(vars, y+z)
-Contractor in 3 dimensions:
-  - forward pass contracts to 1 dimensions
-  - variables: Symbol[:x, :y, :z]
-  - expression: y() + z()
-```
-We can make objects (of type `Separator` or `Contractor`)by just using function name (Note: you have to specify variables explicitly as discussed above when you make objects by using function name). We can also use polynomial function to make objects.
+## Quick start
 
 ```julia
-julia> vars=@variables x y
-(x(), y())
+using IntervalArithmetic, IntervalArithmetic.Symbols
+using IntervalConstraintProgramming
+using IntervalBoxes
+using Symbolics
 
-julia> f(a,b)= a+b
-f (generic function with 1 method)
+vars = @variables x, y
 
-julia> C = Contractor(vars,f)
-Contractor in 2 dimensions:
-  - forward pass contracts to 1 dimensions
-  - variables: Symbol[:x, :y]
-  - expression: x() + y()
+# Create separators from constraints
+C1 = constraint(x^2 + 2y^2 >= 1, vars)
+C2 = constraint(x^2 + y^2 + x * y <= 3, vars)
 
-julia> f(a,b) = a+b <1
-f (generic function with 1 method)
+# Combine with intersection
+C = C1 ⊓ C2
 
-julia> S=Separator(vars, f)
-Separator:
-  - variables: x, y
-  - expression: x() + y() == [-Inf, 1]  
-
-julia> using DynamicPolynomials       #using polynomial functions
-
-julia> pvars = @polyvar x y
-(x, y)
-
-julia> f(a,b) = a + b < 1
-p (generic function with 1 method)
-
-julia> S=Separator(pvars, f)
-Separator:
-  - variables: x, y
-  - expression: x() + y() == [-Inf, 1]
-```
-#### BasicContractor
-Objects of type `Contractor` have four fields (variables, forward, backward and expression), among them data of two fields (forward, backward) are useful (i.e forward and backward functions) for further usage of that object, thats why it is preferred to use an object of type `BasicContractor` in place of `Contractor` which only contain these two fields for less usage of memory by unloading all the extra stuff.(Note: Like object of `Contractor` type,`BasicContractor`'s object will also have all the properties which are discussed above).
-
-```julia
-julia> @variables x y
-(x(), y())
-
-julia> C = BasicContractor(x^2 + y^2)
- Basic version of Contractor
+# Pave the domain
+X = IntervalBox(-5..5, 2)
+inner, boundary = pave(X, C, 0.05)
 ```
 
-
-
-## Set inversion: finding the feasible set
-
-To make progress, we must recursively bisect and apply the contractors, keeping
-track of the region proved to be inside the feasible set, and the region that is
-on the boundary ("both inside and outside"). This is done by the `pave` function,
-that takes a separator, a domain to search inside, and an optional tolerance:
+## Plotting the result
 
 ```julia
-julia> using Plots
+using Plots
 
-julia> x = y = -100..100
-
-julia> S = @constraint 1 <= x^2 + y^2 <= 3
-
-julia> paving = pave(S, X, 0.125);
+plot(collect.(inner), aspectratio=1, lw=0, label="inner")
+plot!(collect.(boundary), aspectratio=1, lw=0, label="boundary")
 ```
 
-`pave` returns an object of type `Paving`. This contains: the separator itself;
-an `inner` approximation, of type `SubPaving`, which is an alias for a `Vector` of `IntervalBox`es;
-a `SubPaving` representing the boxes on the boundary that could not be assigned either to the inside or outside of the set;
-and the tolerance.
-
-We may draw the result using a plot recipe from `IntervalArithmetic`. Either a
-single `IntervalBox`, or a `Vector` of `IntervalBox`es (which a `SubPaving` is)
-maybe be drawn using `plot` from `Plots.jl`:
-```julia
-julia> plot(paving.inner, c="green")
-julia> plot!(paving.boundary, c="gray")
-```
-
-The output should look something like this:
-
-![Ring](ring.png)
-
-
-The green boxes have been **rigorously** proved to be contained within the feasible set,
-and the white boxes to be outside the set. The grey boxes show those that lie on the boundary, whose status is unknown.
-
-### 3D
-
-The package works in any number of dimensions, although it suffers from the usual exponential slowdown in higher dimensions ("combinatorial explosion"); in 3D, it is still relatively fast.
-
-There are sample 3D calculations in the `examples` directory, in particular in the [solid torus notebook](examples/Solid torus.ipynb), which uses [`GLVisualize.gl`](https://github.com/JuliaGL/GLVisualize.jl) to provide an interactive visualization that may be rotated and zoomed. The output for the solid torus looks like this:
-
-![Coloured solid torus](solid_torus.png)
-
+The inner (blue) region is **rigorously** guaranteed to lie *inside* the feasible set.
+The white region is guaranteed to lie *outside* the set.
+The boundary (red) region lies on the boundary, where the status is unknown at the given tolerance.
 
 ## Set operations
-Separators may be combined using the operators `!` (complement), `⊓` and `⊔` to make
-more complicated sets; see the [notebook](https://github.com/JuliaIntervals/IntervalConstraintProgrammingNotebooks/blob/master/Basic%20examples%20of%20separators.ipynb) for several examples. Further examples can be found in the repository [IntervalConstraintProgrammingNotebooks](https://github.com/JuliaIntervals/IntervalConstraintProgrammingNotebooks).
 
-## Author
+Separators may be combined using the operators `!` (complement), `⊓` (intersection) and `⊔` (union) to build
+more complicated sets. See the [Contractors and Separators](@ref) page for details.
 
-- [David P. Sanders](http://sistemas.fciencias.unam.mx/~dsanders)
-    - [Julia lab, MIT](http://julia.mit.edu/)
-    - Departamento de Física, Facultad de Ciencias, Universidad Nacional Autónoma de México (UNAM)
+## Learn more
 
+- [Contractors and Separators](contractors_and_separators.md) -- detailed explanation of the key concepts
+- [Architecture](architecture.md) -- how the internal pipeline works
+- [API](api.md) -- function and type reference
 
-## References:
-- *Applied Interval Analysis*, Luc Jaulin, Michel Kieffer, Olivier Didrit, Eric Walter (2001)
-- Introduction to the Algebra of Separators with Application to Path Planning, Luc Jaulin and Benoît Desrochers, *Engineering Applications of Artificial Intelligence* **33**, 141–147 (2014)
+## References
 
-## Acknowledements
-Financial support is acknowledged from DGAPA-UNAM PAPIME grants PE-105911 and PE-107114, and DGAPA-UNAM PAPIIT grant IN-117214, and from a CONACYT-Mexico sabbatical fellowship. The author thanks Alan Edelman and the Julia group for hospitality during his sabbatical visit. He also thanks Luc Jaulin and Jordan Ninin for the [IAMOOC](http://iamooc.ensta-bretagne.fr/) online course, which introduced him to this subject.
+- *Applied Interval Analysis*, Luc Jaulin, Michel Kieffer, Olivier Didrit, Eric Walter (2001)
+- Introduction to the Algebra of Separators with Application to Path Planning, Luc Jaulin and Benoît Desrochers, *Engineering Applications of Artificial Intelligence* **33**, 141--147 (2014)
diff --git a/future.md b/future.md
new file mode 100644
index 0000000..76ba180
--- /dev/null
+++ b/future.md
@@ -0,0 +1,18 @@
+# Future work
+
+## `@constraint` macro for chained comparisons
+
+Julia lowers `0 <= x^2+y^2 <= 1` to `(0 <= x^2+y^2) && (x^2+y^2 <= 1)`, and `&&` requires a `Bool`, so this fails with symbolic expressions.
+
+A `@constraint` macro could intercept the AST before lowering and convert chained comparisons into `&` (which works symbolically) or directly into the `∈` form:
+
+```julia
+@constraint 0 <= x^2 + y^2 <= 1    # → x^2+y^2 ∈ interval(0, 1)
+@constraint x^2 + y^2 <= 1          # → simple case, works as-is
+```
+
+The package already exports `@constraint` but it is not yet defined.
+
+## ReversePropagation method overwrite warnings
+
+ReversePropagation emits many "Method definition overwritten" warnings with Symbolics v7. Harmless but noisy — should be addressed upstream in ReversePropagation.
diff --git a/src/IntervalConstraintProgramming.jl b/src/IntervalConstraintProgramming.jl
index cca9f62..f6c6c6a 100644
--- a/src/IntervalConstraintProgramming.jl
+++ b/src/IntervalConstraintProgramming.jl
@@ -24,12 +24,16 @@ import Base:
 
 import IntervalArithmetic: mid, interval, emptyinterval, isinf, isinterior, hull, mince
 
-
 @register_symbolic ¬(x)
-@register_symbolic x ∈ y::Interval
 @register_symbolic x ∨ y
 @register_symbolic x ∧ y
 
+# We cannot register `x ∈ y::Interval` as a symbolic operation because
+# SymbolicUtils hash-consing requires isequal/hash, which Interval deliberately
+# does not define (IEEE 1788). Instead, decompose into comparisons that the
+# package already handles.
+Base.in(x::Num, y::Interval) = (x >= Num(inf(y))) & (x <= Num(sup(y)))
+
 export
     # BasicContractor,
     # @contractor,
@@ -50,6 +54,9 @@ export
 
 export ¬
 
+# DAG-based propagation types
+export DAGContractor, DAGSeparator, SharedDAG, CompiledDAG, add_expression!, compile
+
 const reverse_operations = IntervalContractors.reverse_operations
 
 
@@ -58,5 +65,12 @@ include("contractor.jl")
 include("set_operations.jl")
 include("pave.jl")
 
+# DAG-based alternative to code-generation approach
+include("dag/nodes.jl")
+include("dag/build.jl")
+include("dag/propagate.jl")
+include("dag/compile.jl")
+include("dag/contractor.jl")
+
 
 end # module
diff --git a/src/contractor.jl b/src/contractor.jl
index 4aefbf2..3181bbd 100644
--- a/src/contractor.jl
+++ b/src/contractor.jl
@@ -25,7 +25,13 @@ end
 
 # Base.show(io::IO, S::Separator) = print(io, "Separator($(S.ex) ∈ $(S.constraint), vars = $(join(S.vars, ", ")))")
 
-Base.show(io::IO, S::AbstractSeparator) = print(io, "Separator($(S.ex), vars=$(join(S.vars, ", ")))")
+function Base.show(io::IO, S::AbstractSeparator)
+    if hasproperty(S, :constraint)
+        print(io, "Separator($(S.ex) ∈ $(S.constraint), vars=$(join(S.vars, ", ")))")
+    else
+        print(io, "Separator($(S.ex), vars=$(join(S.vars, ", ")))")
+    end
+end
 
 function Separator(orig_expr, vars)
     ex, constraint = analyse(orig_expr)
@@ -33,7 +39,7 @@ function Separator(orig_expr, vars)
     return Separator(ex, vars, constraint)
 end
 
-Separator(ex, vars, constraint::Interval) = Separator(vars, ex ∈ constraint, constraint, make_function(ex, vars), Contractor(ex, vars))
+Separator(ex, vars, constraint::Interval) = Separator(vars, ex, constraint, make_function(ex, vars), Contractor(ex, vars))
 
 function separate_infinite_box(S::Separator, X::IntervalBox)
     # for an box that extends to infinity we cannot evaluate at a corner
diff --git a/src/dag/build.jl b/src/dag/build.jl
new file mode 100644
index 0000000..f756bd6
--- /dev/null
+++ b/src/dag/build.jl
@@ -0,0 +1,83 @@
+const SU = Symbolics.SymbolicUtils
+
+"""
+    add_expression!(dag::SharedDAG, expr::Num, constraint::Interval) → ConstraintEntry
+
+Add an expression with its constraint to an existing SharedDAG. Reuses
+variable nodes and any common subexpressions already in the DAG.
+
+Returns the ConstraintEntry that was added.
+"""
+function add_expression!(dag::SharedDAG, expr::Num, constraint::Interval)
+    sym_expr = Symbolics.value(expr)
+    root = _build_node!(sym_expr, dag)
+    entry = ConstraintEntry(root, constraint)
+    push!(dag.constraints, entry)
+    return entry
+end
+
+"""
+    build_dag(expr::Num, vars::Vector{Num}, constraint::Interval) → SharedDAG
+
+Convenience: build a SharedDAG with a single expression and constraint.
+"""
+function build_dag(expr::Num, vars::Vector{Num}, constraint::Interval)
+    dag = SharedDAG(vars)
+    add_expression!(dag, expr, constraint)
+    return dag
+end
+
+"""
+Recursively build DAG nodes from a symbolic expression, reusing existing
+nodes in the SharedDAG via the node_cache (CSE).
+"""
+function _build_node!(expr, dag::SharedDAG)
+    h = objectid(expr)
+    haskey(dag.node_cache, h) && return dag.node_cache[h]
+
+    node = if Symbolics.issym(expr)
+        name = Symbolics.nameof(expr)
+        if haskey(dag.var_map, name)
+            dag.var_map[name]
+        else
+            error("Unknown variable in expression: $name")
+        end
+
+    elseif SU.iscall(expr)
+        _build_operation!(expr, dag)
+
+    else
+        # Literal number
+        val = Symbolics.value(expr)
+        cn = ConstantNode(val)
+        push!(dag.nodes, cn)
+        cn
+    end
+
+    dag.node_cache[h] = node
+    return node
+end
+
+"""
+Build an OperationNode, decomposing n-ary `+` and `*` into binary chains.
+"""
+function _build_operation!(expr, dag::SharedDAG)
+    op = Symbolics.operation(expr)
+    args = Symbolics.arguments(expr)
+
+    if (op === (+) || op === (*)) && length(args) > 2
+        node = _build_node!(args[1], dag)
+        for i in 2:length(args)
+            right = _build_node!(args[i], dag)
+            parent = OperationNode(op, AbstractNode[node, right])
+            push!(dag.nodes, parent)
+            node = parent
+        end
+        return node
+    end
+
+    child_nodes = AbstractNode[_build_node!(a, dag) for a in args]
+    node = OperationNode(op, child_nodes)
+    push!(dag.nodes, node)
+    return node
+end
diff --git a/src/dag/compile.jl b/src/dag/compile.jl
new file mode 100644
index 0000000..890a6ba
--- /dev/null
+++ b/src/dag/compile.jl
@@ -0,0 +1,316 @@
+struct Instruction
+    op::Symbol
+    out::Int32      # output slot
+    arg1::Int32     # first input slot
+    arg2::Int32     # second input slot (0 for unary; integer exponent for :pow_int)
+end
+
+"""
+    CompiledDAG
+
+A compiled representation of a SharedDAG optimized for fast propagation.
+Converts the node graph into a flat instruction array operating on a
+pre-allocated workspace of intervals. Each DAG node (variable, constant,
+or operation) gets a slot in the workspace. Shared nodes (CSE) map to the
+same slot, so backward contraction from multiple parents accumulates
+correctly via intersection.
+
+This eliminates:
+- Dynamic dispatch on `node.op::Function` (replaced by Symbol comparison)
+- Per-node allocation (pre-allocated flat workspace)
+- Node traversal overhead (flat array iteration)
+"""
+struct CompiledDAG
+    instructions::Vector{Instruction}
+    workspace::Vector{Interval}
+    n_vars::Int
+    constant_slots::Vector{Int32}
+    constant_values::Vector{Interval}
+    root_slots::Vector{Int32}
+    constraints::Vector{Interval}
+end
+
+const _OP_SYMBOL = Dict{Function, Symbol}(
+    (+) => :add, (-) => :sub, (*) => :mul, (/) => :div, (^) => :pow,
+    min => :min, max => :max,
+    sin => :sin, cos => :cos, tan => :tan,
+    asin => :asin, acos => :acos, atan => :atan,
+    sinh => :sinh, cosh => :cosh, tanh => :tanh,
+    asinh => :asinh, acosh => :acosh, atanh => :atanh,
+    exp => :exp, exp2 => :exp2, exp10 => :exp10, expm1 => :expm1,
+    log => :log, log2 => :log2, log10 => :log10, log1p => :log1p,
+    sqrt => :sqrt, abs => :abs, sign => :sign,
+)
+
+"""
+    compile(dag::SharedDAG) → CompiledDAG
+
+Compile a SharedDAG into a flat instruction schedule. Each node in the DAG
+gets a workspace slot (variables first, then constants, then operations).
+Shared nodes get a single slot, preserving CSE.
+"""
+function compile(dag::SharedDAG)
+    node_to_slot = Dict{AbstractNode, Int32}()
+    slot = Int32(0)
+
+    for v in dag.variables
+        slot += Int32(1)
+        node_to_slot[v] = slot
+    end
+    n_vars = Int(slot)
+
+    constant_slots = Int32[]
+    constant_values = Interval[]
+
+    for node in dag.nodes
+        haskey(node_to_slot, node) && continue
+        slot += Int32(1)
+        node_to_slot[node] = slot
+        if node isa ConstantNode
+            push!(constant_slots, slot)
+            push!(constant_values, node.range)
+        end
+    end
+
+    n_slots = Int(slot)
+    workspace = fill(entireinterval(), n_slots)
+
+    for (i, cs) in enumerate(constant_slots)
+        workspace[cs] = constant_values[i]
+    end
+
+    instructions = Instruction[]
+    for node in dag.nodes
+        node isa OperationNode || continue
+        out = node_to_slot[node]
+        children = node.children
+
+        if length(children) == 1
+            a1 = node_to_slot[children[1]]
+            op_sym = node.op === (-) ? :neg : get(_OP_SYMBOL, node.op, :unknown)
+            push!(instructions, Instruction(op_sym, out, a1, Int32(0)))
+
+        elseif length(children) == 2
+            a1 = node_to_slot[children[1]]
+            a2 = node_to_slot[children[2]]
+
+            if node.op === (^) && children[2] isa ConstantNode
+                n = children[2].value
+                if n == round(Int, n)
+                    push!(instructions, Instruction(:pow_int, out, a1, Int32(round(Int, n))))
+                    continue
+                end
+            end
+
+            op_sym = get(_OP_SYMBOL, node.op, :unknown)
+            push!(instructions, Instruction(op_sym, out, a1, a2))
+        end
+    end
+
+    root_slots = Int32[node_to_slot[e.root] for e in dag.constraints]
+    constraints = Interval[e.constraint for e in dag.constraints]
+
+    return CompiledDAG(instructions, workspace, n_vars,
+                       constant_slots, constant_values, root_slots, constraints)
+end
+
+# --- Compiled forward pass ---
+
+function forward!(cd::CompiledDAG, X::IntervalBox)
+    ws = cd.workspace
+    @inbounds for i in 1:cd.n_vars
+        ws[i] = X[i]
+    end
+    @inbounds for (i, cs) in enumerate(cd.constant_slots)
+        ws[cs] = cd.constant_values[i]
+    end
+    @inbounds for inst in cd.instructions
+        ws[inst.out] = _fwd(inst.op, ws, inst.arg1, inst.arg2)
+    end
+end
+
+function _fwd(op::Symbol, ws::Vector{Interval}, a1::Int32, a2::Int32)
+    @inbounds begin
+        op === :add && return ws[a1] + ws[a2]
+        op === :sub && return ws[a1] - ws[a2]
+        op === :mul && return ws[a1] * ws[a2]
+        op === :div && return ws[a1] / ws[a2]
+        op === :pow && return ws[a1] ^ ws[a2]
+        op === :pow_int && return ws[a1] ^ Int(a2)
+        op === :neg && return -ws[a1]
+        op === :min && return IntervalArithmetic.min(ws[a1], ws[a2])
+        op === :max && return IntervalArithmetic.max(ws[a1], ws[a2])
+        op === :sin && return sin(ws[a1])
+        op === :cos && return cos(ws[a1])
+        op === :tan && return tan(ws[a1])
+        op === :asin && return asin(ws[a1])
+        op === :acos && return acos(ws[a1])
+        op === :atan && return atan(ws[a1])
+        op === :sinh && return sinh(ws[a1])
+        op === :cosh && return cosh(ws[a1])
+        op === :tanh && return tanh(ws[a1])
+        op === :asinh && return asinh(ws[a1])
+        op === :acosh && return acosh(ws[a1])
+        op === :atanh && return atanh(ws[a1])
+        op === :exp && return exp(ws[a1])
+        op === :exp2 && return exp2(ws[a1])
+        op === :exp10 && return exp10(ws[a1])
+        op === :expm1 && return expm1(ws[a1])
+        op === :log && return log(ws[a1])
+        op === :log2 && return log2(ws[a1])
+        op === :log10 && return log10(ws[a1])
+        op === :log1p && return log1p(ws[a1])
+        op === :sqrt && return sqrt(ws[a1])
+        op === :abs && return abs(ws[a1])
+        op === :sign && return sign(ws[a1])
+    end
+    return entireinterval()
+end
+
+# --- Compiled backward pass ---
+
+function backward!(cd::CompiledDAG)
+    ws = cd.workspace
+    @inbounds for i in eachindex(cd.root_slots)
+        rs = cd.root_slots[i]
+        ws[rs] = IntervalArithmetic.intersect_interval(ws[rs], cd.constraints[i])
+    end
+
+    if all(@inbounds(isempty_interval(ws[rs])) for rs in cd.root_slots)
+        @inbounds for i in 1:cd.n_vars
+            ws[i] = emptyinterval()
+        end
+        return
+    end
+
+    @inbounds for i in length(cd.instructions):-1:1
+        _bwd!(cd.instructions[i], ws)
+    end
+end
+
+function backward!(cd::CompiledDAG, constraint::Interval)
+    ws = cd.workspace
+    @inbounds begin
+        rs = cd.root_slots[1]
+        ws[rs] = IntervalArithmetic.intersect_interval(ws[rs], constraint)
+    end
+
+    if @inbounds isempty_interval(ws[cd.root_slots[1]])
+        @inbounds for i in 1:cd.n_vars
+            ws[i] = emptyinterval()
+        end
+        return
+    end
+
+    @inbounds for i in length(cd.instructions):-1:1
+        _bwd!(cd.instructions[i], ws)
+    end
+end
+
+function _bwd!(inst::Instruction, ws::Vector{Interval})
+    op = inst.op
+    out = inst.out
+    a1 = inst.arg1
+    a2 = inst.arg2
+    @inbounds z = ws[out]
+
+    if isempty_interval(z)
+        @inbounds ws[a1] = emptyinterval()
+        @inbounds a2 > 0 && op !== :pow_int && (ws[a2] = emptyinterval())
+        return
+    end
+
+    @inbounds begin
+        if op === :add
+            _, l, r = plus_rev(z, ws[a1], ws[a2])
+            ws[a1] = IntervalArithmetic.intersect_interval(ws[a1], l)
+            ws[a2] = IntervalArithmetic.intersect_interval(ws[a2], r)
+        elseif op === :sub
+            _, l, r = minus_rev(z, ws[a1], ws[a2])
+            ws[a1] = IntervalArithmetic.intersect_interval(ws[a1], l)
+            ws[a2] = IntervalArithmetic.intersect_interval(ws[a2], r)
+        elseif op === :mul
+            _, l, r = mul_rev(z, ws[a1], ws[a2])
+            ws[a1] = IntervalArithmetic.intersect_interval(ws[a1], l)
+            ws[a2] = IntervalArithmetic.intersect_interval(ws[a2], r)
+        elseif op === :div
+            _, l, r = div_rev(z, ws[a1], ws[a2])
+            ws[a1] = IntervalArithmetic.intersect_interval(ws[a1], l)
+            ws[a2] = IntervalArithmetic.intersect_interval(ws[a2], r)
+        elseif op === :pow
+            _, l, r = power_rev(z, ws[a1], ws[a2])
+            ws[a1] = IntervalArithmetic.intersect_interval(ws[a1], l)
+            ws[a2] = IntervalArithmetic.intersect_interval(ws[a2], r)
+        elseif op === :pow_int
+            _, x_new, _ = power_rev(z, ws[a1], Int(a2))
+            ws[a1] = IntervalArithmetic.intersect_interval(ws[a1], x_new)
+        elseif op === :neg
+            _, x = minus_rev(z, ws[a1])
+            ws[a1] = IntervalArithmetic.intersect_interval(ws[a1], x)
+        elseif op === :min
+            _, l, r = min_rev(z, ws[a1], ws[a2])
+            ws[a1] = IntervalArithmetic.intersect_interval(ws[a1], l)
+            ws[a2] = IntervalArithmetic.intersect_interval(ws[a2], r)
+        elseif op === :max
+            _, l, r = max_rev(z, ws[a1], ws[a2])
+            ws[a1] = IntervalArithmetic.intersect_interval(ws[a1], l)
+            ws[a2] = IntervalArithmetic.intersect_interval(ws[a2], r)
+        elseif op === :sin;   _, x = sin_rev(z, ws[a1]);   ws[a1] = IntervalArithmetic.intersect_interval(ws[a1], x)
+        elseif op === :cos;   _, x = cos_rev(z, ws[a1]);   ws[a1] = IntervalArithmetic.intersect_interval(ws[a1], x)
+        elseif op === :tan;   _, x = tan_rev(z, ws[a1]);   ws[a1] = IntervalArithmetic.intersect_interval(ws[a1], x)
+        elseif op === :asin;  _, x = asin_rev(z, ws[a1]);  ws[a1] = IntervalArithmetic.intersect_interval(ws[a1], x)
+        elseif op === :acos;  _, x = acos_rev(z, ws[a1]);  ws[a1] = IntervalArithmetic.intersect_interval(ws[a1], x)
+        elseif op === :atan;  _, x = atan_rev(z, ws[a1]);  ws[a1] = IntervalArithmetic.intersect_interval(ws[a1], x)
+        elseif op === :sinh;  _, x = sinh_rev(z, ws[a1]);  ws[a1] = IntervalArithmetic.intersect_interval(ws[a1], x)
+        elseif op === :cosh;  _, x = cosh_rev(z, ws[a1]);  ws[a1] = IntervalArithmetic.intersect_interval(ws[a1], x)
+        elseif op === :tanh;  _, x = tanh_rev(z, ws[a1]);  ws[a1] = IntervalArithmetic.intersect_interval(ws[a1], x)
+        elseif op === :asinh; _, x = asinh_rev(z, ws[a1]); ws[a1] = IntervalArithmetic.intersect_interval(ws[a1], x)
+        elseif op === :acosh; _, x = acosh_rev(z, ws[a1]); ws[a1] = IntervalArithmetic.intersect_interval(ws[a1], x)
+        elseif op === :atanh; _, x = atanh_rev(z, ws[a1]); ws[a1] = IntervalArithmetic.intersect_interval(ws[a1], x)
+        elseif op === :exp;   _, x = exp_rev(z, ws[a1]);   ws[a1] = IntervalArithmetic.intersect_interval(ws[a1], x)
+        elseif op === :exp2;  _, x = exp2_rev(z, ws[a1]);  ws[a1] = IntervalArithmetic.intersect_interval(ws[a1], x)
+        elseif op === :exp10; _, x = exp10_rev(z, ws[a1]); ws[a1] = IntervalArithmetic.intersect_interval(ws[a1], x)
+        elseif op === :expm1; _, x = expm1_rev(z, ws[a1]); ws[a1] = IntervalArithmetic.intersect_interval(ws[a1], x)
+        elseif op === :log;   _, x = log_rev(z, ws[a1]);   ws[a1] = IntervalArithmetic.intersect_interval(ws[a1], x)
+        elseif op === :log2;  _, x = log2_rev(z, ws[a1]);  ws[a1] = IntervalArithmetic.intersect_interval(ws[a1], x)
+        elseif op === :log10; _, x = log10_rev(z, ws[a1]); ws[a1] = IntervalArithmetic.intersect_interval(ws[a1], x)
+        elseif op === :log1p; _, x = log1p_rev(z, ws[a1]); ws[a1] = IntervalArithmetic.intersect_interval(ws[a1], x)
+        elseif op === :sqrt;  _, x = sqrt_rev(z, ws[a1]);  ws[a1] = IntervalArithmetic.intersect_interval(ws[a1], x)
+        elseif op === :abs;   _, x = abs_rev(z, ws[a1]);   ws[a1] = IntervalArithmetic.intersect_interval(ws[a1], x)
+        elseif op === :sign;  _, x = sign_rev(z, ws[a1]);  ws[a1] = IntervalArithmetic.intersect_interval(ws[a1], x)
+        end
+    end
+end
+
+# --- Compiled propagation ---
+
+function propagate!(cd::CompiledDAG, X::IntervalBox; max_iter::Int=10, tol::Float64=1e-10)
+    ws = cd.workspace
+    nv = cd.n_vars
+    for _ in 1:max_iter
+        forward!(cd, X)
+        backward!(cd)
+        X_new = IntervalBox(ntuple(i -> @inbounds(ws[i]), nv))
+        isempty(X_new) && return X_new
+        max_change = maximum(diam(X[i]) - diam(X_new[i]) for i in 1:nv)
+        X = X_new
+        max_change < tol && break
+    end
+    return X
+end
+
+function propagate!(cd::CompiledDAG, X::IntervalBox, constraint::Interval;
+                    max_iter::Int=10, tol::Float64=1e-10)
+    ws = cd.workspace
+    nv = cd.n_vars
+    for _ in 1:max_iter
+        forward!(cd, X)
+        backward!(cd, constraint)
+        X_new = IntervalBox(ntuple(i -> @inbounds(ws[i]), nv))
+        isempty(X_new) && return X_new
+        max_change = maximum(diam(X[i]) - diam(X_new[i]) for i in 1:nv)
+        X = X_new
+        max_change < tol && break
+    end
+    return X
+end
diff --git a/src/dag/contractor.jl b/src/dag/contractor.jl
new file mode 100644
index 0000000..a660f43
--- /dev/null
+++ b/src/dag/contractor.jl
@@ -0,0 +1,138 @@
+"""
+    DAGContractor <: AbstractContractor
+
+A contractor based on compiled DAG propagation (HC4Revise).
+The DAG is built once, compiled into a flat instruction schedule, and
+reused across all calls — no dynamic dispatch, no per-node allocation.
+
+Constructor forms:
+- `DAGContractor(expr, vars)` — builds and compiles its own DAG
+- `DAGContractor(dag, expr, vars)` — adds to an existing SharedDAG and compiles
+"""
+struct DAGContractor{V} <: AbstractContractor
+    vars::V
+    dag::SharedDAG
+    compiled::CompiledDAG
+    ex::Any
+end
+
+function DAGContractor(ex::Num, vars::Vector{Num})
+    dag = SharedDAG(vars)
+    add_expression!(dag, ex, interval(0.0))
+    cd = compile(dag)
+    return DAGContractor(vars, dag, cd, ex)
+end
+
+function DAGContractor(dag::SharedDAG, ex::Num, vars::Vector{Num})
+    add_expression!(dag, ex, interval(0.0))
+    cd = compile(dag)
+    return DAGContractor(vars, dag, cd, ex)
+end
+
+function (C::DAGContractor)(X::IntervalBox, constraint::Interval=interval(0.0))
+    return propagate!(C.compiled, X, constraint)
+end
+
+"""
+    DAGSeparator <: AbstractSeparator
+
+A separator that uses compiled DAG-based propagation.
+Drop-in replacement for `Separator`.
+
+Constructor forms:
+- `DAGSeparator(expr, vars, constraint)` — builds and compiles its own DAG
+- `DAGSeparator(dag, expr, vars, constraint)` — uses an existing SharedDAG
+- `DAGSeparator(symbolic_inequality, vars)` — parses the inequality
+"""
+struct DAGSeparator{V, E, C, F} <: AbstractSeparator
+    vars::V
+    ex::E
+    constraint::C
+    f::F
+    contractor::DAGContractor
+end
+
+function DAGSeparator(expr::Num, vars::Vector{Num}, constraint::Interval)
+    contractor = DAGContractor(expr, vars)
+    f = make_function(expr, vars)
+    return DAGSeparator(vars, expr, constraint, f, contractor)
+end
+
+function DAGSeparator(dag::SharedDAG, expr::Num, vars::Vector{Num}, constraint::Interval)
+    contractor = DAGContractor(dag, expr, vars)
+    f = make_function(expr, vars)
+    return DAGSeparator(vars, expr, constraint, f, contractor)
+end
+
+function (SS::DAGSeparator)(X)
+    if any(x -> isinf(diam(x)), X)
+        return _separate_infinite_dag(SS, X)
+    end
+
+    boundary = SS.contractor(X)
+
+    lb = IntervalBox(inf.(X))
+    ub = IntervalBox(sup.(X))
+    inner = boundary
+    outer = boundary
+
+    lb_image = SS.f(lb)
+    if !isempty_interval(lb_image) && issubset_interval(lb_image, SS.constraint)
+        inner = inner ⊔ lb
+    else
+        outer = outer ⊔ lb
+    end
+
+    ub_image = SS.f(ub)
+    if !isempty_interval(ub_image) && issubset_interval(ub_image, SS.constraint)
+        inner = inner ⊔ ub
+    else
+        outer = outer ⊔ ub
+    end
+
+    return boundary, inner, outer
+end
+
+function _separate_infinite_dag(SS::DAGSeparator, X::IntervalBox)
+    C = SS.contractor
+    a, b = inf(SS.constraint), sup(SS.constraint)
+    inner = C(X, interval(a, b))
+
+    local outer
+    if a == -Inf
+        outer = C(X, interval(b, Inf))
+    elseif b == Inf
+        outer = C(X, interval(-Inf, a))
+    else
+        outer = C(X, interval(-Inf, a)) ⊔ C(X, interval(b, Inf))
+    end
+
+    boundary = inner ⊓ outer
+    return (boundary, inner, outer)
+end
+
+"""
+    DAGSeparator(orig_expr, vars)
+
+Build a DAG-based separator from a symbolic constraint expression,
+handling boolean combinations (`&`, `|`, `!`) by composing separators.
+"""
+function DAGSeparator(orig_expr, vars)
+    ex_val = Symbolics.value(orig_expr)
+
+    if SU.iscall(ex_val)
+        op = Symbolics.operation(ex_val)
+        args = Symbolics.arguments(ex_val)
+
+        if op === (&)
+            return DAGSeparator(Num(args[1]), vars) ⊓ DAGSeparator(Num(args[2]), vars)
+        elseif op === (|)
+            return DAGSeparator(Num(args[1]), vars) ⊔ DAGSeparator(Num(args[2]), vars)
+        elseif op === (!)
+            return ¬(DAGSeparator(Num(args[1]), vars))
+        end
+    end
+
+    ex, constraint = analyse(orig_expr)
+    return DAGSeparator(ex, vars, constraint)
+end
diff --git a/src/dag/nodes.jl b/src/dag/nodes.jl
new file mode 100644
index 0000000..e399d5f
--- /dev/null
+++ b/src/dag/nodes.jl
@@ -0,0 +1,94 @@
+abstract type AbstractNode end
+
+mutable struct VariableNode <: AbstractNode
+    name::Symbol
+    index::Int          # position in the input IntervalBox
+    range::Interval     # current interval bound
+end
+
+mutable struct ConstantNode <: AbstractNode
+    value::Float64
+    range::Interval     # [value, value], set once
+end
+
+ConstantNode(v::Real) = ConstantNode(Float64(v), interval(Float64(v)))
+
+mutable struct OperationNode <: AbstractNode
+    op::Function            # +, *, sin, ^, etc.
+    children::Vector{AbstractNode}
+    range::Interval         # current interval bound (set during forward pass)
+end
+
+OperationNode(op, children) = OperationNode(op, children, entireinterval())
+
+is_variable(n::VariableNode) = true
+is_variable(::AbstractNode) = false
+
+is_constant(n::ConstantNode) = true
+is_constant(::AbstractNode) = false
+
+is_operation(n::OperationNode) = true
+is_operation(::AbstractNode) = false
+
+function Base.show(io::IO, n::VariableNode)
+    print(io, "Var($(n.name), $(n.range))")
+end
+
+function Base.show(io::IO, n::ConstantNode)
+    print(io, "Const($(n.value))")
+end
+
+function Base.show(io::IO, n::OperationNode)
+    print(io, "Op($(n.op), $(length(n.children)) children, $(n.range))")
+end
+
+"""
+A root-constraint pair: an expression root node together with the interval
+constraint that its output must satisfy.
+"""
+struct ConstraintEntry
+    root::AbstractNode
+    constraint::Interval
+end
+
+"""
+    SharedDAG
+
+A persistent DAG that can hold multiple expressions sharing variable nodes
+and common subexpressions. Each expression is registered as a `ConstraintEntry`
+(root node + constraint interval).
+
+Built incrementally via `add_expression!`. The `nodes` vector is kept in
+topological order (leaves first), so forward propagation is a single pass
+and backward propagation is the reverse.
+"""
+mutable struct SharedDAG
+    variables::Vector{VariableNode}
+    nodes::Vector{AbstractNode}         # topological order
+    constraints::Vector{ConstraintEntry}
+    var_map::Dict{Symbol, VariableNode}
+    node_cache::Dict{UInt, AbstractNode}
+end
+
+function SharedDAG(variables::Vector{Num})
+    unwrapped = Symbolics.value.(variables)
+    var_nodes = VariableNode[]
+    var_map = Dict{Symbol, VariableNode}()
+    all_nodes = AbstractNode[]
+
+    for (i, v) in enumerate(unwrapped)
+        name = Symbolics.nameof(v)
+        vn = VariableNode(name, i, entireinterval())
+        push!(var_nodes, vn)
+        push!(all_nodes, vn)
+        var_map[name] = vn
+    end
+
+    return SharedDAG(var_nodes, all_nodes, ConstraintEntry[], var_map, Dict{UInt, AbstractNode}())
+end
+
+function Base.show(io::IO, dag::SharedDAG)
+    vars = join([string(v.name) for v in dag.variables], ", ")
+    nc = length(dag.constraints)
+    print(io, "SharedDAG($(length(dag.nodes)) nodes, vars=($vars), $nc constraints)")
+end
diff --git a/src/dag/propagate.jl b/src/dag/propagate.jl
new file mode 100644
index 0000000..9e66021
--- /dev/null
+++ b/src/dag/propagate.jl
@@ -0,0 +1,231 @@
+"""
+    forward!(dag::SharedDAG, X::IntervalBox)
+
+Forward propagation: set variable ranges from X, then evaluate all operation
+nodes in topological order.
+"""
+function forward!(dag::SharedDAG, X::IntervalBox)
+    for v in dag.variables
+        v.range = X[v.index]
+    end
+    for node in dag.nodes
+        if node isa OperationNode
+            node.range = _eval_forward(node.op, node.children)
+        end
+    end
+end
+
+function _eval_forward(op::Function, children::Vector{AbstractNode})
+    if length(children) == 1
+        return op(children[1].range)
+    elseif length(children) == 2
+        return op(children[1].range, children[2].range)
+    else
+        return op([c.range for c in children]...)
+    end
+end
+
+function _eval_forward(::typeof(^), children::Vector{AbstractNode})
+    base = children[1].range
+    exp_node = children[2]
+    if exp_node isa ConstantNode
+        n = exp_node.value
+        if n == round(Int, n)
+            return base ^ round(Int, n)
+        end
+    end
+    return base ^ exp_node.range
+end
+
+"""
+    backward!(dag::SharedDAG)
+
+Backward propagation for all constraints. For each constraint entry, intersect
+the root's range with the constraint, then propagate backward through the DAG.
+
+When multiple constraints share nodes, each backward pass further narrows
+shared variable and subexpression ranges.
+"""
+function backward!(dag::SharedDAG)
+    for entry in dag.constraints
+        entry.root.range = IntervalArithmetic.intersect_interval(entry.root.range, entry.constraint)
+    end
+
+    all_empty = all(isempty_interval(e.root.range) for e in dag.constraints)
+    if all_empty
+        for v in dag.variables
+            v.range = emptyinterval()
+        end
+        return
+    end
+
+    for i in length(dag.nodes):-1:1
+        node = dag.nodes[i]
+        node isa OperationNode || continue
+        _contract_backward!(node)
+    end
+end
+
+"""
+    backward!(dag::SharedDAG, constraint::Interval)
+
+Backward propagation with a single override constraint applied to the first
+(or only) root. Used by DAGContractor when contracting with a specific constraint.
+"""
+function backward!(dag::SharedDAG, constraint::Interval)
+    entry = dag.constraints[1]
+    entry.root.range = IntervalArithmetic.intersect_interval(entry.root.range, constraint)
+
+    if isempty_interval(entry.root.range)
+        for v in dag.variables
+            v.range = emptyinterval()
+        end
+        return
+    end
+
+    for i in length(dag.nodes):-1:1
+        node = dag.nodes[i]
+        node isa OperationNode || continue
+        _contract_backward!(node)
+    end
+end
+
+function _contract_backward!(node::OperationNode)
+    z = node.range
+    children = node.children
+    if isempty_interval(z)
+        for c in children
+            c.range = emptyinterval()
+        end
+        return
+    end
+    if length(children) == 1
+        _contract_unary!(node.op, z, children[1])
+    elseif length(children) == 2
+        _contract_binary!(node.op, z, children[1], children[2])
+    end
+end
+
+# --- Binary reverse contractors ---
+
+function _contract_binary!(::typeof(+), z, left, right)
+    _, l_new, r_new = plus_rev(z, left.range, right.range)
+    left.range = IntervalArithmetic.intersect_interval(left.range, l_new)
+    right.range = IntervalArithmetic.intersect_interval(right.range, r_new)
+end
+
+function _contract_binary!(::typeof(-), z, left, right)
+    _, l_new, r_new = minus_rev(z, left.range, right.range)
+    left.range = IntervalArithmetic.intersect_interval(left.range, l_new)
+    right.range = IntervalArithmetic.intersect_interval(right.range, r_new)
+end
+
+function _contract_binary!(::typeof(*), z, left, right)
+    _, l_new, r_new = mul_rev(z, left.range, right.range)
+    left.range = IntervalArithmetic.intersect_interval(left.range, l_new)
+    right.range = IntervalArithmetic.intersect_interval(right.range, r_new)
+end
+
+function _contract_binary!(::typeof(/), z, left, right)
+    _, l_new, r_new = div_rev(z, left.range, right.range)
+    left.range = IntervalArithmetic.intersect_interval(left.range, l_new)
+    right.range = IntervalArithmetic.intersect_interval(right.range, r_new)
+end
+
+function _contract_binary!(::typeof(^), z, base_node, exp_node)
+    if exp_node isa ConstantNode
+        n = exp_node.value
+        if n == round(Int, n)
+            _, x_new, _ = power_rev(z, base_node.range, round(Int, n))
+            base_node.range = IntervalArithmetic.intersect_interval(base_node.range, x_new)
+            return
+        end
+    end
+    _, x_new, y_new = power_rev(z, base_node.range, exp_node.range)
+    base_node.range = IntervalArithmetic.intersect_interval(base_node.range, x_new)
+    exp_node.range = IntervalArithmetic.intersect_interval(exp_node.range, y_new)
+end
+
+function _contract_binary!(::typeof(min), z, left, right)
+    _, l_new, r_new = min_rev(z, left.range, right.range)
+    left.range = IntervalArithmetic.intersect_interval(left.range, l_new)
+    right.range = IntervalArithmetic.intersect_interval(right.range, r_new)
+end
+
+function _contract_binary!(::typeof(max), z, left, right)
+    _, l_new, r_new = max_rev(z, left.range, right.range)
+    left.range = IntervalArithmetic.intersect_interval(left.range, l_new)
+    right.range = IntervalArithmetic.intersect_interval(right.range, r_new)
+end
+
+# --- Unary reverse contractors ---
+
+for (jl_op, rev_fn) in [
+    (:sin, :sin_rev), (:cos, :cos_rev), (:tan, :tan_rev),
+    (:asin, :asin_rev), (:acos, :acos_rev), (:atan, :atan_rev),
+    (:sinh, :sinh_rev), (:cosh, :cosh_rev), (:tanh, :tanh_rev),
+    (:asinh, :asinh_rev), (:acosh, :acosh_rev), (:atanh, :atanh_rev),
+    (:exp, :exp_rev), (:exp2, :exp2_rev), (:exp10, :exp10_rev),
+    (:expm1, :expm1_rev),
+    (:log, :log_rev), (:log2, :log2_rev), (:log10, :log10_rev),
+    (:log1p, :log1p_rev),
+    (:sqrt, :sqrt_rev), (:abs, :abs_rev), (:sign, :sign_rev),
+]
+    @eval function _contract_unary!(::typeof($jl_op), z, child)
+        _, x_new = $rev_fn(z, child.range)
+        child.range = IntervalArithmetic.intersect_interval(child.range, x_new)
+    end
+end
+
+function _contract_unary!(::typeof(-), z, child)
+    _, x_new = minus_rev(z, child.range)
+    child.range = IntervalArithmetic.intersect_interval(child.range, x_new)
+end
+
+function _contract_unary!(op, z, child)
+    @warn "No reverse contractor for unary $op — skipping" maxlog=1
+end
+
+function _contract_binary!(op, z, left, right)
+    @warn "No reverse contractor for binary $op — skipping" maxlog=1
+end
+
+"""
+    propagate!(dag::SharedDAG, X::IntervalBox; max_iter=10, tol=1e-10) → IntervalBox
+
+Iterate forward-backward propagation on the DAG until convergence.
+All constraints in the DAG are propagated each iteration, so narrowing
+from one constraint benefits the others via shared variable nodes.
+"""
+function propagate!(dag::SharedDAG, X::IntervalBox; max_iter::Int=10, tol::Float64=1e-10)
+    for _ in 1:max_iter
+        old_ranges = IntervalBox(Tuple(v.range for v in dag.variables))
+        forward!(dag, X)
+        backward!(dag)
+        X = IntervalBox(Tuple(v.range for v in dag.variables))
+        isempty(X) && return X
+        max_change = maximum(diam(old) - diam(new) for (old, new) in zip(old_ranges, X))
+        max_change < tol && break
+    end
+    return X
+end
+
+"""
+    propagate!(dag::SharedDAG, X::IntervalBox, constraint::Interval; ...) → IntervalBox
+
+Propagate with a single override constraint (ignoring stored constraints).
+Used by DAGContractor for boundary/complement contraction.
+"""
+function propagate!(dag::SharedDAG, X::IntervalBox, constraint::Interval;
+                    max_iter::Int=10, tol::Float64=1e-10)
+    for _ in 1:max_iter
+        old_ranges = IntervalBox(Tuple(v.range for v in dag.variables))
+        forward!(dag, X)
+        backward!(dag, constraint)
+        X = IntervalBox(Tuple(v.range for v in dag.variables))
+        isempty(X) && return X
+        max_change = maximum(diam(old) - diam(new) for (old, new) in zip(old_ranges, X))
+        max_change < tol && break
+    end
+    return X
+end
diff --git a/src/utils.jl b/src/utils.jl
index abff343..64b76cd 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -22,19 +22,19 @@ function analyse(ex)
 	if op ∈ (≤, <)
 		constraint = interval(-Inf, 0)
 		Num(lhs - rhs), constraint
-		
+
 	elseif op ∈ (≥, >)
 		constraint = interval(0, +Inf)
 		Num(lhs - rhs), constraint
-		
+
 	elseif op == (==)
 		constraint = interval(0, 0)
 		Num(lhs - rhs), constraint
-	
+
 	else
 		return ex, interval(0, 0)   # implicit 0
 	end
-		
+
 end
 
 
@@ -59,23 +59,23 @@ function separator(ex, vars)
 	lhs, rhs = arguments(ex2)
 
 	if op == &
-		return separator(lhs, vars) ∩ separator(rhs, vars)
+		return separator(lhs, vars) ⊓ separator(rhs, vars)
 
 	elseif op == |
-		return separator(lhs, vars) ∪ separator(rhs, vars)
-	
+		return separator(lhs, vars) ⊔ separator(rhs, vars)
+
 	elseif op ∈ (≤, <)
 		constraint = interval(-Inf, 0)
 		Separator(Num(lhs - rhs), vars, constraint)
-		
+
 	elseif op ∈ (≥, >)
 		constraint = interval(0, +Inf)
 		Separator(Num(lhs - rhs), vars, constraint)
-		
+
 	elseif op == (==)
 		constraint = interval(0, 0)
 		Separator(Num(lhs - rhs), vars, constraint)
-	
+
 	else
 		return Separator(ex, vars, interval(0, 0))   # implicit "== 0"
 	end
diff --git a/test/test_dag.jl b/test/test_dag.jl
new file mode 100644
index 0000000..2ca8c48
--- /dev/null
+++ b/test/test_dag.jl
@@ -0,0 +1,181 @@
+using IntervalArithmetic, IntervalArithmetic.Symbols
+using IntervalConstraintProgramming
+using IntervalConstraintProgramming: build_dag, forward!, backward!, propagate!,
+    SharedDAG, add_expression!, ConstraintEntry
+using Symbolics
+using IntervalBoxes
+using Test
+
+const IType{T} = Union{Interval{T}, BareInterval{T}}
+eq(a::IType, b::IType) = isequal_interval(bareinterval(a), bareinterval(b))
+eq(a::IntervalBox, b::IntervalBox) = all(eq.(a, b))
+eq(a::Vector, b::Vector) = length(a) == length(b) && all(eq.(a, b))
+
+@testset "SharedDAG construction" begin
+    vars = @variables x, y
+
+    dag = SharedDAG(vars)
+    @test length(dag.variables) == 2
+    @test length(dag.nodes) == 2  # just variable nodes
+    @test dag.variables[1].name == :x
+    @test dag.variables[2].name == :y
+
+    add_expression!(dag, x^2 + y^2, interval(-Inf, 0.0))
+    @test length(dag.constraints) == 1
+    @test length(dag.nodes) > 2
+end
+
+@testset "SharedDAG node reuse (CSE)" begin
+    vars = @variables x, y
+
+    dag = SharedDAG(vars)
+    add_expression!(dag, x^2 + y^2 - 1, interval(-Inf, 0.0))
+    add_expression!(dag, x^2 - y, interval(-Inf, 0.0))
+
+    pow_nodes = filter(n -> n isa IntervalConstraintProgramming.OperationNode && n.op === (^), dag.nodes)
+    @test length(pow_nodes) == 2  # x^2 and y^2, not a duplicate x^2
+
+    @test length(dag.constraints) == 2
+end
+
+@testset "SharedDAG n-ary decomposition" begin
+    vars = @variables x, y
+
+    dag = build_dag(x + y + 1, vars, interval(0.0))
+    ops = filter(n -> n isa IntervalConstraintProgramming.OperationNode, dag.nodes)
+    for op in ops
+        @test length(op.children) <= 2
+    end
+end
+
+@testset "Forward propagation" begin
+    vars = @variables x, y
+
+    dag = build_dag(x^2 + y^2, vars, interval(-Inf, 1.0))
+    X = IntervalBox(interval(0.0, 1.0), interval(0.0, 1.0))
+    forward!(dag, X)
+    root_range = dag.constraints[1].root.range
+    @test inf(root_range) <= 0.0 + 1e-10
+    @test sup(root_range) >= 2.0 - 1e-10
+end
+
+@testset "Single-expression propagation" begin
+    vars = @variables x, y
+
+    ex, constraint = IntervalConstraintProgramming.analyse(x^2 + y^2 <= 1)
+    dag = build_dag(ex, vars, constraint)
+    X = IntervalBox(interval(-10.0, 10.0), 2)
+    result = propagate!(dag, X)
+    @test all(inf.(result) .>= -1.0 - 1e-10)
+    @test all(sup.(result) .<= 1.0 + 1e-10)
+end
+
+@testset "Multi-constraint propagation" begin
+    vars = @variables x, y
+
+    dag = SharedDAG(vars)
+    add_expression!(dag, x^2 + y^2 - 1, interval(-Inf, 0.0))
+    add_expression!(dag, x - y, interval(0.0, Inf))
+
+    X = IntervalBox(interval(-10, 10), 2)
+    result = propagate!(dag, X)
+    @test all(inf.(result) .>= -1.0 - 1e-10)
+    @test all(sup.(result) .<= 1.0 + 1e-10)
+
+    # y ≥ x² narrows y ≥ 0
+    dag2 = SharedDAG(vars)
+    add_expression!(dag2, x^2 + y^2 - 1, interval(-Inf, 0.0))
+    add_expression!(dag2, x^2 - y, interval(-Inf, 0.0))
+    result2 = propagate!(dag2, IntervalBox(interval(-10, 10), 2))
+    @test inf(result2[2]) >= -1e-10  # y ≥ 0
+end
+
+@testset "DAGContractor" begin
+    vars = @variables x, y
+
+    C = DAGContractor(x^2 + y^2, vars)
+    X = IntervalBox(-Inf..Inf, 2)
+    result = C(X, -Inf..1)
+    @test eq(result, IntervalBox(-1..1, 2))
+
+    vars3 = @variables x, y, z
+    X3 = IntervalBox(interval(0.5, 1.5), 3)
+    C1 = DAGContractor(x + y, vars3)
+    @test eq(C1(X3, -Inf..1), IntervalBox(interval(0.5, 0.5), interval(0.5, 0.5), interval(0.5, 1.5)))
+end
+
+@testset "DAGContractor with shared DAG" begin
+    vars = @variables x, y
+
+    dag = SharedDAG(vars)
+    C1 = DAGContractor(dag, x^2 + y^2, vars)
+    C2 = DAGContractor(dag, x - y, vars)
+
+    @test C1.dag === C2.dag
+    @test length(dag.constraints) == 2
+end
+
+@testset "DAGSeparator" begin
+    vars = @variables x, y
+
+    II = interval(-100, 100)
+    X = IntervalBox(II, 2)
+    S = DAGSeparator(x^2 + y^2 <= 1, vars)
+
+    boundary, inner, outer = S(X)
+    @test eq(inner, IntervalBox(-1..1, 2))
+    @test eq(outer, IntervalBox(II, 2))
+
+    X_inf = IntervalBox(-Inf..Inf, -Inf..Inf)
+    boundary2, inner2, outer2 = S(X_inf)
+    @test eq(inner2, IntervalBox(-1..1, 2))
+    @test eq(outer2, IntervalBox(-Inf..Inf, 2))
+end
+
+@testset "DAGSeparator with shared DAG" begin
+    vars = @variables x, y
+
+    dag = SharedDAG(vars)
+    ex1, c1 = IntervalConstraintProgramming.analyse(x^2 + y^2 <= 1)
+    ex2, c2 = IntervalConstraintProgramming.analyse(x >= 0)
+    S1 = DAGSeparator(dag, ex1, vars, c1)
+    S2 = DAGSeparator(dag, ex2, vars, c2)
+
+    @test S1.contractor.dag === S2.contractor.dag
+    @test length(dag.constraints) == 2
+end
+
+@testset "DAG pave matches original" begin
+    vars = @variables x, y
+
+    S_orig = Separator(x^2 + y^2 <= 1, vars)
+    S_dag = DAGSeparator(x^2 + y^2 <= 1, vars)
+
+    X = IntervalBox(interval(-3, 3), 2)
+    inner_orig, boundary_orig = pave(X, S_orig, 0.5)
+    inner_dag, boundary_dag = pave(X, S_dag, 0.5)
+
+    @test length(inner_orig) == length(inner_dag)
+    @test length(boundary_orig) == length(boundary_dag)
+end
+
+@testset "DAG with combined separators" begin
+    vars = @variables x, y
+
+    S1 = DAGSeparator(x > 0, vars)
+    S2 = DAGSeparator(y > 0, vars)
+    S = S1 ⊓ S2
+
+    X = IntervalBox(interval(-3, 3), 2)
+    boundary, inner, outer = S(X)
+    @test eq(inner, IntervalBox(0..3, 2))
+end
+
+@testset "DAG 3D torus" begin
+    vars = @variables x, y, z
+
+    S = DAGSeparator(3 - sqrt(x^2 + y^2)^2 + z^2 <= 1, vars)
+    X = IntervalBox(interval(-10, 10), 3)
+    inner, boundary = pave(X, S, 1.0)
+    @test inner[1] isa IntervalBox{3, Float64, Interval{Float64}}
+end