Checking the four colour theorem with arrays

The four colour theorem says that any map drawn on a plane can be coloured with four colours so that no two countries sharing a border get the same colour. It was the first major theorem proved with a computer, and every proof since has kept a large mechanical core: a list of small graphs, each of which has to be checked, one by one, by a program.

In 2023, for a bachelor project in computer-assisted proofs, four classmates (Abel, Tomas, Lex and Jordy) and I wrote one of those programs in Python: github.com/BonucciAndrea/DurreHeesch, with the mathematics written up in our report (PDF). It runs the Dürre–Heesch reducibility check on Steinberger's 2822 configurations, and it did the job: every one of them checked out as D-reducible, after a twelve-hour run on a high-end PC. One configuration at a time it is not fast, though: about three minutes for a 13-vertex ring on today's machine, and hours once the rings reach 15 or 16. This post rewrites it in Amber. It is not a line-by-line translation. Every step of the algorithm becomes an operation on whole columns: all ring colourings at once, every Kempe interchange at once. On the configurations both versions have finished, the Amber one is over 100 times faster at every ring size. At the bottom of the page is the checker itself, running in your browser.

The proof in one paragraph

Suppose some map needs five colours, and take the smallest one. Two facts about it give a contradiction. First, unavoidability: any such map must contain at least one graph from a fixed list. Steinberger found a list of 2822. Second, reducibility: a smallest counterexample cannot contain any graph from the list, because if it did you could shrink the map, colour the smaller one with four colours, and extend that colouring back. Every configuration on the list must therefore pass a reducibility test. This post is about the test Steinberger's list was built for, called D-reducibility. The other half, unavoidability, is a separate discharging argument, and it is not what this code checks.

A configuration is a small planar graph: a ring of 6 to 16 vertices around some interior vertices. The question is about the ring. Take every way to colour the ring with four colours. Some of those colourings extend straight into the interior. The rest have to be rescued: swap colours along Kempe chains until the ring colouring turns into one that is already known to be good. If every ring colouring is eventually good, the configuration is D-reducible.

What a configuration looks like

U_2822.conf lists each configuration as a header followed by one line per vertex. The first is Birkhoff's diamond, 10 vertices with a 6-ring:

1
10 6 16 0          vertices, ring size, directly extendable colourings
0
 1  3    2  9  6   vertex 1, degree 3, neighbours 2 9 6
 2  4    3  8  9  1
 ...

In Amber a configuration is one list, (V;R;E;adj), with 0-based neighbour lists. The browser version at the bottom of this page carries five of them. The native script reads all 2822.

Step 1: every colouring of the ring, at once

Renaming the colours doesn't change anything, so only essential colourings are counted: colour 0 comes first, and each new colour is the next unused one. For a ring of n there are (3ⁿ⁻¹+5)/8 of them when n is even and (3ⁿ⁻¹−1)/8 when it is odd. That is 1,793,614 for n = 16.

The Python builds them recursively, one list at a time. Amber grows all of them together: each step triples the column (every colouring gets each of the three colours that differ from its last one), then keeps the rows that respect the first-use order. The state is three vectors, a base-4 code, the last colour and the highest colour so far:

step:{c:x 0;l:x 1;m:x 2;nl:(4!l+1),(4!l+2),(4!l+3);m3:m,m,m;k:&~nl>m3+1;(((4*c,c,c)+nl)@k;nl@k;(m3|nl)@k)}
ess:{[n]s:(n-2) step/(,1;,1;,1);c:s[0]@&~0=s 1;{4!_x%y}[c]'|*\1,(n-1)#4}

That builds all 1.79 million ring-16 colourings in 0.12 s. The Python's setup, which builds these colourings and their block tables, takes 53 s.

The other half of the trick is how a colouring is found later. Python stores the good colourings in a tree and walks it, one colour at a time. Amber gives each colouring a number: the steps between neighbouring colours are each 1, 2 or 3, so the colouring reads as a base-3 number. "Is this colouring good?" becomes one array lookup, G@P@rank, which is also one lookup for a million colourings at once.

Step 2: which colourings extend inward

Python's ExtendCheck takes one ring colouring and backtracks through the interior. Amber runs a join instead. It starts with one row per ring colouring. Adding an interior vertex multiplies each row by the colours that vertex can still take: the ones not used by its ring neighbours or by interior neighbours already placed. A vertex whose neighbours are all placed can no longer constrain anything, so its column is dropped and duplicate rows are merged. When the last interior vertex is in, the ring colourings that still have a row are the ones that extend.

Every step is a filter, a gather or a distinct over a column. For a ring-12 configuration the whole step takes 15 ms. On every configuration checked, the number it finds matches the count stored in U_2822.conf exactly, and the test suite checks it against a plain JavaScript backtracking search as well.

Step 3: the Kempe chain game

Pick one of the three ways to split the four colours into two pairs, say {0,1} and {2,3}. Around the ring, consecutive vertices whose colours lie in the same pair form a Kempe sector. Inside the configuration, sectors of the same kind may be joined by Kempe chains; the sectors joined together form a block. Which groupings are possible is a theorem (6.6.4 in Fritsch & Fritsch, The Four-Color Theorem): a grouping of the sectors is a block decomposition exactly when every block holds sectors of one type, no two blocks interleave around the ring, and two blocks that touch at one point touch at exactly one more. Our report (§2.2) turns each decomposition into a sequence of up and down arrows, one per sector, with as many ups as downs and never more downs than ups in any prefix. Those are Dyck words, so a ring with m sectors has Catalan(m/2) decompositions: 1, 2, 5, 14, 42, 132 for 2 to 12 sectors. The code generates them from the arrows, and the test suite checks the result against the theorem directly.

Swapping the two colours of each pair throughout a block gives another proper colouring. A ring colouring becomes good if there is some partition such that, whatever the true block decomposition is, some combination of block swaps turns it into a colouring that is already good. Repeat until nothing changes. The report calls the sets this builds Φ₀ ⊂ Φ₁ ⊂ …; if the last one holds every ring colouring, the configuration is D-reducible. Birkhoff's diamond gets there in five rounds, 16 → 21 → 23 → 27 → 30 → 31 good colourings, which is the report's Φ = Φ₅.

One shortcut is safe and saves half the work. Swapping inside every block of one type just renames two colours everywhere, which gives an equivalent colouring. So one block of each type can always be left alone; the code never swaps the first two blocks.

The swap has a neat form. Colours are 0–3, and the partition that pairs 0 with q+1 swaps exactly the colours that differ by XOR q+1. A Kempe interchange is an XOR on a subset of the ring:

X:0 1 2 3 1 0 3 2 2 3 0 1 3 2 1 0          / XOR table, X@(4*a)+b
c:0 1 2 0 1 3                               / a ring colouring
chain:0 0 1 0 0 1                           / vertices 2 and 5: two {2,3} sectors
c+chain*(X@(4*c)+1)-c                       / swap them under q = 0
0 1 3 0 1 2

That is what makes the game vectorisable. The colourings are grouped by partition and by number of sectors. For each group, every allowed interchange of every block decomposition is laid out as one flat column of (interchange, colouring) cells. Each cell is XORed and relabelled, and its rank is looked up in the good-set. A block decomposition passes if any of its interchanges lands on a good colouring. The colouring passes if all of its decompositions pass. So the whole check is a gather, an any-reduce and an all-reduce.

There is a negative control. A single vertex of degree 6 inside a 6-ring is not D-reducible, and both versions say so: 11 of its 31 ring colourings extend directly, and the Kempe game can't rescue a single one of the other 20.

Making it fast

The first working version was only 10–25 times faster than Python. Four changes took it to 100×:

• A lookup, not a tree. The base-3 rank turns every good-set test into one gather.
• Relabelling as a state machine. After a swap, the colours must be renamed back into first-use order before the lookup. A 2,500-state table does the renaming and the rank in one pass, at two lookups per vertex.
• Cache the interchanges. The interchanges a colouring can reach don't depend on which colourings are good. After the first round their ranks are kept, and later rounds are a single gather against the updated good-set.
• Chunks that fit in cache. The flat columns were first built 4 million cells at a time. Cutting that to 125,000 made the whole program twice as fast. Nothing else changed; the working set stopped spilling out of the CPU cache.

Python vs Amber

Both run on one core of the same machine (WSL, gcc 13, CPython 3.12). The Python rows are the repository's code unchanged; the Amber rows are the native dh.k. Times are summed over the configurations of each ring size among configurations 1–38, the range the Python run finished.

Beyond configuration 38, one configuration per ring size. Python took 17 minutes on the ring-14 one; rings 15 and 16 would take it hours each, so only Amber has been run on those:

Two caveats. First, this compares an array language with plain Python lists, not with NumPy; a careful NumPy port would close part of the gap. Second, the two programs don't do the same amount of work. Python's tree search prunes: it abandons an interchange as soon as the colours read so far leave the tree. Amber builds every allowed interchange in full and throws most of them away. At small rings that costs nothing; at rings 15 and 16 there are many more interchanges, and Amber's time grows about 10× per ring size where Python's grows about 6×. The vectorised version wins by doing more work, much faster.

For scale, the report's own timing table has configuration 27 at 2 min 7 s on a home computer, and checking all 2822 configurations took twelve hours on a high-end PC running several at once. The same configuration takes Amber 0.29 s on one core.

What parallelism adds

The configurations are independent, so the obvious move is to run several at once. With separate Amber processes, configurations 1–200 drop from 128 s to 57 s on four processes. Beyond four it stops helping, for two reasons. The slowest single configuration sets a floor. And the kernel is memory-bound: it streams large integer columns, so extra cores wait on the same memory bus.

Parallelising inside one configuration, with interchange chunks on peach, gave 2.35× on a ring-14 configuration with four threads, and more threads didn't help, for the same memory reason. That result needed one change to the engine. peach hands out work in 1024-item batches, so a job of a few dozen heavy tasks all lands on one thread. The measurement above used a build where the batch size shrinks for small jobs.

The whole program

Everything above is one file of 87 lines, with no dependencies beyond the Amber interpreter. It reads U_2822.conf, builds the colourings and block tables, and defines run1 i, which checks configuration i (0-based) and returns its number, ring size, directly extendable count, the count the file expects, the verdict and the time in ms.

/ ============================================================================
/ dh.k  -  Durre-Heesch D-reducibility check of Steinberger's 2822 configurations
/ Amber port of github.com/BonucciAndrea/DurreHeesch (Dreducing_graphs_final.py)
/ .
/ The Python walks one ring colouring at a time. Here every step works on whole
/ columns: all essential colourings of a ring at once, every Kempe interchange
/ of a batch at once. A set of colourings is a list of n digit vectors (digit j
/ = colour of ring vertex j), one element per colouring.
/ ============================================================================

DIR:"./"                                  / the folder that holds U_2822.conf

/ ---- parse: blocks split on blank lines; line 1 = "V R E 0", then V adjacency rows
L:{x@&~x="\r"}'0:DIR,"U_2822.conf"
cut:&0=#'L
blk:{L@x+!y-x}'[0,1+-1_cut;cut]
blk:blk@&3<#'blk
conf:{h:. x 1;(h 0;h 1;h 2;{-1+2_. x}'x@3+!h 0)}'blk

/ ---- essential colourings of an n-ring (colour 0 then 1, 2 before 3, last <> 0)
step:{c:x 0;l:x 1;m:x 2;nl:(4!l+1),(4!l+2),(4!l+3);m3:m,m,m;k:&~nl>m3+1;(((4*c,c,c)+nl)@k;nl@k;(m3|nl)@k)}
ess:{[n]s:(n-2) step/(,1;,1;,1);c:s[0]@&~0=s 1;{4!_x%y}[c]'|*\1,(n-1)#4}

/ ---- rank of a colouring: its successive colour steps (each 1..3) read in base 3
didx:{[D]n:#D;+/(|*\1,(n-2)#3)*{-1+4!4+x}'(1_D)-(-1_D)}

/ ---- block decompositions: ArrowGen + ArrowToBlock exactly as the Python builds them
bits:{[s]{2!_x%y}[!*/s#2]'|*\1,(s-1)#2}
arrows:{[s]B:bits s;ub:(+\B)-B;db:(!s)-ub;ok:&/(B|~0=db-ub)&(~B)|ub<s%2;+B@\:&ok}
a2b:{{[l1;a]$[a;l1,*&~@[(1+#l1)#0;l1;:;1];l1,$[0=n:l1?*|l1;*|l1;l1 n-1]]}/[,0;-1_x]}
/ every allowed interchange of one decomposition, as a mask over sectors (blocks 0,1 stay)
masks:{[bl]k:1+|/bl;$[k<3;!0;+/{[bl;W;i]$[bl[i]<2;0*W;(*/i#2)*2!_W%*/(bl[i]-2)#2]}[bl;!*/(k-2)#2]'!#bl]}
DEC:{$[x<2;();masks'a2b'arrows x]}'!17
P2:*\1,16#2
/ PAR:1 runs the interchange chunks on peach workers (needs an Amber whose peach grain
/ adapts to small inputs, see README). CELLS = cells per chunk: small enough to stay in cache
PAR:0;CELLS:125000
peach:{`pe(x;y)}                        / as in std.k (a bare script run does not load std.k)
X:0 1 2 3 1 0 3 2 2 3 0 1 3 2 1 0

/ ---- ring-size cache: D = essential colourings, P = rank -> position
CN:0
ring:{[n]$[n=CN;0;[D::ess n;P::@[(*/(n-2)#3)#-1;didx D;:;!#*D];CN::n]];0}

/ ---- which ring colourings extend inside: a join over interior vertices, one
/ column per live interior vertex, rows deduplicated as soon as a vertex is done
stepv:{[v]nb:A v;rn:nb@&nb<R;pa:ACT?nb@&(nb>R-1)&nb<v;
 sel:{[rn;pa;c]m:{[c;m;r]m&~c=(D r)@RID}[c]/[(#RID)#1;rn];&{[c;m;p]m&~c=ACL p}[c]/[m;pa]}[rn;pa]'!4;
 ix:,/sel;RID::RID ix;ACL::(ACL@\:ix),,,/(#'sel)#'!4;ACT::ACT,v;
 keep:&{[v;u]v<|/A u}[v]'ACT;ACT::ACT keep;ACL::ACL keep;k:#ACT;
 $[k;[w:|*\1,(k-1)#4;u:?(RID*4*w 0)++/w*ACL;RID::_u%4*w 0;ACL::{[u;p]4!_u%p}[u]'w];RID::?RID];0}
extend:{RID::!#*D;ACT::!0;ACL::();stepv'R+!V-R;@[(#*D)#0;RID;:;1]}

/ ---- Kempe game. For the rows g (of cand) whose partition-q ring has sn sectors:
/ build every allowed interchange as one flat column, relabel, look it up in G
/ relabel-to-first-occurrence as a state machine: state = which label each colour got
/ (base 5, 0 = none yet); one lookup per vertex gives the label and the next state
Z:!625;ZM:{[z;c](5!_z%*/c#5)-1}[Z]'!4;ZK:+/~0>ZM
NXT:,/+{[d]m:ZM d;Z+(m<0)*(ZK+1)**/d#5}'!4
LAB:,/+{[d]m:ZM d;m+(m<0)*ZK-m}'!4
/ same machine with the previous label folded into the state, emitting the rank digit
Y:!2500;YK:4*_Y%4;YP:4!Y
NX:,/+{[d]k:YK+d;(4*NXT k)+LAB k}'!4
RD:,/+{[d]l:LAB YK+d;(4!4+l-YP)-1}'!4
/ digits of one chunk: vertex j of every (interchange,row) cell, swapped or not, via
/ SWT[cell's mask slot, sector, colour]; the fold relabels and ranks in one pass
kv:{[ds;swt;ft;s;j]k:(4*s 0)+swt@ft+(#ft)#ds j;(NX k;(3*s 1)+RD k)}
/ kcell is pure (reads globals, writes none), so chunks can run on peach workers
kcell:{[q;r]m:#r;T:#TT;ft:64*_(!m*T)%m;ds:{[r;d;i]d[r]+4*i r}[r]'[DC;SID];
 a:!T*64;dd:4!a;b:2!_(TT@_a%64)%P2@_(64!a)%4;swt:dd+b*(X@(4*dd)+q+1)-dd;
 z:0*ft;(CAND r;P@(kv[ds;swt;ft]/[(z+4;z);1+!CN-1])1;MP;m)}
kstore:{[e]$[BUD>#e 1;[BUD::BUD-#e 1;CACHE::CACHE,,e];NOC::1];keval e}
keval:{[e]MPX::e 2;ok:kok[G@e 1;e 3];G::@[G;(e 0)@&ok;:;1];0}
kok:{[gf;m]&/{[gf;m;p]|/{[gf;m;x]gf@(x*m)+!m}[gf;m]'p}[gf;m]'MPX}
kgo:{[q;g;ms]u:?,/ms;TT::u;MP::TT?/:ms;ch:(0N;1|_CELLS%#TT)#g;kstore'$[PAR&1<#ch;peach[kcell[q];ch];kcell[q]'ch];0}
kgroup:{[q;sn;g]ms:DEC sn;$[(0=#ms)|0=&/1,#'ms;0;kgo[q;g;ms]]}
kround:{[q]CAND::&~G;$[#CAND;kr2 q;0]}
kr2:{[q]DC::D@\:CAND;
 out:{[q;x](x=0)|x=q+1}[q]'DC;SID::+\(,(#CAND)#0),~(1_out)=-1_out;
 S:1+*|SID;lo:*|out;SID::{[S;lo;x]x*~lo&x=S-1}[S;lo]'SID;s:S-lo;
 {[q;s;sn]g:&s=sn;$[#g;kgroup[q;sn;g];0]}[q;s]'2*1+!_CN%2;0}

/ ---- one configuration: 1 = D-reducible
run1:{[i]c:conf i;V::c 0;R::c 1;A::c 3;t0:`t[];ring R;G::extend[];e:+/G;
 CACHE::();NOC::0;BUD::200000000;kround'0 1 2;
 {$[NOC;kround'0 1 2;keval'CACHE];+/G}/0;
 ok:(&/G)&e=c 2;(i+1;R;e;c 2;ok;_(`t[]-t0)%1000)}

Put U_2822.conf from the Python repository in the same folder as dh.k, start Amber in that folder (or point DIR on line 11 at it), load the file with \l dh.k, and run, for example, {run1 x}'!30 to check the first thirty configurations. Download the file as dh.k; the version the playground below loads is assets/dh/dh.k, the same code with five configurations written inline instead of read from disk.

Try it

This is the real checker: the Amber interpreter compiled to WebAssembly, running in your browser with the file above loaded. Pick a preset, press Run (or Ctrl+Enter), and edit anything you like. Each line is one statement; lines starting with / are comments.

• Ring colourings builds every essential colouring of an 8-ring (274 of them) and prints the first five, one per row.
• Check a configuration runs the whole algorithm on configuration 27, a 12-ring around 8 interior vertices. It takes a couple of seconds in the browser.
• A graph that fails checks the degree-6 wheel. Its neighbour lists are written out so you can change them and check your own graph.
• Kempe swap = XOR is the swap from the section above.
• Block decompositions counts the block decompositions for 2 to 12 sectors.

check prints one line per field:

A configuration is written (V;R;E;neighbours): V vertices numbered from 0, the first R of them forming the ring in order, E the expected extendable count (or 0N), and then one neighbour list per vertex. The wheel is (7;6;0N;(1 5 6;2 0 6;…;0 1 2 3 4 5)): ring vertex 0 touches its ring neighbours 1 and 5 and the centre 6, and the centre touches every ring vertex. If you edit it, keep the graph a plausible configuration: list every edge from both ends, keep the ring a cycle of at least 6 vertices, and keep every inner face a triangle. The checker will run on anything, but its verdict only means D-reducibility for a real configuration.

Names you can use: c1, c11, c18 and c27 (configurations 1, 11, 18 and 27 of the 2822, rings 6, 10, 11 and 12), wheel, check, and the pieces used above: ess n (all essential colourings of an n-ring, one digit vector per ring vertex), X (the XOR table) and DEC (the interchange masks for each sector count). Try check c11, #*ess 12, or turn the wheel into Birkhoff's diamond by pasting check c1.