diff --git a/papers/face_monochromatic_pairs/experiments/check_30_residual.py b/papers/face_monochromatic_pairs/experiments/check_30_residual.py
new file mode 100644
index 0000000..40c962d
--- /dev/null
+++ b/papers/face_monochromatic_pairs/experiments/check_30_residual.py
@@ -0,0 +1,256 @@
+"""Verify the structural characterization of the 30 residual
+chord-apex+Kempe colourings on which the G'-pentagon fallback
+empirically fails (= |S| = 8 AND all G'-pentagons are hit by S).
+
+Hypothesis: these are exactly the colourings on triangulations where
+v's cyclic neighbour-degree sequence has the form (5, 7, 7, 5, 5)
+(or cyclic shifts), and the reduction index i is chosen so that the
+two "7"s land on the flank positions (n_i, n_{i+1}) = (7, 7).
+
+For each of the 30 colourings, report:
+  - parent triangulation order n_G,
+  - cyclic degree sequence of v's 5 neighbours,
+  - reduction index i,
+  - p_G and # G'-pentagons hit by S,
+  - the actual deciding face's length and type.
+
+If the hypothesis holds, all 30 are on (5, 7, 7, 5, 5)-type
+configurations.
+
+Run with: sage experiments/check_30_residual.py
+"""
+import os
+import sys
+import time
+
+from sage.all import Graph
+from sage.graphs.graph_generators import graphs
+
+HERE = os.path.dirname(os.path.abspath(__file__))
+sys.path.insert(0, HERE)
+
+from check_conj_3_8_scaled import (
+    apply_reduction,
+    proper_3_edge_colorings,
+    matches_chord_apex_kempe,
+    kempe_cycle_set,
+    edge_idx,
+)
+from check_heawood_on_kempe import dual_of, vertices_of_kempe
+
+
+def is_g_prime_pentagon(f, named):
+    if len(f) != 5: return False
+    fset = {frozenset(e) for e in f}
+    return not (named['side_0'] in fset or named['side_1'] in fset
+                or named['spike'] in fset or named['merged'] in fset)
+
+
+def classify_face(f, named):
+    fset = {frozenset(e) for e in f}
+    has_s0 = named['side_0'] in fset
+    has_s1 = named['side_1'] in fset
+    has_sp = named['spike'] in fset
+    has_m = named['merged'] in fset
+    if has_s0 and has_sp and not has_s1 and not has_m: return 'flank-lower'
+    if has_sp and has_s1 and not has_s0 and not has_m: return 'flank-upper'
+    if has_s0 and has_s1 and has_m and not has_sp: return 'outer'
+    if has_m and not (has_s0 or has_s1 or has_sp): return 'merged'
+    if not (has_s0 or has_s1 or has_sp or has_m): return 'G-prime'
+    return 'mixed'
+
+
+def test_one(D, G_parent):
+    D.is_planar(set_embedding=True)
+    G_parent.is_planar(set_embedding=True)
+    emb_G = G_parent.get_embedding()
+    examples = []   # 30 residual colourings
+    # We need to map dual face → parent vertex (= the face's "dual vertex")
+    # For each dual face we want to identify the parent vertex v.
+    # The dual face's vertices correspond to triangular faces of G_parent
+    # incident to v.
+    # Build a face→parent-vertex map.
+    parent_faces = G_parent.faces()
+    face_to_dual_vertex_idx = {}
+    # The dual vertices are indexed by face number.
+    # Mapping: parent_face[fi] is a face, and dual vertex fi corresponds.
+    # We need the reverse: given a degree-5 vertex v in G_parent, its
+    # corresponding dual face = the face whose vertices are the 5
+    # parent_face-indices that contain v.
+    v_to_dual_face = {}
+    for v in G_parent.vertex_iterator():
+        if G_parent.degree(v) != 5: continue
+        # Find the 5 parent_face-indices that contain v
+        contains = []
+        for fi, pf in enumerate(parent_faces):
+            for e in pf:
+                if v in e:
+                    contains.append(fi)
+                    break
+        if len(contains) != 5: continue
+        v_to_dual_face[v] = contains   # 5 dual vertex indices forming F_v boundary
+
+    for face in D.faces():
+        if len(face) != 5: continue
+        face_verts = [e[0] for e in face]
+        # Identify which parent v this face corresponds to
+        v_parent = None
+        for v, dual_verts in v_to_dual_face.items():
+            if set(face_verts) == set(dual_verts):
+                v_parent = v
+                break
+        if v_parent is None: continue
+        # Get the cyclic degree sequence of v_parent's neighbours
+        nbrs_cyclic = emb_G[v_parent]
+        cyc_degs = [G_parent.degree(u) for u in nbrs_cyclic]
+
+        for i_red in range(5):
+            res = apply_reduction(D, face, i_red, 9999)
+            if res is None: continue
+            H = res['H']; named = res['named']
+            H.is_planar(set_embedding=True)
+            edges, colorings = proper_3_edge_colorings(H)
+            cand = [c for c in colorings
+                    if matches_chord_apex_kempe(edges, c, named)]
+            v_n = 9999
+            for col in cand:
+                # Identify bad sub-case (ii.B)
+                target = {named['side_0'], named['spike']}
+                lower_flank = None
+                for f in H.faces():
+                    if target.issubset({frozenset(e) for e in f}):
+                        lower_flank = f; break
+                if lower_flank is None or len(lower_flank) != 5: continue
+                arc_verts = [e[0] for e in lower_flank]
+                if v_n not in arc_verts: continue
+                k = arc_verts.index(v_n)
+                cyc = arc_verts[k:] + arc_verts[:k]
+                A_i = next(iter(named['side_0'] - {v_n}))
+                A_ip1 = next(iter(named['spike'] - {v_n}))
+                if cyc[1] == A_i and cyc[4] == A_ip1:
+                    P_1, P_2 = cyc[2], cyc[3]
+                elif cyc[1] == A_ip1 and cyc[4] == A_i:
+                    P_2, P_1 = cyc[2], cyc[3]
+                else: continue
+                merged_idx = edge_idx(edges, named['merged'])
+                c_col = col[merged_idx]
+                c_0_col = col[edge_idx(edges, named['side_0'])]
+                c_1_col = col[edge_idx(edges, named['side_1'])]
+                e_AiP1 = edge_idx(edges, frozenset((A_i, P_1)))
+                e_P1P2 = edge_idx(edges, frozenset((P_1, P_2)))
+                if e_AiP1 is None or e_P1P2 is None: continue
+                if col[e_AiP1] != c_1_col or col[e_P1P2] != c_0_col:
+                    continue
+                a = c_col
+                other = [x for x in range(3) if x != a]
+                kc_b = kempe_cycle_set(edges, col, merged_idx, (a, other[0]))
+                kc_c = kempe_cycle_set(edges, col, merged_idx, (a, other[1]))
+                V_b = vertices_of_kempe(edges, kc_b)
+                V_c = vertices_of_kempe(edges, kc_c)
+                V_union = V_b | V_c
+                S = set(H.vertices()) - V_union
+                if P_1 in V_union: continue
+                if len(S) != 8: continue
+
+                # Check if hit = p_G = 8 (= all G'-pentagons hit)
+                p_total = 0
+                p_hit = 0
+                for f in H.faces():
+                    if not is_g_prime_pentagon(f, named): continue
+                    p_total += 1
+                    verts = {u for (u, v) in f} | {v for (u, v) in f}
+                    if verts & S: p_hit += 1
+                if p_total != 8 or p_hit != 8: continue
+
+                # This is a residual case. Find the actual deciding face.
+                deciding_faces = []
+                for f in H.faces():
+                    verts = {u for (u, v) in f} | {v for (u, v) in f}
+                    if not verts.issubset(V_union): continue
+                    L = len(f)
+                    if L % 3 == 0: continue
+                    deciding_faces.append((classify_face(f, named), L))
+
+                # Compute deg sequence around v_parent (cyclic)
+                examples.append({
+                    'cyc_degs': cyc_degs,
+                    'i_red': i_red,
+                    'deciding_faces': deciding_faces,
+                    'p_total': p_total,
+                    'p_hit': p_hit,
+                    'S_size': len(S),
+                })
+    return examples
+
+
+def main(max_n=20, time_budget_per_n=1800):
+    print("Verifying the 30 residual chord-apex+Kempe colourings.\n")
+    grand_examples = []
+    for n in range(12, max_n + 1):
+        start = time.time()
+        try:
+            triangulations = list(graphs.triangulations(n, minimum_degree=5))
+        except Exception as ex:
+            print(f"n={n}: cannot enumerate ({ex})")
+            continue
+        n_count = 0
+        for tri_idx, G in enumerate(triangulations):
+            if time.time() - start > time_budget_per_n:
+                print(f"  n={n}: timeout at tri {tri_idx}")
+                break
+            G.is_planar(set_embedding=True)
+            D = dual_of(G)
+            exs = test_one(D, G)
+            for ex in exs:
+                ex['n_G'] = n
+                ex['tri_idx'] = tri_idx
+            n_count += len(exs)
+            grand_examples.extend(exs)
+        elapsed = time.time() - start
+        print(f"n={n}: {n_count} residual colourings [{elapsed:.0f}s]")
+        sys.stdout.flush()
+    print()
+    print("=" * 70)
+    print(f"Total residual colourings: {len(grand_examples)}")
+    # Count by cyclic deg sequence (canonical = min rotation)
+    def canon(t):
+        t = tuple(t)
+        best = t
+        n = len(t)
+        for r in range(1, n):
+            rot = t[r:] + t[:r]
+            if rot < best: best = rot
+        # Also reflections
+        rev = t[::-1]
+        for r in range(n):
+            rot = rev[r:] + rev[:r]
+            if rot < best: best = rot
+        return best
+    canon_dist = {}
+    for ex in grand_examples:
+        key = canon(ex['cyc_degs'])
+        canon_dist[key] = canon_dist.get(key, 0) + 1
+    print("\nCyclic deg sequence distribution (canonical = lex-min "
+          "over rotations + reflections):")
+    for k in sorted(canon_dist, key=lambda x: -canon_dist[x]):
+        c = canon_dist[k]
+        print(f"  {k}: {c}")
+    # Decision face distribution
+    df_dist = {}
+    for ex in grand_examples:
+        for (t, L) in ex['deciding_faces']:
+            df_dist[(t, L)] = df_dist.get((t, L), 0) + 1
+    print("\nDeciding face (type, length) distribution:")
+    for k in sorted(df_dist, key=lambda x: -df_dist[x]):
+        t, L = k
+        print(f"  {t}, |f|={L}: {df_dist[k]}")
+    # Print first 5 examples
+    print("\nFirst 5 examples:")
+    for ex in grand_examples[:5]:
+        print(f"  n={ex['n_G']}, tri#{ex['tri_idx']}, i={ex['i_red']}, "
+              f"cyc_degs={ex['cyc_degs']}, "
+              f"deciding={ex['deciding_faces']}")
+
+
+if __name__ == '__main__':
+    main()
diff --git a/papers/face_monochromatic_pairs/experiments/check_30_residual_v2.py b/papers/face_monochromatic_pairs/experiments/check_30_residual_v2.py
new file mode 100644
index 0000000..4086549
--- /dev/null
+++ b/papers/face_monochromatic_pairs/experiments/check_30_residual_v2.py
@@ -0,0 +1,214 @@
+"""Simpler version: enumerate the residual (|S|=8, hit=p_total=8)
+colourings without requiring v_parent identification. For each,
+determine:
+  - the n_k sequence of the reduction (from the reduced dual's
+    F_k structure),
+  - whether any G'-face (length ≢ 0 mod 3) is uncovered.
+
+The earlier check_S_face_structure.py showed at most 30 |S|=8 cases
+with hit = 8, but didn't constrain p_total. Of those, ≤30 have
+p_total = 8 (= the residual).
+
+Run with: sage experiments/check_30_residual_v2.py
+"""
+import os
+import sys
+import time
+
+from sage.all import Graph
+from sage.graphs.graph_generators import graphs
+
+HERE = os.path.dirname(os.path.abspath(__file__))
+sys.path.insert(0, HERE)
+
+from check_conj_3_8_scaled import (
+    apply_reduction,
+    proper_3_edge_colorings,
+    matches_chord_apex_kempe,
+    kempe_cycle_set,
+    edge_idx,
+)
+from check_heawood_on_kempe import dual_of, vertices_of_kempe
+
+
+def is_g_prime_pentagon(f, named):
+    if len(f) != 5: return False
+    fset = {frozenset(e) for e in f}
+    return not (named['side_0'] in fset or named['side_1'] in fset
+                or named['spike'] in fset or named['merged'] in fset)
+
+
+def is_g_prime_face(f, named):
+    fset = {frozenset(e) for e in f}
+    return not (named['side_0'] in fset or named['side_1'] in fset
+                or named['spike'] in fset or named['merged'] in fset)
+
+
+def test_one(D):
+    D.is_planar(set_embedding=True)
+    residual_colourings = []
+    other_bad = []
+    for face in D.faces():
+        if len(face) != 5: continue
+        for i_red in range(5):
+            res = apply_reduction(D, face, i_red, 9999)
+            if res is None: continue
+            H = res['H']; named = res['named']
+            H.is_planar(set_embedding=True)
+            edges, colorings = proper_3_edge_colorings(H)
+            cand = [c for c in colorings
+                    if matches_chord_apex_kempe(edges, c, named)]
+            v_n = 9999
+            # Compute the n_k sequence for this reduction
+            # n_k = degree of the face F_k of original G' = length of
+            # the corresponding face in the reduced dual after subtracting
+            # the reduction effect.
+            # Simpler: the flank face F^♭_{i, i+1} has length n_i - 1.
+            # So n_i = length of F^♭_{i,i+1} + 1.
+            # Get all faces:
+            named_face_lengths = {}
+            for f in H.faces():
+                fset = {frozenset(e) for e in f}
+                if named['side_0'] in fset and named['spike'] in fset:
+                    named_face_lengths['flank_lower'] = len(f)
+                if named['spike'] in fset and named['side_1'] in fset:
+                    named_face_lengths['flank_upper'] = len(f)
+                if (named['side_0'] in fset and named['side_1'] in fset
+                    and named['merged'] in fset):
+                    named_face_lengths['outer'] = len(f)
+                if (named['merged'] in fset and named['side_0'] not in fset
+                    and named['side_1'] not in fset
+                    and named['spike'] not in fset):
+                    named_face_lengths['merged'] = len(f)
+            n_i = named_face_lengths.get('flank_lower', 0) + 1   # = n_i, where the flank covers
+            n_ip1 = named_face_lengths.get('flank_upper', 0) + 1
+            # n_{i+3} = F_merged length + 2
+            n_ip3 = named_face_lengths.get('merged', 0) + 2
+            # n_{i+2} + n_{i+4} from outer
+            outer_len = named_face_lengths.get('outer', 0)
+            # outer_len = n_{i+2} + n_{i+4} - 3
+
+            for col in cand:
+                target = {named['side_0'], named['spike']}
+                lower_flank = None
+                for f in H.faces():
+                    if target.issubset({frozenset(e) for e in f}):
+                        lower_flank = f; break
+                if lower_flank is None or len(lower_flank) != 5: continue
+                arc_verts = [e[0] for e in lower_flank]
+                if v_n not in arc_verts: continue
+                k = arc_verts.index(v_n)
+                cyc = arc_verts[k:] + arc_verts[:k]
+                A_i = next(iter(named['side_0'] - {v_n}))
+                A_ip1 = next(iter(named['spike'] - {v_n}))
+                if cyc[1] == A_i and cyc[4] == A_ip1:
+                    P_1, P_2 = cyc[2], cyc[3]
+                elif cyc[1] == A_ip1 and cyc[4] == A_i:
+                    P_2, P_1 = cyc[2], cyc[3]
+                else: continue
+                merged_idx = edge_idx(edges, named['merged'])
+                c_col = col[merged_idx]
+                c_0_col = col[edge_idx(edges, named['side_0'])]
+                c_1_col = col[edge_idx(edges, named['side_1'])]
+                e_AiP1 = edge_idx(edges, frozenset((A_i, P_1)))
+                e_P1P2 = edge_idx(edges, frozenset((P_1, P_2)))
+                if e_AiP1 is None or e_P1P2 is None: continue
+                if col[e_AiP1] != c_1_col or col[e_P1P2] != c_0_col:
+                    continue
+                a = c_col
+                other = [x for x in range(3) if x != a]
+                kc_b = kempe_cycle_set(edges, col, merged_idx, (a, other[0]))
+                kc_c = kempe_cycle_set(edges, col, merged_idx, (a, other[1]))
+                V_b = vertices_of_kempe(edges, kc_b)
+                V_c = vertices_of_kempe(edges, kc_c)
+                V_union = V_b | V_c
+                S = set(H.vertices()) - V_union
+                if P_1 in V_union: continue
+                # bad colouring
+                # Count G'-pentagons total and hit
+                p_total = 0
+                p_hit = 0
+                non_pent_uncovered = []
+                for f in H.faces():
+                    if not is_g_prime_face(f, named): continue
+                    L = len(f)
+                    verts = {u for (u, v) in f} | {v for (u, v) in f}
+                    if L == 5:
+                        p_total += 1
+                        if verts & S: p_hit += 1
+                    else:
+                        if L % 3 != 0 and verts.issubset(V_union):
+                            non_pent_uncovered.append(L)
+                # Is this a "residual" case?
+                S_size = len(S)
+                if S_size == 8 and p_total == p_hit:
+                    residual_colourings.append({
+                        'n_i': n_i, 'n_ip1': n_ip1, 'n_ip3': n_ip3,
+                        'outer_len': outer_len,
+                        'p_total': p_total, 'p_hit': p_hit,
+                        'S_size': S_size,
+                        'non_pent_uncovered': non_pent_uncovered,
+                    })
+                elif S_size == 8 and p_hit == p_total - 1 and p_total >= 7:
+                    # Border case: only 1 pentagon uncovered
+                    pass
+
+    return residual_colourings
+
+
+def main(max_n=20, time_budget_per_n=1800):
+    print("Detailed analysis of |S|=8, p_hit = p_total residual "
+          "chord-apex+Kempe colourings.\n")
+    grand_residuals = []
+    for n in range(12, max_n + 1):
+        start = time.time()
+        try:
+            triangulations = list(graphs.triangulations(n, minimum_degree=5))
+        except Exception as ex:
+            print(f"n={n}: cannot enumerate ({ex})")
+            continue
+        n_count = 0
+        for tri_idx, G in enumerate(triangulations):
+            if time.time() - start > time_budget_per_n:
+                print(f"  n={n}: timeout at tri {tri_idx}")
+                break
+            G.is_planar(set_embedding=True)
+            D = dual_of(G)
+            resids = test_one(D)
+            for r in resids:
+                r['n_G'] = n
+                r['tri_idx'] = tri_idx
+            n_count += len(resids)
+            grand_residuals.extend(resids)
+        elapsed = time.time() - start
+        print(f"n={n}: {n_count} residual colourings [{elapsed:.0f}s]")
+        sys.stdout.flush()
+    print()
+    print("=" * 70)
+    print(f"Total residual colourings: {len(grand_residuals)}")
+    if grand_residuals:
+        # n_k sequence distribution
+        seq_dist = {}
+        for r in grand_residuals:
+            key = (r['n_i'], r['n_ip1'], r['n_ip3'], r['outer_len'])
+            seq_dist[key] = seq_dist.get(key, 0) + 1
+        print("\n(n_i, n_{i+1}, n_{i+3}, F_outer length) distribution:")
+        for k, c in sorted(seq_dist.items(), key=lambda x: -x[1]):
+            print(f"  {k}: {c}")
+        # For each, does a non-pentagon G'-face provide a deciding face?
+        has_non_pent = sum(1 for r in grand_residuals if r['non_pent_uncovered'])
+        print(f"\nResidual colourings with at least one length-≢0-mod-3 "
+              f"G'-face uncovered: {has_non_pent} / {len(grand_residuals)} "
+              f"({100*has_non_pent/len(grand_residuals):.2f}%)")
+        # Lengths of those non-pent G'-faces
+        len_dist = {}
+        for r in grand_residuals:
+            for L in r['non_pent_uncovered']:
+                len_dist[L] = len_dist.get(L, 0) + 1
+        print("Lengths of uncovered non-pentagon G'-faces:")
+        for L, c in sorted(len_dist.items()):
+            print(f"  |f| = {L}: {c}")
+
+
+if __name__ == '__main__':
+    main()
diff --git a/papers/face_monochromatic_pairs/experiments/check_gprime_pentagon_always_works.py b/papers/face_monochromatic_pairs/experiments/check_gprime_pentagon_always_works.py
new file mode 100644
index 0000000..c8d2273
--- /dev/null
+++ b/papers/face_monochromatic_pairs/experiments/check_gprime_pentagon_always_works.py
@@ -0,0 +1,150 @@
+"""Verify: across all 1,314 bad chord-apex+Kempe colourings, does
+at least one G'-pentagon (length 5) have its boundary entirely in
+V(K_b) ∪ V(K_c)?
+
+(Earlier characterize_S_vertices.py and check_30_residual_v2.py
+suggested yes, but let's verify directly.)
+
+This is the G'-pentagon fallback conjecture for the empirically
+identified bad cases.
+
+Run with: sage experiments/check_gprime_pentagon_always_works.py
+"""
+import os
+import sys
+import time
+
+from sage.all import Graph
+from sage.graphs.graph_generators import graphs
+
+HERE = os.path.dirname(os.path.abspath(__file__))
+sys.path.insert(0, HERE)
+
+from check_conj_3_8_scaled import (
+    apply_reduction,
+    proper_3_edge_colorings,
+    matches_chord_apex_kempe,
+    kempe_cycle_set,
+    edge_idx,
+)
+from check_heawood_on_kempe import dual_of, vertices_of_kempe
+
+
+def test_one(D):
+    D.is_planar(set_embedding=True)
+    bad_count = 0
+    bad_with_uncov_pent = 0
+    examples_no_uncov = []
+    for face in D.faces():
+        if len(face) != 5: continue
+        for i_red in range(5):
+            res = apply_reduction(D, face, i_red, 9999)
+            if res is None: continue
+            H = res['H']; named = res['named']
+            H.is_planar(set_embedding=True)
+            edges, colorings = proper_3_edge_colorings(H)
+            cand = [c for c in colorings
+                    if matches_chord_apex_kempe(edges, c, named)]
+            v_n = 9999
+            for col in cand:
+                target = {named['side_0'], named['spike']}
+                lower_flank = None
+                for f in H.faces():
+                    if target.issubset({frozenset(e) for e in f}):
+                        lower_flank = f; break
+                if lower_flank is None or len(lower_flank) != 5: continue
+                arc_verts = [e[0] for e in lower_flank]
+                if v_n not in arc_verts: continue
+                k = arc_verts.index(v_n)
+                cyc = arc_verts[k:] + arc_verts[:k]
+                A_i = next(iter(named['side_0'] - {v_n}))
+                A_ip1 = next(iter(named['spike'] - {v_n}))
+                if cyc[1] == A_i and cyc[4] == A_ip1:
+                    P_1, P_2 = cyc[2], cyc[3]
+                elif cyc[1] == A_ip1 and cyc[4] == A_i:
+                    P_2, P_1 = cyc[2], cyc[3]
+                else: continue
+                merged_idx = edge_idx(edges, named['merged'])
+                c_col = col[merged_idx]
+                c_0_col = col[edge_idx(edges, named['side_0'])]
+                c_1_col = col[edge_idx(edges, named['side_1'])]
+                e_AiP1 = edge_idx(edges, frozenset((A_i, P_1)))
+                e_P1P2 = edge_idx(edges, frozenset((P_1, P_2)))
+                if e_AiP1 is None or e_P1P2 is None: continue
+                if col[e_AiP1] != c_1_col or col[e_P1P2] != c_0_col:
+                    continue
+                a = c_col
+                other = [x for x in range(3) if x != a]
+                kc_b = kempe_cycle_set(edges, col, merged_idx, (a, other[0]))
+                kc_c = kempe_cycle_set(edges, col, merged_idx, (a, other[1]))
+                V_b = vertices_of_kempe(edges, kc_b)
+                V_c = vertices_of_kempe(edges, kc_c)
+                V_union = V_b | V_c
+                S = set(H.vertices()) - V_union
+                if P_1 in V_union: continue
+                bad_count += 1
+                # Search for an UNCOVERED G'-pentagon
+                found_uncov_pent = False
+                for f in H.faces():
+                    if len(f) != 5: continue
+                    fset = {frozenset(e) for e in f}
+                    # Must be G'-pentagon (not adjacent to F_v reduction)
+                    if (named['side_0'] in fset or named['side_1'] in fset
+                        or named['spike'] in fset or named['merged'] in fset):
+                        continue
+                    verts = {u for (u, v) in f} | {v for (u, v) in f}
+                    if verts.issubset(V_union):
+                        found_uncov_pent = True
+                        break
+                if found_uncov_pent:
+                    bad_with_uncov_pent += 1
+                else:
+                    if len(examples_no_uncov) < 3:
+                        examples_no_uncov.append({
+                            'S_size': len(S),
+                        })
+    return bad_count, bad_with_uncov_pent, examples_no_uncov
+
+
+def main(max_n=20, time_budget_per_n=1800):
+    print("Direct check: G'-pentagon fallback on bad colourings.\n")
+    grand_bad = 0
+    grand_with = 0
+    grand_ex = []
+    for n in range(12, max_n + 1):
+        start = time.time()
+        try:
+            triangulations = list(graphs.triangulations(n, minimum_degree=5))
+        except Exception as ex:
+            print(f"n={n}: cannot enumerate ({ex})")
+            continue
+        n_bad = 0; n_with = 0
+        for tri_idx, G in enumerate(triangulations):
+            if time.time() - start > time_budget_per_n:
+                print(f"  n={n}: timeout")
+                break
+            G.is_planar(set_embedding=True)
+            D = dual_of(G)
+            nb, nw, exs = test_one(D)
+            n_bad += nb; n_with += nw
+            grand_ex.extend(exs)
+        elapsed = time.time() - start
+        print(f"n={n}: {n_bad} bad, {n_with} with uncovered G'-pentagon "
+              f"[{elapsed:.0f}s]")
+        sys.stdout.flush()
+        grand_bad += n_bad
+        grand_with += n_with
+    print()
+    print("=" * 70)
+    print(f"Total bad colourings: {grand_bad}")
+    print(f"With at least one uncovered G'-pentagon: {grand_with} "
+          f"({100*grand_with/max(grand_bad, 1):.2f}%)")
+    if grand_bad - grand_with > 0:
+        print(f"WITHOUT uncovered G'-pentagon: {grand_bad - grand_with}")
+        print("Examples:")
+        for ex in grand_ex:
+            print(f"  |S| = {ex['S_size']}")
+
+
+if __name__ == '__main__':
+    main()