diff --git a/KnotDiagrams/Project.toml b/KnotDiagrams/Project.toml
new file mode 100644
index 0000000..a2baf9e
--- /dev/null
+++ b/KnotDiagrams/Project.toml
@@ -0,0 +1,17 @@
+name = "KnotDiagrams"
+uuid = "a7e3c1f0-5d4b-4e8a-9f1c-2b3d4e5f6a7b"
+version = "0.1.0"
+
+[deps]
+DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
+
+[extras]
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
+
+[targets]
+test = ["Test"]
+
+[compat]
+DataFrames = "1"
+julia = "1.9"
diff --git a/KnotDiagrams/docs/narrative.txt b/KnotDiagrams/docs/narrative.txt
new file mode 100644
index 0000000..b3bfec9
--- /dev/null
+++ b/KnotDiagrams/docs/narrative.txt
@@ -0,0 +1,105 @@
+Fraying-Model Fragment-Count Distribution Predictions
+=====================================================
+
+This document states testable predictions about the fragment-count
+distribution L_D(k) for knot diagrams, derived from the combinatorics
+of Kauffman bracket state sums.
+
+Definitions
+-----------
+
+For an n-crossing knot diagram D, the 2^n global resolution states are
+obtained by choosing one of two smoothings at each crossing.  Each state
+yields a collection of k disjoint simple closed curves ("fragments").
+
+  L_D(k) = |{ states s : s produces exactly k fragments }|
+
+The tropical summary is:
+  support(L_D) = { k : L_D(k) > 0 }
+  min, max     = extrema of support
+  |support|    = cardinality
+  class        = "spread" if |support| > 1, "collapse" if |support| = 1
+
+Prediction 1 — No Collapse for Nontrivial Knots
+------------------------------------------------
+
+Every prime knot with crossing number >= 3 has |support| > 1 (is "spread").
+
+Rationale: For an alternating prime knot, the all-A state (Seifert state)
+produces a number of circles related to the genus, while the all-B state
+produces a different count.  The minimum and maximum circle counts in
+the state sum always differ for nontrivial knots.
+
+Prediction 2 — Trefoil (3_1) Has Support {1, 2, 3}
+---------------------------------------------------
+
+The trefoil should produce states with 1, 2, or 3 circles.
+
+  L_{3_1}(1) = 3    (three states yield 1 circle)
+  L_{3_1}(2) = 4    (four states yield 2 circles)
+  L_{3_1}(3) = 1    (one state yields 3 circles)
+
+Verification: the all-0 state (smoothing 0 at every crossing) gives 2
+circles, and the all-1 state gives 3 circles.  Three mixed states give
+1 circle each.  This matches the Kauffman bracket structure for the
+trefoil.
+
+Prediction 3 — Topology Predicts the Distribution
+--------------------------------------------------
+
+Diagrams of different knot types at the same crossing number should
+produce different fragment-count histograms L_D(k).  Specifically:
+
+  5_1  vs  5_2   — different histograms, likely different supports
+  6_1  vs  6_2   — different histograms
+  6_1  vs  6_3   — different histograms
+  6_2  vs  6_3   — different histograms
+  7_i  vs  7_j   — different histograms for i ≠ j (all 21 pairs)
+
+This is the central testable claim: the fragment-count distribution is a
+topological invariant of the diagram (up to diagram equivalence), and
+distinct knot types generally have distinct distributions.
+
+Prediction 4 — Torus Knots Have Wider Support
+----------------------------------------------
+
+Torus knots T(2, 2k+1) have wider support (larger |support|) than
+non-torus knots at the same crossing number.
+
+  5_1 (torus)     — wider support than 5_2 (twist knot)
+  7_1 (torus)     — wider support than 7_2, ..., 7_7
+
+Rationale: the Seifert state of a torus knot produces the minimum
+number of Seifert circles (k+1 for T(2,2k+1)), while the all-opposite
+state produces a larger count.  The spread between min and max is
+maximised for torus knots among alternating knots at a given crossing
+number.
+
+Prediction 5 — Support Range Grows with Crossing Number
+--------------------------------------------------------
+
+For the torus knot series T(2, 2k+1):
+  3_1: |support| around 3
+  5_1: |support| around 4-5
+  7_1: |support| around 5-6
+
+The support range (max - min + 1) grows roughly linearly with crossing
+number for torus knots.
+
+Prediction 6 — The Distribution Refines the Jones Polynomial
+-------------------------------------------------------------
+
+Two knots can share the same Jones polynomial but have different
+fragment-count distributions.  The L_D(k) histogram retains state-level
+information that the polynomial (which sums over states with signs and
+weights) discards.  This makes L_D(k) a finer invariant in some cases.
+
+Expected Outcomes
+-----------------
+
+Running the analysis on all prime knots up to 7 crossings should show:
+  - 14 diagrams, all classified as "spread"
+  - 25 pairwise comparisons (1 + 3 + 21 for n=5,6,7)
+  - Majority of pairs have different supports
+  - All pairs have different histograms (stronger claim)
+  - Torus knots (3_1, 5_1, 7_1) have the widest support at each cn
diff --git a/KnotDiagrams/scripts/plot_distributions.jl b/KnotDiagrams/scripts/plot_distributions.jl
new file mode 100644
index 0000000..bcfa8ae
--- /dev/null
+++ b/KnotDiagrams/scripts/plot_distributions.jl
@@ -0,0 +1,92 @@
+#!/usr/bin/env julia
+# SPDX-License-Identifier: PMPL-1.0-or-later
+#
+# Generate histogram plots for fragment-count distributions.
+# Requires CairoMakie.
+#
+# Usage:
+#   julia --project=. scripts/plot_distributions.jl
+
+using KnotDiagrams
+using CairoMakie
+
+function plot_all(fds::AbstractVector{FragmentDistribution};
+                  output_dir::String = joinpath(@__DIR__, "..", "plots"))
+    mkpath(output_dir)
+
+    # Group by crossing number
+    groups = Dict{Int,Vector{FragmentDistribution}}()
+    for fd in fds
+        cn = fd.diagram.crossing_number
+        push!(get!(Vector{FragmentDistribution}, groups, cn), fd)
+    end
+
+    # Per-crossing-number overlaid histograms
+    for cn in sort(collect(keys(groups)))
+        group = groups[cn]
+        plot_group(group, cn, output_dir)
+    end
+
+    # Per-diagram individual plots
+    for fd in fds
+        plot_single(fd, output_dir)
+    end
+
+    n_group = length(groups)
+    n_single = length(fds)
+    println("Saved $(n_group + n_single) plots to $output_dir")
+end
+
+function plot_group(group::Vector{FragmentDistribution}, cn::Int, output_dir::String)
+    all_k = Set{Int}()
+    for fd in group
+        union!(all_k, keys(fd.histogram))
+    end
+    k_range = minimum(all_k):maximum(all_k)
+
+    fig = Figure(size = (800, 500))
+    ax = Axis(fig[1, 1],
+        title = "Fragment-count distribution — $(cn)-crossing knots",
+        xlabel = "Fragment count k",
+        ylabel = "Number of states",
+        xticks = collect(k_range),
+    )
+
+    n_diagrams = length(group)
+    bar_width = 0.8 / n_diagrams
+    offsets = range(-0.4 + bar_width / 2, step = bar_width, length = n_diagrams)
+
+    for (i, fd) in enumerate(group)
+        ks = collect(k_range)
+        counts = [get(fd.histogram, k, 0) for k in ks]
+        barplot!(ax, ks .+ offsets[i], counts,
+            width = bar_width,
+            label = fd.diagram.name,
+            color = Cycled(i),
+        )
+    end
+
+    axislegend(ax, position = :rt)
+    save(joinpath(output_dir, "distribution_n$(cn).png"), fig, px_per_unit = 2)
+end
+
+function plot_single(fd::FragmentDistribution, output_dir::String)
+    k_range = fd.summary.min_k:fd.summary.max_k
+    ks = collect(k_range)
+    counts = [get(fd.histogram, k, 0) for k in ks]
+
+    fig = Figure(size = (600, 400))
+    ax = Axis(fig[1, 1],
+        title = "L_D(k) — $(fd.diagram.name)",
+        xlabel = "Fragment count k",
+        ylabel = "Number of states",
+        xticks = ks,
+    )
+    barplot!(ax, ks, counts, color = :steelblue)
+
+    save(joinpath(output_dir, "distribution_$(fd.diagram.name).png"), fig, px_per_unit = 2)
+end
+
+# Main
+results = analyse_all(verbose = false)
+plot_all(results.distributions)
diff --git a/KnotDiagrams/scripts/run_analysis.jl b/KnotDiagrams/scripts/run_analysis.jl
new file mode 100644
index 0000000..ff703e8
--- /dev/null
+++ b/KnotDiagrams/scripts/run_analysis.jl
@@ -0,0 +1,142 @@
+#!/usr/bin/env julia
+# SPDX-License-Identifier: PMPL-1.0-or-later
+#
+# Full analysis pipeline: compute fragment-count distributions for all
+# prime knots up to 7 crossings and compare within crossing numbers.
+#
+# Usage:
+#   julia --project=. scripts/run_analysis.jl
+#
+# For plots, run separately:
+#   julia --project=. scripts/plot_distributions.jl
+
+using KnotDiagrams
+using DataFrames
+
+function main()
+    println("=" ^ 72)
+    println("  Knot Diagram Fragment-Count Distribution Analysis")
+    println("=" ^ 72)
+    println()
+
+    # Run full analysis
+    results = analyse_all(verbose = true)
+    println()
+
+    # Print summary table
+    println("=" ^ 72)
+    println("SUMMARY TABLE")
+    println("=" ^ 72)
+    show(stdout, results.summary_table; allrows = true, allcols = true)
+    println("\n")
+
+    # Print comparison table
+    println("=" ^ 72)
+    println("COMPARISON TABLE")
+    println("=" ^ 72)
+    show(stdout, results.comparison_table; allrows = true, allcols = true)
+    println("\n")
+
+    # Validate against narrative predictions
+    validate_narrative(results)
+
+    return results
+end
+
+"""
+    validate_narrative(results)
+
+Check results against the narrative predictions in docs/narrative.txt.
+"""
+function validate_narrative(results)
+    println("=" ^ 72)
+    println("NARRATIVE VALIDATION")
+    println("=" ^ 72)
+
+    passed = 0
+    failed = 0
+
+    # Prediction 1: All prime knots are "spread"
+    print("  [1] All prime knots have |support| > 1 (spread): ")
+    all_spread = all(fd -> fd.summary.classification == :spread, results.distributions)
+    if all_spread
+        println("PASS")
+        passed += 1
+    else
+        collapsed = filter(fd -> fd.summary.classification == :collapse, results.distributions)
+        println("FAIL — collapsed: $(join([fd.diagram.name for fd in collapsed], ", "))")
+        failed += 1
+    end
+
+    # Prediction 2: Trefoil support includes at least {1,2} or {2,3}
+    print("  [2] Trefoil support contains multiple values: ")
+    trefoil = first(filter(fd -> fd.diagram.name == "3_1", results.distributions))
+    if trefoil.summary.support_size >= 2
+        println("PASS (support = $(trefoil.summary.support))")
+        passed += 1
+    else
+        println("FAIL (support = $(trefoil.summary.support))")
+        failed += 1
+    end
+
+    # Prediction 3: Figure-eight has different distribution from trefoil
+    print("  [3] Figure-eight distribution differs from trefoil: ")
+    fig8 = first(filter(fd -> fd.diagram.name == "4_1", results.distributions))
+    if fig8.histogram != trefoil.histogram
+        println("PASS")
+        passed += 1
+    else
+        println("FAIL")
+        failed += 1
+    end
+
+    # Prediction 4: 5_1 and 5_2 have different distributions
+    print("  [4] 5_1 vs 5_2 histograms differ: ")
+    cn5 = filter(r -> r.crossing_number == 5, results.comparisons)
+    if !isempty(cn5) && first(cn5).histograms_differ
+        println("PASS")
+        passed += 1
+    else
+        println("FAIL")
+        failed += 1
+    end
+
+    # Prediction 5: Most pairs at same crossing number have differing supports
+    print("  [5] Majority of same-cn pairs have differing supports: ")
+    n_total = length(results.comparisons)
+    n_differ = count(r -> r.supports_differ, results.comparisons)
+    ratio = n_total > 0 ? n_differ / n_total : 0.0
+    if ratio > 0.5
+        println("PASS ($n_differ / $n_total = $(round(ratio * 100, digits=1))%)")
+        passed += 1
+    else
+        println("FAIL ($n_differ / $n_total = $(round(ratio * 100, digits=1))%)")
+        failed += 1
+    end
+
+    # Prediction 6: Torus knots have wider support
+    print("  [6] Torus knots tend to have wider support: ")
+    fd51 = first(filter(fd -> fd.diagram.name == "5_1", results.distributions))
+    fd52 = first(filter(fd -> fd.diagram.name == "5_2", results.distributions))
+    fd71 = first(filter(fd -> fd.diagram.name == "7_1", results.distributions))
+    other_7 = filter(fd -> fd.diagram.crossing_number == 7 && fd.diagram.name != "7_1",
+                     results.distributions)
+    avg_other = isempty(other_7) ? 0.0 :
+        sum(fd.summary.support_size for fd in other_7) / length(other_7)
+
+    torus_wider = fd51.summary.support_size >= fd52.summary.support_size &&
+                  fd71.summary.support_size >= avg_other
+    if torus_wider
+        println("PASS (5_1: $(fd51.summary.support_size) vs 5_2: $(fd52.summary.support_size), " *
+                "7_1: $(fd71.summary.support_size) vs avg others: $(round(avg_other, digits=1)))")
+        passed += 1
+    else
+        println("INCONCLUSIVE (5_1: $(fd51.summary.support_size) vs 5_2: $(fd52.summary.support_size), " *
+                "7_1: $(fd71.summary.support_size) vs avg others: $(round(avg_other, digits=1)))")
+    end
+
+    println()
+    println("  Passed: $passed / $(passed + failed)")
+end
+
+main()
diff --git a/KnotDiagrams/src/KnotDiagrams.jl b/KnotDiagrams/src/KnotDiagrams.jl
new file mode 100644
index 0000000..1bc19a2
--- /dev/null
+++ b/KnotDiagrams/src/KnotDiagrams.jl
@@ -0,0 +1,113 @@
+# SPDX-License-Identifier: PMPL-1.0-or-later
+"""
+    KnotDiagrams
+
+Compute fragment-count distributions L_D(k) for knot diagrams represented
+as signed Gauss codes.
+
+For an n-crossing diagram, the package enumerates all 2^n global resolution
+states (one binary smoothing choice per crossing), counts the resulting
+closed-curve fragments via union-find, and builds a histogram mapping
+fragment count k to the number of states producing k fragments.
+
+Tropical summaries (support, min, max, spread vs. collapse) and pairwise
+comparisons across diagrams with the same crossing number are provided to
+test whether topology predicts the fragment-count distribution.
+
+# Quick start
+```julia
+using KnotDiagrams
+
+# Analyse all prime knots up to 7 crossings
+results = analyse_all()
+
+# Single diagram
+trefoil = KnotDiagram("3_1", 3, GaussCode([1, -2, 3, -1, 2, -3]))
+fd = fragment_distribution(trefoil)
+println(fd)
+```
+"""
+module KnotDiagrams
+
+using DataFrames
+
+include("union_find.jl")
+include("types.jl")
+include("smoothing.jl")
+include("distribution.jl")
+include("comparison.jl")
+include("knot_table.jl")
+
+# --- Public API ---
+
+export GaussCode, KnotDiagram, Crossing
+export FragmentDistribution, TropicalSummary, ComparisonResult
+export UnionFind, find!, uf_merge!, num_components
+export count_circles, enumerate_states
+export fragment_distribution, tropical_summary
+export distribution_table
+export compare_pair, compare_by_crossing_number, comparison_table
+export prime_knots_up_to_7, torus_knot_gauss, dt_to_gauss
+export analyse_all
+
+"""
+    analyse_all(; verbose::Bool = true) -> NamedTuple
+
+Run the full analysis pipeline on all prime knots up to 7 crossings:
+1. Compute fragment-count distributions
+2. Build summary table
+3. Compare pairs within each crossing number
+4. Return all results
+
+Returns a `NamedTuple` with fields:
+- `diagrams`: the KnotDiagram vector
+- `distributions`: the FragmentDistribution vector
+- `summary_table`: DataFrame of per-diagram summaries
+- `comparisons`: vector of ComparisonResult
+- `comparison_table`: DataFrame of pairwise comparisons
+"""
+function analyse_all(; verbose::Bool = true)
+    diagrams = prime_knots_up_to_7()
+
+    verbose && println("Computing fragment-count distributions for $(length(diagrams)) diagrams...")
+    distributions = [fragment_distribution(d) for d in diagrams]
+
+    if verbose
+        println()
+        for fd in distributions
+            println(fd)
+            println()
+        end
+    end
+
+    summary_df = distribution_table(distributions)
+    comparisons = compare_by_crossing_number(distributions)
+    comparison_df = comparison_table(comparisons)
+
+    if verbose
+        println("=" ^ 72)
+        println("PAIRWISE COMPARISONS (same crossing number)")
+        println("=" ^ 72)
+        for r in comparisons
+            println("  ", r)
+        end
+        println()
+
+        n_pairs = length(comparisons)
+        n_differ_support = count(r -> r.supports_differ, comparisons)
+        n_differ_hist = count(r -> r.histograms_differ, comparisons)
+        println("Total pairs compared: $n_pairs")
+        println("Pairs with differing supports: $n_differ_support / $n_pairs")
+        println("Pairs with differing histograms: $n_differ_hist / $n_pairs")
+    end
+
+    return (
+        diagrams = diagrams,
+        distributions = distributions,
+        summary_table = summary_df,
+        comparisons = comparisons,
+        comparison_table = comparison_df,
+    )
+end
+
+end # module KnotDiagrams
diff --git a/KnotDiagrams/src/comparison.jl b/KnotDiagrams/src/comparison.jl
new file mode 100644
index 0000000..51db7c9
--- /dev/null
+++ b/KnotDiagrams/src/comparison.jl
@@ -0,0 +1,94 @@
+# SPDX-License-Identifier: PMPL-1.0-or-later
+#
+# Compare fragment-count distributions across diagrams.
+
+"""
+    ComparisonResult
+
+Result of comparing two knot diagrams with the same crossing number.
+"""
+struct ComparisonResult
+    diagram_a::String
+    diagram_b::String
+    crossing_number::Int
+    support_a::Vector{Int}
+    support_b::Vector{Int}
+    supports_differ::Bool
+    histograms_differ::Bool
+end
+
+"""
+    compare_pair(fd_a::FragmentDistribution, fd_b::FragmentDistribution) -> ComparisonResult
+
+Compare the fragment-count distributions of two diagrams.
+"""
+function compare_pair(fd_a::FragmentDistribution, fd_b::FragmentDistribution)
+    sa = fd_a.summary.support
+    sb = fd_b.summary.support
+
+    ComparisonResult(
+        fd_a.diagram.name,
+        fd_b.diagram.name,
+        fd_a.diagram.crossing_number,
+        sa, sb,
+        sa != sb,
+        fd_a.histogram != fd_b.histogram,
+    )
+end
+
+"""
+    compare_by_crossing_number(fds::AbstractVector{FragmentDistribution}) -> Vector{ComparisonResult}
+
+Group diagrams by crossing number, then compare every pair within each group.
+Returns a vector of `ComparisonResult` for all pairs.
+"""
+function compare_by_crossing_number(fds::AbstractVector{FragmentDistribution})
+    # Group by crossing number
+    groups = Dict{Int,Vector{FragmentDistribution}}()
+    for fd in fds
+        cn = fd.diagram.crossing_number
+        push!(get!(Vector{FragmentDistribution}, groups, cn), fd)
+    end
+
+    results = ComparisonResult[]
+    for cn in sort(collect(keys(groups)))
+        group = groups[cn]
+        for i in 1:length(group), j in (i+1):length(group)
+            push!(results, compare_pair(group[i], group[j]))
+        end
+    end
+
+    return results
+end
+
+"""
+    comparison_table(results::Vector{ComparisonResult}) -> DataFrame
+
+Build a comparison DataFrame.
+"""
+function comparison_table(results::Vector{ComparisonResult})
+    rows = map(results) do r
+        (
+            diagram_a = r.diagram_a,
+            diagram_b = r.diagram_b,
+            crossing_number = r.crossing_number,
+            support_a = string(r.support_a),
+            support_b = string(r.support_b),
+            supports_differ = r.supports_differ,
+            histograms_differ = r.histograms_differ,
+        )
+    end
+    return DataFrame(rows)
+end
+
+"""
+    Base.show(io::IO, r::ComparisonResult)
+
+Pretty-print a comparison result.
+"""
+function Base.show(io::IO, r::ComparisonResult)
+    verdict = r.supports_differ ? "DIFFER" : "SAME"
+    print(io, "$(r.diagram_a) vs $(r.diagram_b) [n=$(r.crossing_number)]: ")
+    print(io, "support $(verdict)  ")
+    print(io, "$(r.support_a) vs $(r.support_b)")
+end
diff --git a/KnotDiagrams/src/distribution.jl b/KnotDiagrams/src/distribution.jl
new file mode 100644
index 0000000..a441d2e
--- /dev/null
+++ b/KnotDiagrams/src/distribution.jl
@@ -0,0 +1,87 @@
+# SPDX-License-Identifier: PMPL-1.0-or-later
+#
+# Fragment-count distribution L_D(k) and tropical summary.
+
+"""
+    fragment_distribution(diagram::KnotDiagram) -> FragmentDistribution
+
+Compute the fragment-count distribution for a knot diagram.  This enumerates
+all 2^n resolution states, counts the number of circles (fragments) in each
+state, and builds a histogram mapping fragment count k to the number of states
+that produce exactly k fragments.
+"""
+function fragment_distribution(diagram::KnotDiagram)
+    states = enumerate_states(diagram.gauss_code)
+    histogram = Dict{Int,Int}()
+
+    for (_, circles) in states
+        histogram[circles] = get(histogram, circles, 0) + 1
+    end
+
+    summary = tropical_summary(histogram)
+    total = 1 << diagram.gauss_code.num_crossings
+
+    return FragmentDistribution(diagram, histogram, total, summary)
+end
+
+"""
+    tropical_summary(histogram::Dict{Int,Int}) -> TropicalSummary
+
+Compute the tropical summary of a fragment-count histogram:
+support, min, max, |support|, and classification (spread vs collapse).
+"""
+function tropical_summary(histogram::Dict{Int,Int})
+    support = sort(collect(keys(histogram)))
+    min_k = first(support)
+    max_k = last(support)
+    support_size = length(support)
+    classification = support_size > 1 ? :spread : :collapse
+
+    return TropicalSummary(support, min_k, max_k, support_size, classification)
+end
+
+"""
+    Base.show(io::IO, fd::FragmentDistribution)
+
+Pretty-print a fragment-count distribution.
+"""
+function Base.show(io::IO, fd::FragmentDistribution)
+    println(io, "FragmentDistribution: $(fd.diagram.name)")
+    println(io, "  Crossing number: $(fd.diagram.crossing_number)")
+    println(io, "  Total states: $(fd.total_states)")
+    println(io, "  Histogram (k => count):")
+    for k in sort(collect(keys(fd.histogram)))
+        println(io, "    k=$k : $(fd.histogram[k]) states")
+    end
+    s = fd.summary
+    println(io, "  Tropical summary:")
+    println(io, "    support    = $(s.support)")
+    println(io, "    min(supp)  = $(s.min_k)")
+    println(io, "    max(supp)  = $(s.max_k)")
+    println(io, "    |support|  = $(s.support_size)")
+    print(io,   "    class      = $(s.classification)")
+end
+
+"""
+    distribution_table(fds::AbstractVector{FragmentDistribution}) -> DataFrame
+
+Build a summary DataFrame from a collection of fragment distributions.
+Requires DataFrames.jl.
+"""
+function distribution_table(fds::AbstractVector{FragmentDistribution})
+    rows = map(fds) do fd
+        s = fd.summary
+        (
+            name = fd.diagram.name,
+            crossing_number = fd.diagram.crossing_number,
+            total_states = fd.total_states,
+            histogram = join(["$k=>$(fd.histogram[k])" for k in sort(collect(keys(fd.histogram)))], ", "),
+            support = string(s.support),
+            min_k = s.min_k,
+            max_k = s.max_k,
+            support_size = s.support_size,
+            classification = string(s.classification),
+        )
+    end
+    return DataFrame(rows)
+end
diff --git a/KnotDiagrams/src/knot_table.jl b/KnotDiagrams/src/knot_table.jl
new file mode 100644
index 0000000..c7c4412
--- /dev/null
+++ b/KnotDiagrams/src/knot_table.jl
@@ -0,0 +1,120 @@
+# SPDX-License-Identifier: PMPL-1.0-or-later
+#
+# Hardcoded Gauss codes for prime knots up to 7 crossings.
+#
+# Convention: signed Gauss code where +c = overcrossing at crossing c,
+# -c = undercrossing at crossing c.  Each crossing label appears exactly
+# once positive and once negative.
+#
+# Sources:
+#   - 3_1, 4_1, 5_1: specified in project requirements
+#   - Others: derived from Dowker-Thistlethwaite notation via
+#     standard conversion (see `dt_to_gauss` below).
+#
+# DT notation reference (all prime knots up to 7 crossings):
+#   3_1: [4, 6, 2]           5_1: [6, 8, 10, 2, 4]
+#   4_1: [4, 6, 8, 2]        5_2: [4, 8, 10, 2, 6]
+#   6_1: [4, 8, 12, 2, 10, 6]    6_2: [4, 8, 10, 12, 2, 6]
+#   6_3: [4, 8, 10, 2, 12, 6]
+#   7_1: [8, 10, 12, 14, 2, 4, 6]    7_2: [4, 10, 14, 12, 2, 8, 6]
+#   7_3: [4, 10, 12, 14, 2, 6, 8]    7_4: [6, 10, 12, 14, 2, 4, 8]
+#   7_5: [4, 8, 12, 2, 14, 6, 10]    7_6: [4, 8, 12, 14, 2, 10, 6]
+#   7_7: [4, 8, 10, 14, 2, 12, 6]
+
+"""
+    dt_to_gauss(dt::Vector{Int}) -> Vector{Int}
+
+Convert a Dowker-Thistlethwaite notation to a signed Gauss code.
+
+In DT notation for an n-crossing knot, `dt[i]` is the even position
+paired with odd position `2i - 1`.  If `dt[i] > 0`, the odd-numbered
+encounter is the overcrossing; if `dt[i] < 0`, the even-numbered
+encounter is the overcrossing.
+
+Returns a signed Gauss code (1-indexed crossing labels).
+"""
+function dt_to_gauss(dt::Vector{Int})
+    n = length(dt)
+
+    # Build position → (crossing_label, is_over) map
+    pos_to_crossing = Dict{Int, Tuple{Int, Bool}}()
+    for i in 1:n
+        odd_pos = 2i - 1
+        even_pos = abs(dt[i])
+        is_positive = dt[i] > 0
+
+        # Crossing label = i
+        # If dt[i] > 0: odd position is over, even is under
+        # If dt[i] < 0: even position is over, odd is under
+        pos_to_crossing[odd_pos] = (i, is_positive)       # odd: over if positive
+        pos_to_crossing[even_pos] = (i, !is_positive)     # even: under if positive
+    end
+
+    # Traverse positions 1..2n and build signed Gauss code
+    gauss = Vector{Int}(undef, 2n)
+    for pos in 1:2n
+        label, is_over = pos_to_crossing[pos]
+        gauss[pos] = is_over ? label : -label
+    end
+
+    return gauss
+end
+
+"""
+    prime_knots_up_to_7() -> Vector{KnotDiagram}
+
+Return `KnotDiagram`s for all prime knots up to 7 crossings.
+"""
+function prime_knots_up_to_7()
+    diagrams = KnotDiagram[]
+
+    # --- 3 crossings ---
+    push!(diagrams, KnotDiagram("3_1", 3, GaussCode([1, -2, 3, -1, 2, -3])))
+
+    # --- 4 crossings ---
+    push!(diagrams, KnotDiagram("4_1", 4, GaussCode([1, -2, 3, -4, 2, -1, 4, -3])))
+
+    # --- 5 crossings ---
+    push!(diagrams, KnotDiagram("5_1", 5, GaussCode([1, -2, 3, -4, 5, -1, 2, -3, 4, -5])))
+    push!(diagrams, KnotDiagram("5_2", 5, GaussCode(dt_to_gauss([4, 8, 10, 2, 6]))))
+
+    # --- 6 crossings ---
+    push!(diagrams, KnotDiagram("6_1", 6, GaussCode(dt_to_gauss([4, 8, 12, 2, 10, 6]))))
+    push!(diagrams, KnotDiagram("6_2", 6, GaussCode(dt_to_gauss([4, 8, 10, 12, 2, 6]))))
+    push!(diagrams, KnotDiagram("6_3", 6, GaussCode(dt_to_gauss([4, 8, 10, 2, 12, 6]))))
+
+    # --- 7 crossings ---
+    push!(diagrams, KnotDiagram("7_1", 7, GaussCode(dt_to_gauss([8, 10, 12, 14, 2, 4, 6]))))
+    push!(diagrams, KnotDiagram("7_2", 7, GaussCode(dt_to_gauss([4, 10, 14, 12, 2, 8, 6]))))
+    push!(diagrams, KnotDiagram("7_3", 7, GaussCode(dt_to_gauss([4, 10, 12, 14, 2, 6, 8]))))
+    push!(diagrams, KnotDiagram("7_4", 7, GaussCode(dt_to_gauss([6, 10, 12, 14, 2, 4, 8]))))
+    push!(diagrams, KnotDiagram("7_5", 7, GaussCode(dt_to_gauss([4, 8, 12, 2, 14, 6, 10]))))
+    push!(diagrams, KnotDiagram("7_6", 7, GaussCode(dt_to_gauss([4, 8, 12, 14, 2, 10, 6]))))
+    push!(diagrams, KnotDiagram("7_7", 7, GaussCode(dt_to_gauss([4, 8, 10, 14, 2, 12, 6]))))
+
+    return diagrams
+end
+
+"""
+    torus_knot_gauss(q::Int) -> Vector{Int}
+
+Generate a signed Gauss code for the torus knot T(2, q) where q is odd.
+These follow a regular alternating pattern: the first q entries alternate
+`+1, -2, +3, ...` and the second q entries alternate `-1, +2, -3, ...`.
+"""
+function torus_knot_gauss(q::Int)
+    @assert isodd(q) && q >= 3 "q must be odd and >= 3"
+    n = q  # number of crossings
+
+    code = Vector{Int}(undef, 2n)
+    # First half: 1, -2, 3, -4, ..., ±n
+    for i in 1:n
+        code[i] = iseven(i) ? -i : i
+    end
+    # Second half: -1, 2, -3, 4, ..., ∓n
+    for i in 1:n
+        code[n + i] = iseven(i) ? i : -i
+    end
+
+    return code
+end
diff --git a/KnotDiagrams/src/smoothing.jl b/KnotDiagrams/src/smoothing.jl
new file mode 100644
index 0000000..ba9d65b
--- /dev/null
+++ b/KnotDiagrams/src/smoothing.jl
@@ -0,0 +1,73 @@
+# SPDX-License-Identifier: PMPL-1.0-or-later
+#
+# State enumeration and circle counting via smoothing.
+#
+# For an n-crossing diagram there are 2^n "resolution states", one per
+# binary choice at each crossing.  At crossing c (over-strand at position
+# p, under-strand at position q in the Gauss code):
+#
+#   Smoothing 0 — pair incoming-over with outgoing-under and vice-versa:
+#       union(arc[p-1], arc[q])   and   union(arc[q-1], arc[p])
+#
+#   Smoothing 1 — pair incoming-over with incoming-under and outgoing-over
+#       with outgoing-under:
+#       union(arc[p-1], arc[q-1]) and   union(arc[p], arc[q])
+#
+# Arc indexing: for a Gauss code of length 2n there are 2n arcs.  Arc i
+# (1-indexed) is the strand segment from position i to position i+1
+# (wrapping around: arc 2n goes from position 2n back to position 1).
+
+"""
+    count_circles(gc::GaussCode, state::UInt) -> Int
+
+Given a `GaussCode` and a resolution state (an unsigned integer whose
+i-th bit selects smoothing 0 or 1 at crossing i), return the number
+of disjoint circles in the fully-smoothed diagram.
+"""
+function count_circles(gc::GaussCode, state::UInt)
+    m = gc.num_arcs  # = 2n
+    uf = UnionFind(m)
+
+    for crossing in gc.crossings
+        p = crossing.over_pos
+        q = crossing.under_pos
+
+        # Arc indices (1-indexed, wrapping)
+        arc_before_p = mod1(p - 1, m)  # arc entering position p
+        arc_after_p  = p                # arc leaving position p  (= arc index p)
+        arc_before_q = mod1(q - 1, m)
+        arc_after_q  = q
+
+        bit = (state >> (crossing.label - 1)) & 1
+        if bit == 0
+            # Smoothing 0: cross-connect
+            uf_merge!(uf, arc_before_p, arc_after_q)
+            uf_merge!(uf, arc_before_q, arc_after_p)
+        else
+            # Smoothing 1: parallel-connect
+            uf_merge!(uf, arc_before_p, arc_before_q)
+            uf_merge!(uf, arc_after_p, arc_after_q)
+        end
+    end
+
+    return num_components(uf)
+end
+
+"""
+    enumerate_states(gc::GaussCode) -> Vector{Tuple{UInt, Int}}
+
+Enumerate all 2^n resolution states for the given Gauss code.
+Returns a vector of `(state, num_circles)` pairs.
+"""
+function enumerate_states(gc::GaussCode)
+    n = gc.num_crossings
+    total = 1 << n  # 2^n
+    results = Vector{Tuple{UInt,Int}}(undef, total)
+
+    for s in UInt(0):UInt(total - 1)
+        circles = count_circles(gc, s)
+        results[s + 1] = (s, circles)
+    end
+
+    return results
+end
diff --git a/KnotDiagrams/src/types.jl b/KnotDiagrams/src/types.jl
new file mode 100644
index 0000000..209387d
--- /dev/null
+++ b/KnotDiagrams/src/types.jl
@@ -0,0 +1,109 @@
+# SPDX-License-Identifier: PMPL-1.0-or-later
+#
+# Core types for knot diagram analysis.
+
+"""
+    Crossing
+
+A crossing in a knot diagram, identified by its two positions in the
+Gauss code: `over_pos` (where the strand passes over) and `under_pos`
+(where it passes under). Positions are 1-indexed.
+"""
+struct Crossing
+    label::Int
+    over_pos::Int
+    under_pos::Int
+end
+
+"""
+    GaussCode
+
+A signed Gauss code representation of a knot diagram.
+
+Convention: the code is a vector of signed integers. Positive `+c` means
+the strand passes **over** crossing `c`; negative `-c` means it passes
+**under** crossing `c`. Each crossing label appears exactly twice — once
+positive, once negative.
+
+# Examples
+```julia
+trefoil = GaussCode([1, -2, 3, -1, 2, -3])
+figure_eight = GaussCode([1, -2, 3, -4, 2, -1, 4, -3])
+```
+"""
+struct GaussCode
+    code::Vector{Int}
+    crossings::Vector{Crossing}
+    num_crossings::Int
+    num_arcs::Int
+
+    function GaussCode(code::Vector{Int})
+        n = length(code)
+        @assert iseven(n) "Gauss code must have even length"
+
+        num_crossings = n ÷ 2
+        labels = sort(unique(abs.(code)))
+        @assert length(labels) == num_crossings "Each crossing must appear exactly twice"
+        @assert labels == collect(1:num_crossings) "Crossing labels must be 1:n"
+
+        # Find over/under positions for each crossing
+        over_pos = zeros(Int, num_crossings)
+        under_pos = zeros(Int, num_crossings)
+
+        for (pos, val) in enumerate(code)
+            c = abs(val)
+            if val > 0
+                @assert over_pos[c] == 0 "Crossing $c has two overcrossings"
+                over_pos[c] = pos
+            else
+                @assert under_pos[c] == 0 "Crossing $c has two undercrossings"
+                under_pos[c] = pos
+            end
+        end
+
+        crossings = [Crossing(c, over_pos[c], under_pos[c]) for c in 1:num_crossings]
+
+        for cr in crossings
+            @assert cr.over_pos != 0 && cr.under_pos != 0 "Crossing $(cr.label) missing over or under"
+        end
+
+        new(code, crossings, num_crossings, n)
+    end
+end
+
+"""
+    KnotDiagram
+
+A named knot diagram with its Gauss code and crossing number.
+"""
+struct KnotDiagram
+    name::String
+    crossing_number::Int
+    gauss_code::GaussCode
+end
+
+"""
+    TropicalSummary
+
+Summary statistics for a fragment-count distribution.
+"""
+struct TropicalSummary
+    support::Vector{Int}
+    min_k::Int
+    max_k::Int
+    support_size::Int
+    classification::Symbol  # :spread or :collapse
+end
+
+"""
+    FragmentDistribution
+
+The fragment-count distribution L_D(k) for a knot diagram: a histogram
+mapping fragment count k to the number of states producing k fragments.
+"""
+struct FragmentDistribution
+    diagram::KnotDiagram
+    histogram::Dict{Int,Int}     # k => count of states with k fragments
+    total_states::Int
+    summary::TropicalSummary
+end
diff --git a/KnotDiagrams/src/union_find.jl b/KnotDiagrams/src/union_find.jl
new file mode 100644
index 0000000..bef0565
--- /dev/null
+++ b/KnotDiagrams/src/union_find.jl
@@ -0,0 +1,70 @@
+# SPDX-License-Identifier: PMPL-1.0-or-later
+#
+# Union-Find (Disjoint Set Union) data structure for counting connected
+# components after smoothing a knot diagram. This is the same algorithmic
+# core used by bracket_polynomial in Skein.jl/src/polynomials.jl.
+
+"""
+    UnionFind
+
+Weighted union-find with path compression.
+Tracks the number of distinct components for efficient fragment counting.
+"""
+mutable struct UnionFind
+    parent::Vector{Int}
+    rank::Vector{Int}
+    num_components::Int
+end
+
+"""
+    UnionFind(n::Int) -> UnionFind
+
+Create a union-find structure with `n` elements (1-indexed), each in its own set.
+"""
+function UnionFind(n::Int)
+    UnionFind(collect(1:n), zeros(Int, n), n)
+end
+
+"""
+    find!(uf::UnionFind, x::Int) -> Int
+
+Find the root representative of the set containing `x`, with path compression.
+"""
+function find!(uf::UnionFind, x::Int)
+    while uf.parent[x] != x
+        uf.parent[x] = uf.parent[uf.parent[x]]  # path halving
+        x = uf.parent[x]
+    end
+    return x
+end
+
+"""
+    uf_merge!(uf::UnionFind, x::Int, y::Int) -> Bool
+
+Merge the sets containing `x` and `y`. Returns `true` if they were in
+different sets (i.e., a merge actually occurred).
+
+Named `uf_merge!` rather than `union!` to avoid shadowing `Base.union!`.
+"""
+function uf_merge!(uf::UnionFind, x::Int, y::Int)
+    rx = find!(uf, x)
+    ry = find!(uf, y)
+    rx == ry && return false
+    if uf.rank[rx] < uf.rank[ry]
+        uf.parent[rx] = ry
+    elseif uf.rank[rx] > uf.rank[ry]
+        uf.parent[ry] = rx
+    else
+        uf.parent[ry] = rx
+        uf.rank[rx] += 1
+    end
+    uf.num_components -= 1
+    return true
+end
+
+"""
+    num_components(uf::UnionFind) -> Int
+
+Return the current number of disjoint sets.
+"""
+num_components(uf::UnionFind) = uf.num_components
diff --git a/KnotDiagrams/test/runtests.jl b/KnotDiagrams/test/runtests.jl
new file mode 100644
index 0000000..bdd9cb2
--- /dev/null
+++ b/KnotDiagrams/test/runtests.jl
@@ -0,0 +1,219 @@
+# SPDX-License-Identifier: PMPL-1.0-or-later
+using Test
+using KnotDiagrams
+
+@testset "KnotDiagrams" begin
+
+    @testset "UnionFind" begin
+        uf = UnionFind(6)
+        @test num_components(uf) == 6
+
+        @test uf_merge!(uf, 1, 2)
+        @test num_components(uf) == 5
+        @test find!(uf, 1) == find!(uf, 2)
+
+        @test uf_merge!(uf, 3, 4)
+        @test num_components(uf) == 4
+
+        @test uf_merge!(uf, 1, 3)
+        @test num_components(uf) == 3
+        @test find!(uf, 2) == find!(uf, 4)
+
+        # Merging already-connected elements should not change count
+        @test !uf_merge!(uf, 2, 3)
+        @test num_components(uf) == 3
+    end
+
+    @testset "GaussCode construction" begin
+        # Valid Gauss code (trefoil)
+        gc = GaussCode([1, -2, 3, -1, 2, -3])
+        @test gc.num_crossings == 3
+        @test gc.num_arcs == 6
+
+        # Crossing positions
+        c1 = gc.crossings[1]
+        @test c1.over_pos == 1 && c1.under_pos == 4
+
+        c2 = gc.crossings[2]
+        @test c2.over_pos == 5 && c2.under_pos == 2
+
+        c3 = gc.crossings[3]
+        @test c3.over_pos == 3 && c3.under_pos == 6
+
+        # Figure-eight
+        gc4 = GaussCode([1, -2, 3, -4, 2, -1, 4, -3])
+        @test gc4.num_crossings == 4
+        @test gc4.num_arcs == 8
+
+        # Invalid: odd length
+        @test_throws AssertionError GaussCode([1, -2, 3])
+    end
+
+    @testset "Trefoil circle counts" begin
+        gc = GaussCode([1, -2, 3, -1, 2, -3])
+
+        # Manually verified state (0,0,0) -> 2 circles
+        @test count_circles(gc, UInt(0b000)) == 2
+
+        # Manually verified state (1,1,1) -> 3 circles
+        @test count_circles(gc, UInt(0b111)) == 3
+
+        # Enumerate all 8 states
+        states = enumerate_states(gc)
+        @test length(states) == 8
+
+        # Collect circle counts
+        circle_counts = sort([c for (_, c) in states])
+        # Manually verified: [1, 1, 1, 2, 2, 2, 2, 3]
+        @test circle_counts == [1, 1, 1, 2, 2, 2, 2, 3]
+    end
+
+    @testset "Trefoil fragment distribution" begin
+        diagram = KnotDiagram("3_1", 3, GaussCode([1, -2, 3, -1, 2, -3]))
+        fd = fragment_distribution(diagram)
+
+        @test fd.total_states == 8
+        @test fd.histogram[1] == 3
+        @test fd.histogram[2] == 4
+        @test fd.histogram[3] == 1
+
+        s = fd.summary
+        @test s.support == [1, 2, 3]
+        @test s.min_k == 1
+        @test s.max_k == 3
+        @test s.support_size == 3
+        @test s.classification == :spread
+    end
+
+    @testset "Figure-eight fragment distribution" begin
+        diagram = KnotDiagram("4_1", 4, GaussCode([1, -2, 3, -4, 2, -1, 4, -3]))
+        fd = fragment_distribution(diagram)
+
+        @test fd.total_states == 16
+
+        s = fd.summary
+        # Figure-eight should be spread
+        @test s.classification == :spread
+        @test s.support_size > 1
+        # Should have support different from trefoil (different crossing number,
+        # but verify it's nontrivial)
+        @test s.min_k >= 1
+    end
+
+    @testset "Cinquefoil fragment distribution" begin
+        diagram = KnotDiagram("5_1", 5, GaussCode([1, -2, 3, -4, 5, -1, 2, -3, 4, -5]))
+        fd = fragment_distribution(diagram)
+
+        @test fd.total_states == 32
+
+        s = fd.summary
+        @test s.classification == :spread
+        @test s.support_size > 1
+    end
+
+    @testset "DT to Gauss conversion" begin
+        # Trefoil: DT [4, 6, 2]
+        gauss = dt_to_gauss([4, 6, 2])
+        gc = GaussCode(gauss)
+        @test gc.num_crossings == 3
+
+        # The converted code should produce the same fragment distribution
+        # as the hand-coded trefoil
+        d1 = KnotDiagram("3_1_dt", 3, gc)
+        d2 = KnotDiagram("3_1", 3, GaussCode([1, -2, 3, -1, 2, -3]))
+        fd1 = fragment_distribution(d1)
+        fd2 = fragment_distribution(d2)
+        @test fd1.histogram == fd2.histogram
+
+        # Figure-eight: DT [4, 6, 8, 2]
+        gauss4 = dt_to_gauss([4, 6, 8, 2])
+        gc4 = GaussCode(gauss4)
+        @test gc4.num_crossings == 4
+
+        d3 = KnotDiagram("4_1_dt", 4, gc4)
+        d4 = KnotDiagram("4_1", 4, GaussCode([1, -2, 3, -4, 2, -1, 4, -3]))
+        fd3 = fragment_distribution(d3)
+        fd4 = fragment_distribution(d4)
+        @test fd3.histogram == fd4.histogram
+    end
+
+    @testset "Torus knot Gauss code generator" begin
+        # T(2,3) = trefoil
+        code3 = torus_knot_gauss(3)
+        @test code3 == [1, -2, 3, -1, 2, -3]
+
+        # T(2,5) = cinquefoil
+        code5 = torus_knot_gauss(5)
+        @test code5 == [1, -2, 3, -4, 5, -1, 2, -3, 4, -5]
+
+        # T(2,7)
+        code7 = torus_knot_gauss(7)
+        gc7 = GaussCode(code7)
+        @test gc7.num_crossings == 7
+    end
+
+    @testset "Knot table validity" begin
+        diagrams = prime_knots_up_to_7()
+        @test length(diagrams) == 14  # 1 + 1 + 2 + 3 + 7
+
+        # All diagrams should have valid Gauss codes
+        for d in diagrams
+            @test d.gauss_code.num_crossings == d.crossing_number
+        end
+    end
+
+    @testset "Compare by crossing number" begin
+        diagrams = prime_knots_up_to_7()
+        distributions = [fragment_distribution(d) for d in diagrams]
+        comparisons = compare_by_crossing_number(distributions)
+
+        # Should have comparisons for crossing numbers with > 1 knot
+        @test length(comparisons) > 0
+
+        # 5-crossing: 1 pair (5_1 vs 5_2)
+        cn5_pairs = filter(r -> r.crossing_number == 5, comparisons)
+        @test length(cn5_pairs) == 1
+
+        # 6-crossing: 3 pairs (6_1 vs 6_2, 6_1 vs 6_3, 6_2 vs 6_3)
+        cn6_pairs = filter(r -> r.crossing_number == 6, comparisons)
+        @test length(cn6_pairs) == 3
+
+        # 7-crossing: C(7,2) = 21 pairs
+        cn7_pairs = filter(r -> r.crossing_number == 7, comparisons)
+        @test length(cn7_pairs) == 21
+    end
+
+    @testset "Distribution and comparison tables" begin
+        diagrams = prime_knots_up_to_7()
+        distributions = [fragment_distribution(d) for d in diagrams]
+
+        df = distribution_table(distributions)
+        @test size(df, 1) == length(diagrams)
+        @test "name" in names(df)
+        @test "support" in names(df)
+
+        comparisons = compare_by_crossing_number(distributions)
+        cdf = comparison_table(comparisons)
+        @test size(cdf, 1) == length(comparisons)
+        @test "supports_differ" in names(cdf)
+    end
+
+    @testset "Spread vs collapse — all prime knots are spread" begin
+        diagrams = prime_knots_up_to_7()
+        for d in diagrams
+            fd = fragment_distribution(d)
+            @test fd.summary.classification == :spread "$(d.name) should be spread"
+        end
+    end
+
+    @testset "Different knot types have different distributions" begin
+        # 5_1 vs 5_2 should differ
+        d51 = KnotDiagram("5_1", 5, GaussCode([1, -2, 3, -4, 5, -1, 2, -3, 4, -5]))
+        d52 = KnotDiagram("5_2", 5, GaussCode(dt_to_gauss([4, 8, 10, 2, 6])))
+        fd51 = fragment_distribution(d51)
+        fd52 = fragment_distribution(d52)
+        result = compare_pair(fd51, fd52)
+        @test result.histograms_differ
+    end
+
+end