stblake · msollami · Jul 4, 2026
diff --git a/Mathilda_spec.md b/Mathilda_spec.md
@@ -51,6 +51,7 @@ Each category lives in [`docs/spec/builtins/`](docs/spec/builtins/):
 | Pattern matching (`MatchQ`, `Cases`, `DeleteCases`, `Position`, `Count`, ...) | [`builtins/pattern-matching.md`](docs/spec/builtins/pattern-matching.md) |
 | File I/O (`Get`, `Put`, `PutAppend`, `>>`, `>>>`) | [`builtins/file-io.md`](docs/spec/builtins/file-io.md) |
 | Graphics (`Graphics`, `Show`, `Plot`, `Point`, `Line`, `Rectangle`, `Circle`, `Disk`, `Polygon`, `Text`, ...) | [`builtins/graphics.md`](docs/spec/builtins/graphics.md) |
+| Graphs (`Graph`, `VertexList`, `EdgeList`, `AdjacencyMatrix`, `FindShortestPath`, `ConnectedComponents`, `PageRank`, `GraphPlot`, ...) | [`builtins/graphs.md`](docs/spec/builtins/graphs.md) |
 | FLINT context (`` FLINT`PolynomialGCD ``, `` FLINT`Factor ``, `` FLINT`Det ``, `` FLINT`Zeta ``, ...) — direct access to the FLINT-backed kernels | [`builtins/flint.md`](docs/spec/builtins/flint.md) |
 
 ## Changelog

diff --git a/SPEC.md b/SPEC.md
@@ -63,6 +63,14 @@ src/
   poly/             Polynomial subsystem (univariate, multivariate, factoring,
                     algebraic-number factoring, polynomial solving)
   linalg/           Dense linear algebra; eigen kernels split by algorithm
+  graph/            Graph subsystem: Graph[] construction/validation, queries
+                    (VertexList/EdgeList/degree/AdjacencyList), matrix views
+                    (Adjacency/Incidence/Kirchhoff/Distance), generators
+                    (Complete/Cycle/Path/Petersen/Random/...), algorithms
+                    (shortest path, components, spanning tree, connectivity,
+                    Euler/Hamilton, cliques, colouring), centrality (PageRank,
+                    Katz, betweenness, clustering) and GraphPlot/Graph3D
+                    (one builtin per file, mirrors linalg/)
   simp/             Simplify, trig simplification, trig rationalisation
   calculus/         D / Dt / Derivative, Series, Limit, Integrate, Risch-Norman
   special_functions/ Higher transcendental & special functions: Gamma, LogGamma,
@@ -265,7 +273,8 @@ return arg0;
 `*_init()` function: `parfrac`, `modular`, `facint`, `comparisons`, `boolean`,
 `list`, `replace`, `patterns`, `cond`, `iter`, `complex`, `trig`, `simp`,
 `hyperbolic`, `logexp`, `piecewise`, `attr`, `purefunc`, `stats`, `poly`,
-`facpoly`, `rat`, `expand`, `info`, `datetime`, `linalg`, `load`, `graphics`.
+`facpoly`, `rat`, `expand`, `info`, `datetime`, `linalg`, `load`, `graphics`,
+`graph`.
 
 After C-side init, `main()` loads `src/internal/init.m`, which `Get[]`s the
 remaining `.m` bootstrap files.

diff --git a/docs/spec/builtins/graphs.md b/docs/spec/builtins/graphs.md
diff --git a/docs/spec/changelog/2026-06-29.md b/docs/spec/changelog/2026-06-29.md
@@ -2,6 +2,69 @@
 
 Feature additions and fixes recorded during this week.
 
+## Graph subsystem: `src/graph/` — construction, queries, matrices, generators, algorithms, drawing (2026-07-04)
+
+New `src/graph/` module (mirrors `src/linalg/`: one builtin per translation
+unit, registered in `graph.c` via `graph_init()`). Graphs are ordinary `Expr`
+trees — `Graph[List[verts], List[edges]]` with each edge a `DirectedEdge[u, v]`
+or `UndirectedEdge[u, v]` — so every generic tool (`Part`, `Map`, `ReplaceAll`,
+pattern matching) works on them unchanged. `Rule`/`->` and `TwoWayRule`/`<->`
+are accepted as parse-time sugar and normalised on construction; `<->` is a new
+operator in `parse.c` at Rule precedence. The constructor validates structure
+(2-arg edges, endpoints present in the vertex list, no self-loops, no parallel
+edges) and returns the canonical form, leaving malformed input unevaluated.
+
+Builtins added, by group:
+
+- **Construction / predicates:** `Graph`, `GraphQ`, `DirectedGraphQ`,
+  `UndirectedGraphQ`, `MixedGraphQ`, `EmptyGraphQ`, `CompleteGraphQ`,
+  `PathGraphQ`, `TreeGraphQ`, `RegularGraphQ`, `AcyclicGraphQ`, `BipartiteQ`,
+  `EulerianGraphQ`, `HamiltonianGraphQ`, `ConnectedGraphQ`,
+  `StronglyConnectedGraphQ`, `IndexGraph`, `Subgraph`, `NeighborhoodGraph`.
+- **Queries:** `VertexList`, `EdgeList`, `VertexCount`, `EdgeCount`,
+  `VertexDegree`, `VertexInDegree`, `VertexOutDegree`, `DegreeSequence`,
+  `AdjacencyList`, `IncidenceList`, `GraphDensity`, `GraphReciprocity`,
+  `GraphAssortativity`.
+- **Matrix views:** `AdjacencyMatrix`, `IncidenceMatrix`, `KirchhoffMatrix`
+  (Laplacian), `DistanceMatrix`.
+- **Structural ops:** `GraphComplement`, `GraphUnion`, `GraphDisjointUnion`,
+  `GraphIntersection`, `GraphDifference`, `GraphJoin`, `GraphProduct`,
+  `GraphPower`, `ReverseGraph`, `LineGraph`, `TransitiveClosureGraph`,
+  `TransitiveReductionGraph`, `VertexAdd`/`VertexDelete`/`VertexContract`,
+  `EdgeAdd`/`EdgeDelete`/`EdgeContract`.
+- **Generators:** `CompleteGraph`, `CompleteKaryTree`, `CycleGraph`,
+  `PathGraph`, `CirculantGraph`, `CocktailPartyGraph`, `TuranGraph`,
+  `KneserGraph`, `GeneralizedPetersenGraph`, `PrismGraph`, `AntiprismGraph`,
+  `LadderGraph`, `GearGraph`, `HelmGraph`, `SunletGraph`, `FriendshipGraph`,
+  `WheelGraph`, dodecahedral/icosahedral, and `RandomGraph`.
+- **Algorithms:** `FindShortestPath`, `GraphDistance`, `ConnectedComponents`,
+  `WeaklyConnectedComponents`, `ConnectedComponentCount`, `FindSpanningTree`,
+  `VertexConnectivity`/`EdgeConnectivity`, `KCoreComponents`, `VertexCoreness`,
+  `FindEulerianCycle`, `FindHamiltonianCycle`/`FindHamiltonianPath`,
+  `FindClique`, `FindIndependentVertexSet`, `FindVertexCover`, `FindEdgeCover`,
+  `FindDominatingSet`, `FindGraphMatching`, `FindVertexColoring`,
+  `ChromaticNumber`, `ChromaticPolynomial`.
+- **Centrality:** `PageRank`, `KatzCentrality`, `BetweennessCentrality`,
+  `EdgeBetweennessCentrality`, `ClosenessCentrality`, `DegreeCentrality`,
+  `LocalClusteringCoefficient`, `MeanClusteringCoefficient`,
+  `GlobalClusteringCoefficient`.
+- **Drawing:** `GraphPlot` (2D node-link diagram → `Graphics[...]`),
+  `HighlightGraph`, `Graph3D`/`GraphPlot3D` (→ `Graphics3D[...]`). A bare valid
+  `Graph`/`Graph3D` result auto-renders as a diagram in the REPL and notebook
+  sidecar, the same way `Graphics` does.
+
+`PageRank` and `KatzCentrality` are computed **exactly** by assembling a
+rational transition/resolvent matrix and calling the existing `LinearSolve`,
+so results are exact rationals rather than floating-point approximations.
+Shared helpers (adjacency build, validity checks, vertex indexing, layout,
+Graphics/Graphics3D emission) live in `graph_util.c`, `layout.c`, `render3d.c`.
+Graph-diagram serialisation was added to `graphics_json.c`
+(`graphics3d_to_plotly_json`, diagram mode for `Disk`/`Point`/`Text`).
+`print.c` renders edges infix (`u -> v` / `u <-> v`) and prints a bare `Graph`
+as a terse `Graph[<n vertices, m edges>]` summary (InputForm/FullForm still
+round-trip through the literal constructor). Every builtin is `Protected`, has
+a docstring, and is covered by `tests/test_graph.c`.
+
 ## RootReduce: WL-faithful algebraic-number canonicalisation (G1–G5) (2026-07-04)
 
 Brought `RootReduce` up to the core Wolfram-Language contract. Previously it only

diff --git a/makefile b/makefile
@@ -15,7 +15,7 @@ else
 endif
 
 CC = gcc
-CFLAGS = -O3 -std=c99 -Wall -Wextra -g -I./src -I./src/list -I./src/linalg -I./src/numbertheory -I./src/poly -I./src/simp -I./src/calculus -I./src/sum -I./src/product -I./src/special_functions -I./src/numerical_calculus -I./src/numerical_roots -I./src/graphics -I/usr/local/include
+CFLAGS = -O3 -std=c99 -Wall -Wextra -g -I./src -I./src/list -I./src/linalg -I./src/numbertheory -I./src/poly -I./src/simp -I./src/calculus -I./src/sum -I./src/product -I./src/special_functions -I./src/numerical_calculus -I./src/numerical_roots -I./src/graphics -I./src/graph -I/usr/local/include
 
 # Readline is available on macOS and Linux but not on Windows (MinGW).
 # Build with USE_READLINE=0 to disable it explicitly (e.g. for cross-builds
@@ -146,7 +146,7 @@ ifeq ($(USE_FLINT), 1)
 endif
 
 SRC_DIR = src
-SRC = $(wildcard $(SRC_DIR)/*.c) $(wildcard $(SRC_DIR)/list/*.c) $(wildcard $(SRC_DIR)/linalg/*.c) $(wildcard $(SRC_DIR)/numbertheory/*.c) $(wildcard $(SRC_DIR)/poly/*.c) $(wildcard $(SRC_DIR)/simp/*.c) $(wildcard $(SRC_DIR)/calculus/*.c) $(wildcard $(SRC_DIR)/sum/*.c) $(wildcard $(SRC_DIR)/product/*.c) $(wildcard $(SRC_DIR)/special_functions/*.c) $(wildcard $(SRC_DIR)/numerical_calculus/*.c) $(wildcard $(SRC_DIR)/numerical_roots/*.c) $(wildcard $(SRC_DIR)/graphics/*.c)
+SRC = $(wildcard $(SRC_DIR)/*.c) $(wildcard $(SRC_DIR)/list/*.c) $(wildcard $(SRC_DIR)/linalg/*.c) $(wildcard $(SRC_DIR)/numbertheory/*.c) $(wildcard $(SRC_DIR)/poly/*.c) $(wildcard $(SRC_DIR)/simp/*.c) $(wildcard $(SRC_DIR)/calculus/*.c) $(wildcard $(SRC_DIR)/sum/*.c) $(wildcard $(SRC_DIR)/product/*.c) $(wildcard $(SRC_DIR)/special_functions/*.c) $(wildcard $(SRC_DIR)/numerical_calculus/*.c) $(wildcard $(SRC_DIR)/numerical_roots/*.c) $(wildcard $(SRC_DIR)/graphics/*.c) $(wildcard $(SRC_DIR)/graph/*.c)
 ifneq ($(USE_GRAPHICS), 1)
 SRC := $(filter-out $(SRC_DIR)/graphics/render.c $(SRC_DIR)/graphics/hershey_font.c, $(SRC))
 endif

diff --git a/src/core.c b/src/core.c
@@ -654,6 +654,8 @@ void core_init(void) {
     zero_test_init();
     void graphics_init(void);
     graphics_init();
+    void graph_init(void);
+    graph_init();
 
     /* Options/SetOptions/OptionValue + the default-options registry. Runs last
      * so every option-name symbol used by the registry is already interned. */

diff --git a/src/graph/acyclic.c b/src/graph/acyclic.c
@@ -0,0 +1,100 @@
+/* acyclic.c - AcyclicGraphQ[g] and TopologicalSort[g].
+ *
+ * Both rest on Kahn's algorithm (topological sort by repeatedly removing
+ * in-degree-0 vertices): a directed graph is acyclic iff every vertex is
+ * eventually removed. O(V+E) over the shared integer adjacency.
+ *
+ *   TopologicalSort[g]  -> a vertex order with every edge pointing forward, or
+ *                          $Failed if g has a directed cycle (or undirected
+ *                          edges, which act as 2-cycles). Directed graphs only,
+ *                          matching the Wolfram Language.
+ *   AcyclicGraphQ[g]    -> True iff g has no cycle. For an all-directed graph
+ *                          that means a DAG (Kahn removes all vertices); for a
+ *                          graph with undirected edges it means the underlying
+ *                          undirected graph is a forest (|E| = |V| - #components).
+ *
+ * Memory (SPEC section 4): returns freshly-allocated results; frees res.
+ */
+
+#include "graph.h"
+#include "expr.h"
+#include "sym_names.h"
+#include <stdlib.h>
+
+/* Kahn's algorithm. Returns a malloc'd array of vertex indices in topological
+ * order (caller frees) and writes the number ordered to *count. When *count
+ * equals a->n the graph is a DAG; otherwise a cycle blocked the remainder. */
+static int* topo_order(const GraphAdj* a, int* count) {
+    int n = a->n;
+    int* indeg = malloc((size_t)(n > 0 ? n : 1) * sizeof(int));
+    int* order = malloc((size_t)(n > 0 ? n : 1) * sizeof(int));
+    int* q     = malloc((size_t)(n > 0 ? n : 1) * sizeof(int));
+    if (!indeg || !order || !q) { free(indeg); free(order); free(q); *count = -1; return NULL; }
+    for (int i = 0; i < n; i++) indeg[i] = a->indeg[i];
+    int head = 0, tail = 0;
+    for (int i = 0; i < n; i++) if (indeg[i] == 0) q[tail++] = i;
+    int k = 0;
+    while (head < tail) {
+        int u = q[head++];
+        order[k++] = u;
+        for (int j = 0; j < a->outdeg[u]; j++) {
+            int w = a->out[u][j];
+            if (--indeg[w] == 0) q[tail++] = w;
+        }
+    }
+    free(indeg); free(q);
+    *count = k;
+    return order;
+}
+
+Expr* builtin_topological_sort(Expr* res) {
+    if (res->data.function.arg_count != 1) return NULL;
+    GraphAdj* a = graph_build_adj(res->data.function.args[0]);
+    if (!a) return NULL;
+    int cnt = 0;
+    int* order = topo_order(a, &cnt);
+    if (!order) { graph_adj_free(a); return NULL; }
+
+    Expr* out;
+    if (cnt == a->n) {
+        Expr** items = (a->n > 0) ? calloc((size_t)a->n, sizeof(Expr*)) : NULL;
+        for (int i = 0; i < a->n; i++)
+            items[i] = expr_copy(a->verts->data.function.args[order[i]]);
+        out = expr_new_function(expr_new_symbol(SYM_List), items, (size_t)a->n);
+        free(items);
+    } else {
+        out = expr_new_symbol(SYM_DollarFailed);   /* cyclic / undirected */
+    }
+    free(order); graph_adj_free(a);
+    return out;
+}
+
+Expr* builtin_acyclic_graph_q(Expr* res) {
+    if (res->data.function.arg_count != 1) return NULL;
+    const Expr* g = res->data.function.args[0];
+    if (!graph_is_valid(g)) return expr_new_symbol(SYM_False);
+
+    const Expr* edges = g->data.function.args[1];
+    size_t ne = edges->data.function.arg_count;
+    int all_directed = 1;
+    for (size_t i = 0; i < ne; i++)
+        if (graph_edge_kind(edges->data.function.args[i]) != SYM_DirectedEdge) { all_directed = 0; break; }
+
+    GraphAdj* a = graph_build_adj(g);
+    if (!a) return NULL;
+
+    int acyclic;
+    if (all_directed) {
+        int cnt = 0;
+        int* order = topo_order(a, &cnt);
+        if (!order) { graph_adj_free(a); return NULL; }
+        acyclic = (cnt == a->n);          /* DAG iff every vertex removed */
+        free(order);
+    } else {
+        /* Underlying undirected graph is a forest iff E = V - components. */
+        int comps = graph_count_components(a, NULL, NULL);
+        acyclic = ((int)ne == a->n - comps);
+    }
+    graph_adj_free(a);
+    return expr_new_symbol(acyclic ? SYM_True : SYM_False);
+}
diff --git a/src/graph/adjgraph.c b/src/graph/adjgraph.c
@@ -0,0 +1,76 @@
+/* adjgraph.c - AdjacencyGraph[m]: build a graph from a 0/1 adjacency matrix.
+ *
+ * The inverse of AdjacencyMatrix. Vertices are the integers 1..n. A symmetric
+ * matrix yields an undirected graph (one UndirectedEdge per i<j with m[i][j]=1);
+ * an asymmetric matrix yields a directed graph (a DirectedEdge for each off-
+ * diagonal m[i][j]=1). Diagonal entries (self-loops) are ignored. The result is
+ * returned as a Graph[...] expression and canonicalized/validated by the
+ * evaluator (builtin_graph).
+ *
+ * Round-trips with AdjacencyMatrix when the source graph's vertices are 1..n.
+ *
+ * Memory (SPEC section 4): returns a freshly-allocated Graph; frees res.
+ */
+
+#include "graph.h"
+#include "expr.h"
+#include "sym_names.h"
+#include <stdlib.h>
+
+/* Integer value of matrix entry, or -1 if not a plain integer. */
+static int entry_int(const Expr* row, size_t j) {
+    const Expr* e = row->data.function.args[j];
+    if (e->type != EXPR_INTEGER) return -1;
+    return (int)e->data.integer;
+}
+
+Expr* builtin_adjacency_graph(Expr* res) {
+    if (res->data.function.arg_count != 1) return NULL;
+    const Expr* m = res->data.function.args[0];
+    if (!graph_is_list(m)) return NULL;
+
+    size_t n = m->data.function.arg_count;
+    /* Validate: n x n, entries 0/1. */
+    for (size_t i = 0; i < n; i++) {
+        const Expr* row = m->data.function.args[i];
+        if (!graph_is_list(row) || row->data.function.arg_count != n) return NULL;
+        for (size_t j = 0; j < n; j++) {
+            int v = entry_int(row, j);
+            if (v != 0 && v != 1) return NULL;
+        }
+    }
+
+    /* Symmetry test. */
+    int symmetric = 1;
+    for (size_t i = 0; i < n && symmetric; i++)
+        for (size_t j = 0; j < n; j++)
+            if (entry_int(m->data.function.args[i], j)
+                != entry_int(m->data.function.args[j], i)) { symmetric = 0; break; }
+
+    /* Vertices 1..n. */
+    Expr** verts = (n > 0) ? calloc(n, sizeof(Expr*)) : NULL;
+    for (size_t i = 0; i < n; i++) verts[i] = expr_new_integer((int64_t)i + 1);
+    Expr* vlist = expr_new_function(expr_new_symbol(SYM_List), verts, n);
+    free(verts);
+
+    /* Edges. Upper-bound count then trim via arg_count. */
+    Expr** edges = (n > 0) ? calloc(n * n, sizeof(Expr*)) : NULL;
+    size_t ne = 0;
+    for (size_t i = 0; i < n; i++) {
+        for (size_t j = 0; j < n; j++) {
+            if (i == j) continue;                    /* ignore self-loops */
+            if (entry_int(m->data.function.args[i], j) != 1) continue;
+            if (symmetric && j < i) continue;        /* one undirected edge per pair */
+            Expr* ea[2] = { expr_new_integer((int64_t)i + 1),
+                            expr_new_integer((int64_t)j + 1) };
+            edges[ne++] = expr_new_function(
+                expr_new_symbol(symmetric ? SYM_UndirectedEdge : SYM_DirectedEdge),
+                ea, 2);
+        }
+    }
+    Expr* elist = expr_new_function(expr_new_symbol(SYM_List), edges, ne);
+    free(edges);
+
+    Expr* gargs[2] = { vlist, elist };
+    return expr_new_function(expr_new_symbol(SYM_Graph), gargs, 2);
+}
diff --git a/src/graph/adjlist.c b/src/graph/adjlist.c
@@ -0,0 +1,71 @@
+/* adjlist.c - AdjacencyList[g] and AdjacencyList[g, v].
+ *
+ * Neighbors of v: successors for directed edges (v -> u yields u), and both
+ * endpoints for undirected edges. This matches the row convention of
+ * AdjacencyMatrix (a 1 in row v, column u means v is adjacent to u). Neighbors
+ * are returned in first-appearance order, de-duplicated.
+ *
+ * AdjacencyList[g]     -> {neighbors(v1), neighbors(v2), ...} in vertex order.
+ * AdjacencyList[g, v]  -> neighbors(v).
+ *
+ * Memory (SPEC section 4): returns freshly-allocated lists; evaluator frees res.
+ */
+
+#include "graph.h"
+#include "expr.h"
+#include "sym_names.h"
+#include <stdlib.h>
+
+/* Build List of neighbors of `v` (deduplicated, first-appearance order). */
+static Expr* neighbors_of(const Expr* g, const Expr* v) {
+    const Expr* edges = g->data.function.args[1];
+    size_t ne = edges->data.function.arg_count;
+    Expr** nbr = (ne > 0) ? calloc(ne * 2, sizeof(Expr*)) : NULL;
+    size_t n = 0;
+
+    for (size_t i = 0; i < ne; i++) {
+        const Expr* e = edges->data.function.args[i];
+        const char* kind = graph_edge_kind(e);
+        const Expr* a = e->data.function.args[0];
+        const Expr* b = e->data.function.args[1];
+        const Expr* add = NULL;
+        if (kind == SYM_UndirectedEdge) {
+            if (expr_eq(a, v)) add = b;
+            else if (expr_eq(b, v)) add = a;
+        } else {
+            if (expr_eq(a, v)) add = b;   /* successor only */
+        }
+        if (add) {
+            int seen = 0;
+            for (size_t j = 0; j < n; j++)
+                if (expr_eq(nbr[j], add)) { seen = 1; break; }
+            if (!seen) nbr[n++] = expr_copy((Expr*)add);
+        }
+    }
+    Expr* list = expr_new_function(expr_new_symbol(SYM_List), nbr, n);
+    free(nbr);
+    return list;
+}
+
+Expr* builtin_adjacency_list(Expr* res) {
+    size_t argc = res->data.function.arg_count;
+    if (argc < 1 || argc > 2) return NULL;
+    const Expr* g = res->data.function.args[0];
+    if (!graph_is_valid(g)) return NULL;
+    const Expr* verts = g->data.function.args[0];
+
+    if (argc == 2) {
+        const Expr* v = res->data.function.args[1];
+        if (graph_vertex_index(verts, v) < 0) return NULL;
+        return neighbors_of(g, v);
+    }
+
+    size_t nv = verts->data.function.arg_count;
+    Expr** rows = (nv > 0) ? calloc(nv, sizeof(Expr*)) : NULL;
+    if (nv > 0 && !rows) return NULL;
+    for (size_t i = 0; i < nv; i++)
+        rows[i] = neighbors_of(g, verts->data.function.args[i]);
+    Expr* list = expr_new_function(expr_new_symbol(SYM_List), rows, nv);
+    free(rows);
+    return list;
+}