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/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/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)