Nurikabe

Nurikabe is the puzzle of islands in a sea. The grid is a rectangle of cells; some cells carry positive integers. The solver must shade some cells black and leave the rest white, subject to four rules. Every numbered cell is white and serves as the seed of a connected white island of exactly the declared size. Islands of different seeds do not touch each other along an edge. The black cells together form a single connected sea. Nowhere does the sea contain a $2 \times 2$ all-black square.

The puzzle was published by Nikoli in $1991$, with the name Nurikabe drawn from a yokai of Japanese folklore: a white-walled spirit that obstructs travellers’ paths. It joins the family of shading puzzles of which Heyawake, LITS, and Hitori are also members; among them it is the one whose constraints are most balanced between local (no $2 \times 2$, adjacency) and global (per-region connectivity, region size).

The chapter introduces the cleanest CP-SAT pattern for multi-region connectivity: one BFS-distance variable per cell, with the distance equal to zero at the region’s root (a clue cell for islands, a designated cell for the sea), and a reified “has a same-region neighbour with $d - 1$” constraint elsewhere. The same pattern will return in subsequent shading puzzles.

Rules and a small instance

A Nurikabe puzzle is an $R \times C$ grid. A subset of cells carry a positive integer clue indicating the size of the white island whose seed it is. The constraints:

1. Every cell is either black or white.

2. Each numbered cell is white and belongs to a connected white region of exactly the declared size.

3. Two distinct islands share no edge.

4. The set of black cells is connected.

5. No $2 \times 2$ block of cells is entirely black.

Figure 13.1 is puzzle $50$ from the janko.at Nurikabe archive, an $8 \times 8$ instance by Ooya Tate Tadashi. The clues run from $1$ (a singleton island) to $6$.

The programming model

Let $r_{r,c}$ denote the region label of cell $(r, c)$, an integer between $0$ and $N$ where $N$ is the number of clues. The label $0$ marks the sea; labels $1, \ldots, N$ mark the islands in some fixed enumeration of the clue cells. Let $s_{r, c} \in \{0, 1\}$ be the indicator $r_{r, c} = 0$ (“$(r, c)$ is sea”).

The local constraints are direct.

1. Each clue cell $(c_i)$ has $r_{c_i} = i$ and $s_{c_i} = 0$.

2. For each island $i$, $|\{(r, c) : r_{r, c} = i\}| = k_i$ (the clue’s value).

3. For each pair of orthogonally adjacent cells $(p, q)$ with $r_p \ne r_q$, at least one of $s_p, s_q$ is $1$. Equivalently: two whites that are adjacent must share a label.

4. The total number of sea cells is $R \cdot C - \sum_i k_i$.

5. No $2 \times 2$ block is entirely sea: $s_{r, c} + s_{r+1, c} + s_{r, c+1} + s_{r+1, c+1} \le 3$.

Connectivity by BFS-distance

Local constraints alone do not enforce that each region is connected. To do so, we attach a distance variable $d_{r, c} \in \{0, 1, \ldots, R \cdot C - 1\}$ to every cell, intended to hold the BFS distance from the cell to the root of its region.

• For each clue $c_i$, $d_{c_i} = 0$ (the clue is the root of island $i$).

• Exactly one sea cell has $d = 0$ (the sea root); introduce a Boolean $\rho_{r, c}$ for each cell with $\rho_{r, c} = 1 \Rightarrow d_{r, c} = 0 \wedge s_{r, c} = 1$, and $\sum_{r, c} \rho_{r, c} = 1$.

• For every other cell, $d_{r, c} \ge 1$ and there exists an orthogonal neighbour $(r', c')$ with $r_{r', c'} = r_{r, c}$ and $d_{r', c'} + 1 = d_{r, c}$.

The third bullet’s existential is a reified disjunction over the four neighbours; CP-SAT handles it through AddBoolOr guarded by OnlyEnforceIf. The result is that the cells of any region form a tree rooted at the declared root, hence are all reachable from the root, hence connected.

A solver in seventy lines

from ortools.sat.python import cp_model

def solve_nurikabe(R, C, clues):
    """{(r,c): size} -> 2D grid of {0,1}."""
    cs = sorted(clues)
    N = len(cs)
    sea_total = R * C - sum(clues.values())
    m = cp_model.CpModel()
    reg = {(r, c): m.NewIntVar(0, N, "")
           for r in range(R) for c in range(C)}
    sea = {(r, c): m.NewBoolVar("")
           for r in range(R) for c in range(C)}
    for k in reg:
        m.Add(reg[k] == 0).OnlyEnforceIf(sea[k])
        m.Add(reg[k] >= 1).OnlyEnforceIf(sea[k].Not())
    for i, c in enumerate(cs, 1):
        m.Add(reg[c] == i)
    for i, c in enumerate(cs, 1):
        bs = []
        for k in reg:
            b = m.NewBoolVar("")
            m.Add(reg[k] == i).OnlyEnforceIf(b)
            m.Add(reg[k] != i).OnlyEnforceIf(b.Not())
            bs.append(b)
        m.Add(sum(bs) == clues[c])
    m.Add(sum(sea.values()) == sea_total)
    for r in range(R-1):
        for c in range(C-1):
            m.Add(sea[(r,c)] + sea[(r+1,c)]
                  + sea[(r,c+1)] + sea[(r+1,c+1)] <= 3)
    def same_region(p, q):
        b = m.NewBoolVar("")
        m.Add(reg[p] == reg[q]).OnlyEnforceIf(b)
        m.Add(reg[p] != reg[q]).OnlyEnforceIf(b.Not())
        return b
    for r in range(R):
        for c in range(C):
            for dr, dc in [(0,1),(1,0)]:
                nr, nc = r+dr, c+dc
                if not (0 <= nr < R and 0 <= nc < C):
                    continue
                m.AddBoolOr([same_region((r,c), (nr,nc)),
                             sea[(r,c)], sea[(nr,nc)]])
    d = {(r, c): m.NewIntVar(0, R*C-1, "")
         for r in range(R) for c in range(C)}
    for c in cs: m.Add(d[c] == 0)
    rho = {(r, c): m.NewBoolVar("")
           for r in range(R) for c in range(C)}
    for k in rho:
        m.Add(d[k] == 0).OnlyEnforceIf(rho[k])
        m.Add(sea[k] == 1).OnlyEnforceIf(rho[k])
    if sea_total > 0:
        m.Add(sum(rho.values()) == 1)
    def less_dist(p, q):
        b = m.NewBoolVar("")
        m.Add(d[q] + 1 == d[p]).OnlyEnforceIf(b)
        m.Add(d[q] + 1 != d[p]).OnlyEnforceIf(b.Not())
        return b
    cs_set = set(cs)
    for r in range(R):
        for c in range(C):
            if (r, c) in cs_set: continue
            sr = rho[(r, c)]
            m.Add(d[(r,c)] >= 1).OnlyEnforceIf(sr.Not())
            preds = []
            for dr, dc in [(-1,0),(1,0),(0,-1),(0,1)]:
                nr, nc = r+dr, c+dc
                if not (0 <= nr < R and 0 <= nc < C):
                    continue
                same = same_region((r,c), (nr,nc))
                less = less_dist((r,c), (nr,nc))
                pred = m.NewBoolVar("")
                pair = [same, less]
                neg = [same.Not(), less.Not()]
                m.AddBoolAnd(pair).OnlyEnforceIf(pred)
                m.AddBoolOr(neg).OnlyEnforceIf(pred.Not())
                preds.append(pred)
            m.AddBoolOr(preds).OnlyEnforceIf(sr.Not())
    s = cp_model.CpSolver()
    s.Solve(m)
    return [[s.Value(sea[(r, c)])
             for c in range(C)] for r in range(R)]

The toy of Figure 13.1 solves uniquely in about $75$ milliseconds.

A hard instance

Figure 13.3 is puzzle $38$ from the janko.at Nurikabe archive, the hardest $10 \times 10$ instance there, by Koyoppz. With fourteen clues totalling $48$ white cells, the sea has $52$ cells and threads sinuously around several neighbouring islands. CP-SAT returns the unique shading in about $180$ milliseconds.

Sources. Nurikabe was published by Nikoli in $1991$. The name comes from a Japanese folkloric yokai, a white-walled spirit that blocks travellers. The toy is puzzle $50$ in the janko.at Nurikabe archive (Ooya Tate Tadashi); the hard instance is puzzle $38$ (Koyoppz). McPhail established the NP-completeness of Nurikabe in $2003$ via a reduction from $3$-SAT (Carleton University technical report); the BFS-distance connectivity formulation used here is standard in CP-SAT modelling and goes back to the network-design literature.