---
jupyter:
  description: 'The Steiner tree problem with planex2d: the shortest network joining given points, with Steiner points
    where three edges meet at 120 degrees as constraints of a sketch; the Fermat point of a triangle against
    Torricelli''s construction, the two trees of a rectangle, and the regular pentagon against its closed form.'
---

# The Steiner tree problem

What is the shortest network of roads that joins a set of towns? The roads may meet at junctions of their own, Steiner
points, and in the shortest network every Steiner point joins exactly three roads at angles of 120°. For the four
corners of a square of side 1, the shortest network with no junctions but the corners, the minimum spanning tree, is 3
long, and the Steiner tree with two Steiner points is $1 + \sqrt{3} \approx 2.732$ long. For many points, the problem
is NP-hard, and the number of ways to connect them grows faster than exponentially. A soap film between two glass
plates joined by pins forms a Steiner tree by itself, though not always the shortest one.

The example states the rule of 120° as constraints of a sketch, so that the solver places the Steiner points, and
checks the trees against known lengths, for the Fermat point of a triangle, the two trees of a rectangle and the
tree of a regular pentagon.

![The Steiner tree of a regular pentagon, with three Steiner points](../_static/gallery/steiner-tree.png){.result}

```python
import math

import planex as px
```

```python nbsphinx="hidden"
from docs_figures import Row
```

## A tree as a sketch

`tree` draws a line for every edge, from a rough position of each Steiner point. The towns are fixed. At each Steiner
point, the ends of its three edges coincide, and two angles of 120° between its edges fix the third. The solver then
moves the Steiner points to where the angles hold, and `length` adds up the lines:

```python
def tree(towns, steiner, edges):
    """a sketch of the tree: towns {name: (x, y)} fixed, Steiner points {name: rough (x, y)}, edges [(a, b)]; and
    an end of a line at each Steiner point"""
    at = {**towns, **steiner}
    sk = px.Sketch()
    ends, lines = {s: [] for s in steiner}, {s: [] for s in steiner}
    for a, b in edges:
        if b in towns:
            a, b = b, a
        line = f"{a}_{b}"
        pen = sk.Move(at[a], px.Fix) if a in towns else sk.Move(at[a])
        pen.Line(at[b][0] - at[a][0], at[b][1] - at[a][1], id=line)
        for v, end in ((a, "start"), (b, "end")):
            if v in steiner:
                ends[v].append(f"{line}.{end}")
                lines[v].append(line)
    for s in steiner:
        first, *rest = ends[s]
        for other in rest:
            sk.Constrain(first, other, px.Coincident)
        e1, e2, e3 = lines[s]
        sk.Constrain(e1, e2, px.Angle, 120)
        sk.Constrain(e2, e3, px.Angle, 120)
    return sk, {s: ends[s][0] for s in steiner}


def length(sk):
    """the length of all lines of the sketch"""
    return sum(sk.Measure(line)["length"] for line in sk.elements["Line"])
```

## Three towns

For a triangle whose angles are all below 120°, the one Steiner point is the Fermat point. Torricelli's construction
gives the length of the tree. With the equilateral triangle $BCD$ erected on the side $BC$, outside the triangle, the
tree is as long as $AD$, which the print compares (Figure 1):

```python
A, B, C = (0, 0), (100, 0), (30, 70)
three, junctions = tree({"A": A, "B": B, "C": C}, {"S": (40, 25)}, [("A", "S"), ("B", "S"), ("C", "S")])
turn = math.radians(-60)
D = (B[0] + (C[0] - B[0]) * math.cos(turn) - (C[1] - B[1]) * math.sin(turn),
     B[1] + (C[0] - B[0]) * math.sin(turn) + (C[1] - B[1]) * math.cos(turn))  # C turned about B by -60°
print("status", three.status, "| tree", length(three), "| AD", math.dist(A, D))
```

```python nbsphinx="hide-input"
Row((three, "The Fermat point", {"names": False, "fill": False, "dots": [A, B, C],
                                 "marks": [tuple(three.Location(junctions["S"])[:2])],
                                 "dashed": [px.Line(B, D), px.Line(C, D), px.Line(A, D)],
                                 "labels": [("A", (-6, -6)), ("B", (106, -6)), ("C", (26, 76)), ("D", (131, 100))],
                                 "view": ((-10, 138), (-12, 106))}),
    caption="Three towns, their Steiner point and Torricelli's construction, dashed")
```

## Four towns

The four corners of a rectangle can be joined in two ways with two Steiner points, by a bridge along the long side or
by one across. The example is a rectangle of $a = 150$ by $b = 100$. Like a soap film, the solver finds the shortest tree
of the topology it starts from, and in both topologies every angle of the tree is 120°. Their lengths are $a + \sqrt{3}\,b$ and
$b + \sqrt{3}\,a$, and the bridge along the long side is shorter. Both are compared with these lengths (Figure 2):

```python
rectangle = {"a": (0, 0), "b": (150, 0), "c": (150, 100), "d": (0, 100)}
along, _ = tree(rectangle, {"s": (40, 50), "t": (110, 50)}, [("a", "s"), ("d", "s"), ("s", "t"), ("b", "t"), ("c", "t")])
across, _ = tree(rectangle, {"s": (75, 30), "t": (75, 70)}, [("a", "s"), ("b", "s"), ("s", "t"), ("c", "t"), ("d", "t")])
print("along", length(along), "a + √3 b", 150 + math.sqrt(3) * 100)
print("across", length(across), "b + √3 a", 100 + math.sqrt(3) * 150)
print("minimum spanning tree", 150 + 2 * 100)
```

The tree across the rectangle is even longer than the minimum spanning tree, three of the rectangle's sides. A solver
or a soap film that starts in the wrong topology stays in it.

```python nbsphinx="hide-input"
Row((along, f"Along the long side, {length(along):.1f}", {"names": False, "fill": False, "dots": list(rectangle.values())}),
    (across, f"Across, {length(across):.1f}", {"names": False, "fill": False, "dots": list(rectangle.values())}),
    caption="The two trees of a rectangle")
```

## Five towns

The corners of a regular pentagon need three Steiner points. Its Steiner tree is 3.8911568233 times its side long, a
root of $x^8 - 17x^6 + 29x^4 - 13x^2 + 1$, against 4 sides for the minimum spanning tree. For six towns at the corners
of a regular hexagon, no Steiner point helps, and five of its sides are the shortest network, a result Du, Hwang and
Weng proved for regular polygons in 1987. The pentagon's tree for a side of 100 is compared with the closed form
$\sqrt{17 + 7\sqrt{5} + \sqrt{6\,(65 + 29\sqrt{5})}}\,/\,2$ (Figure 3):

```python
R = 50 / math.sin(math.radians(36))  # the circumradius of a side of 100
corners = {f"t{k}": (R * math.cos(math.radians(90 + 72 * k)), R * math.sin(math.radians(90 + 72 * k))) for k in range(5)}


def toward(p, q, f):
    """the middle of p and q, pulled toward the center by f"""
    return (f * (p[0] + q[0]) / 2, f * (p[1] + q[1]) / 2)


starts = {"s1": toward(corners["t1"], corners["t2"], 0.7), "s2": (0, 0.1 * R), "s3": toward(corners["t3"], corners["t4"], 0.7)}
pentagon, junctions = tree(corners, starts, [("t1", "s1"), ("t2", "s1"), ("s1", "s2"), ("t0", "s2"), ("s2", "s3"),
                                             ("t3", "s3"), ("t4", "s3")])
exact = math.sqrt(17 + 7 * math.sqrt(5) + math.sqrt(6 * (65 + 29 * math.sqrt(5)))) / 2
print("status", pentagon.status, "| tree", length(pentagon), "| closed form", 100 * exact, "| spanning tree", 400)
```

```python nbsphinx="hide-input"
Row((pentagon, "Three Steiner points", {"names": False, "fill": False, "dots": list(corners.values()),
                                         "marks": [tuple(pentagon.Location(e)[:2]) for e in junctions.values()]}),
    caption="The Steiner tree of a regular pentagon")
```

See [Mechanisms: the sketch as a solver](../core/constraints/mechanisms.md) for sketches that compute positions.

```python nbsphinx="hidden"
import sys
import matplotlib.pyplot as plt
sys.path.insert(0, "..")
from docs_figures import INK, PLANEX, RED
from gallery import save_figure
fig, ax = plt.subplots(figsize=(6, 6))
outline = [corners[f"t{k}"] for k in range(5)] + [corners["t0"]]
ax.plot(*zip(*outline), color="#c3cbd3", lw=1.2, ls="--")
for line in pentagon.elements["Line"]:
    ax.plot(*zip(*[pentagon.Location(f"{line}.{end}")[:2] for end in ("start", "end")]), color=PLANEX, lw=3)
ax.scatter(*zip(*corners.values()), s=90, color=INK, zorder=3)
ax.scatter(*zip(*[pentagon.Location(e)[:2] for e in junctions.values()]), s=40, color=RED, zorder=3)
ax.set_aspect("equal")
ax.margins(0.04)
ax.axis("off")
save_figure(fig, "steiner-tree")
plt.close(fig)
```