1"""Exact SAT formulations for zero forcing in P(n,k).23(a) round_sat: the round-based encoding of Proposition P4 with horizon T = 2n - s.4(b) fort_cover: fort-cover formulation (Proposition P2) with lazy fort generation,5 fort minimisation, dihedral images of every fort (sound: automorphisms map6 forts to forts) and first-force normalisation (four directed-edge orbit7 representatives). UNSAT at bound b proves Z > b; the final static CNF is8 re-solved with proof logging and the DRAT proof is checked by drat-trim.9(c) boundary_sat: exists X with |X| = p and |dX| <= b? (Krishnan-Ma-Zhang style10 boundary formula adapted to k = 5), with rotation normalisation x_{u0} v x_{v0}.1112All CNFs are built with python-sat (sequential-counter cardinality encodings).13"""14import json15import os16import subprocess17import time1819import common20from pysat.card import CardEnc, EncType21from pysat.formula import CNF22from pysat.solvers import Solver2324from closure_oracle import closure, closure_fast25from graph_constructor import adjacency, adjacency_masks, dihedral_maps, label, u, v2627DRAT_TRIM = os.path.join(common.ROOT, "code", "tools", "drat-trim", "drat-trim")282930# ----------------------------------------------------------------------------31# helpers32# ----------------------------------------------------------------------------33class VarPool:34 def __init__(self):35 self.top = 03637 def new(self):38 self.top += 139 return self.top404142def atmost(lits, bound, pool):43 enc = CardEnc.atmost(lits=lits, bound=bound, top_id=pool.top, encoding=EncType.seqcounter)44 pool.top = max(pool.top, enc.nv)45 return enc.clauses464748def equals(lits, bound, pool):49 enc = CardEnc.equals(lits=lits, bound=bound, top_id=pool.top, encoding=EncType.seqcounter)50 pool.top = max(pool.top, enc.nv)51 return enc.clauses525354def is_fort(adjmask, F, N):55 """F: bitmask. Fort iff nonempty and every w outside F has 0 or >= 2 neighbours in F."""56 if F == 0:57 return False58 outside = ((1 << N) - 1) & ~F59 while outside:60 low = outside & -outside61 w = low.bit_length() - 162 outside ^= low63 x = adjmask[w] & F64 if x and (x & (x - 1)) == 0:65 return False66 return True676869def minimise_fort(adjmask, F, N):70 """Greedy: remove vertices while the set stays a fort (result is a fort, subset of F)."""71 changed = True72 while changed:73 changed = False74 rest = F75 while rest:76 low = rest & -rest77 rest ^= low78 G = F & ~low79 if G and is_fort(adjmask, G, N):80 F = G81 changed = True82 return F838485def mask_to_list(m):86 out = []87 while m:88 low = m & -m89 out.append(low.bit_length() - 1)90 m ^= low91 return out929394def solve_with_proof(clauses, nv, solver_name, cnf_path, proof_path, time_limit=None):95 """Solves a static CNF with proof logging, writes DIMACS + DRAT, runs drat-trim.96 Returns dict with sat, proof_checked, checker_exit, hashes, times."""97 cnf = CNF(from_clauses=clauses)98 cnf.nv = max(cnf.nv, nv)99 cnf.to_file(cnf_path)100 t0 = time.time()101 with Solver(name=solver_name, bootstrap_with=clauses, with_proof=True) as s:102 sat = s.solve()103 model = s.get_model() if sat else None104 proof = None if sat else s.get_proof()105 t_solve = time.time() - t0106 out = {"sat": sat, "solver": solver_name, "solve_seconds": round(t_solve, 3), "cnf_path": os.path.relpath(cnf_path, common.ROOT),107 "cnf_sha256": common.sha256_file(cnf_path), "num_vars": cnf.nv, "num_clauses": len(clauses)}108 if sat:109 out["model"] = model110 return out111 with open(proof_path, "w") as f:112 f.write("\n".join(proof) + "\n")113 out["proof_path"] = os.path.relpath(proof_path, common.ROOT)114 out["proof_sha256"] = common.sha256_file(proof_path)115 out["proof_lines"] = len(proof)116 t1 = time.time()117 cmd = [DRAT_TRIM, cnf_path, proof_path]118 r = subprocess.run(cmd, capture_output=True, text=True)119 out["checker_cmd"] = " ".join(os.path.relpath(c, common.ROOT) if c.startswith("/") else c for c in cmd)120 out["checker_exit"] = r.returncode121 out["checker_verified"] = ("s VERIFIED" in r.stdout)122 out["checker_seconds"] = round(time.time() - t1, 3)123 out["checker_stdout_tail"] = r.stdout.strip().splitlines()[-3:]124 return out125126127# ----------------------------------------------------------------------------128# (a) round-based encoding (P4)129# ----------------------------------------------------------------------------130def round_sat(n, k, s, solver_name="glucose4", T=None, time_limit=None):131 """Is there a zero-forcing set of size <= s? Horizon T = 2n - s (P4).132 sat is None if the time limit was hit (unknown). Default solver Glucose 4.1133 (independent code base from the CaDiCaL used in the fort-cover / boundary paths)."""134 import threading135 adj = adjacency(n, k)136 N = 2 * n137 if T is None:138 T = 2 * n - s139 pool = VarPool()140 x = [pool.new() for _ in range(N)]141 b = [[pool.new() for _ in range(T + 1)] for _ in range(N)]142 f = {}143 for w in range(N):144 for bb in adj[w]:145 for t in range(1, T + 1):146 f[(bb, w, t)] = pool.new()147 cl = []148 for w in range(N):149 cl.append([-b[w][0], x[w]])150 cl.append([b[w][0], -x[w]])151 for t in range(1, T + 1):152 fs = [f[(bb, w, t)] for bb in adj[w]]153 cl.append([-b[w][t], b[w][t - 1]] + fs)154 cl.append([b[w][t], -b[w][t - 1]])155 for fv in fs:156 cl.append([b[w][t], -fv])157 for bb in adj[w]:158 fv = f[(bb, w, t)]159 cl.append([-fv, b[bb][t - 1]])160 for w2 in adj[bb]:161 if w2 != w:162 cl.append([-fv, b[w2][t - 1]])163 cl.append([b[w][T]])164 cl += atmost(x, s, pool)165 t0 = time.time()166 with Solver(name=solver_name, bootstrap_with=cl) as sv:167 if time_limit:168 timer = threading.Timer(time_limit, lambda: sv.interrupt())169 timer.start()170 sat = sv.solve_limited(expect_interrupt=True)171 timer.cancel()172 else:173 sat = sv.solve()174 model = sv.get_model() if sat else None175 S = None176 if sat:177 S = [w for w in range(N) if model[x[w] - 1] > 0]178 return {"sat": sat, "S": S, "num_vars": pool.top, "num_clauses": len(cl), "T": T,179 "solve_seconds": round(time.time() - t0, 3), "solver": solver_name, "timeout": sat is None}180181182# ----------------------------------------------------------------------------183# (b) fort-cover with lazy forts, symmetry and first-force normalisation184# ----------------------------------------------------------------------------185def first_force_representatives(n, k):186 """Four directed-edge orbit representatives under D_n: (source, target)."""187 return [(u(0, n), v(0, n)), (u(0, n), u(1, n)), (v(0, n), u(0, n)), (v(0, n), v(k, n))]188189190def fort_cover(n, k, b, solver_name="cadical195", time_limit=1200.0, normalise=True, symmetric=True,191 seed_forts=None, verbose=False):192 """Decides whether P(n,k) has a zero-forcing set of size <= b.193 Returns dict: sat (True/False/None on timeout), witness S (if sat), forts (list of vertex lists),194 base clauses, nv, iteration statistics."""195 adj = adjacency(n, k)196 adjm = adjacency_masks(n, k)197 N = 2 * n198 full = (1 << N) - 1199 pool = VarPool()200 svar = [pool.new() for _ in range(N)] # svar[w] = w in S201 base = []202 base += atmost(svar, b, pool)203 if normalise:204 sel = []205 for (src, tgt) in first_force_representatives(n, k):206 c = pool.new()207 sel.append(c)208 base.append([-c, svar[src]])209 base.append([-c, -svar[tgt]])210 for w2 in adj[src]:211 if w2 != tgt:212 base.append([-c, svar[w2]])213 base.append(sel)214 maps = dihedral_maps(n) if symmetric else [list(range(N))]215 fort_clauses = []216 forts = []217 seen = set()218219 def add_fort(F):220 """F bitmask; adds clauses for F and its dihedral images."""221 added = 0222 for p in maps:223 img = 0224 for w in mask_to_list(F):225 img |= 1 << p[w]226 if img in seen:227 continue228 seen.add(img)229 forts.append(img)230 fort_clauses.append([svar[w] for w in mask_to_list(img)])231 added += 1232 return added233234 if seed_forts:235 for F in seed_forts:236 add_fort(F)237 t0 = time.time()238 iters = 0239 result = None240 witness = None241 with Solver(name=solver_name, bootstrap_with=base + fort_clauses) as sv:242 while True:243 if time.time() - t0 > time_limit:244 result = None245 break246 iters += 1247 sat = sv.solve()248 if not sat:249 result = False250 break251 model = sv.get_model()252 S = 0253 for w in range(N):254 if model[svar[w] - 1] > 0:255 S |= 1 << w256 B = closure_fast(adjm, S)257 if B == full:258 result = True259 witness = mask_to_list(S)260 break261 W = full & ~B262 assert is_fort(adjm, W, N)263 F = minimise_fort(adjm, W, N)264 assert is_fort(adjm, F, N) and (F & S) == 0265 before = len(fort_clauses)266 add_fort(F)267 for c in fort_clauses[before:]:268 sv.add_clause(c)269 if verbose and iters % 50 == 0:270 print(" iter %d forts %d elapsed %.1fs" % (iters, len(forts), time.time() - t0), flush=True)271 return {"sat": result, "witness": witness, "forts": [mask_to_list(F) for F in forts], "base_clauses": base,272 "fort_clauses": fort_clauses, "nv": pool.top, "iterations": iters, "wall_seconds": round(time.time() - t0, 3),273 "svar": svar, "solver": solver_name, "normalised": normalise, "symmetric": symmetric}274275276# ----------------------------------------------------------------------------277# (c) boundary formulas278# ----------------------------------------------------------------------------279def boundary_clauses(n, k, p, bmax, normalise=True):280 """CNF: exists X, |X| = p, |dX| <= bmax (z_w >= [w in dX]). Returns (clauses, x, z, nv)."""281 adj = adjacency(n, k)282 N = 2 * n283 pool = VarPool()284 x = [pool.new() for _ in range(N)]285 z = [pool.new() for _ in range(N)]286 cl = []287 for w in range(N):288 for t in adj[w]:289 cl.append([x[w], -x[t], z[w]]) # (not x_w and x_t) -> z_w290 cl.append([-z[w], -x[w]]) # z_w -> not x_w (so z counts only boundary vertices)291 cl += equals(x, p, pool)292 cl += atmost(z, bmax, pool)293 if normalise:294 cl.append([x[u(0, n)], x[v(0, n)]])295 return cl, x, z, pool.top296297298def boundary_of(n, k, X):299 adj = adjacency(n, k)300 Xs = set(X)301 return sorted({w for x in X for w in adj[x] if w not in Xs})302303304def boundary_sat(n, k, p, bmax, solver_name="cadical195"):305 cl, x, z, nv = boundary_clauses(n, k, p, bmax)306 t0 = time.time()307 with Solver(name=solver_name, bootstrap_with=cl) as sv:308 sat = sv.solve()309 model = sv.get_model() if sat else None310 X = None311 if sat:312 X = [w for w in range(2 * n) if model[x[w] - 1] > 0]313 assert len(X) == p and len(boundary_of(n, k, X)) <= bmax314 return {"sat": sat, "X": X, "boundary": (boundary_of(n, k, X) if X else None), "num_vars": nv,315 "num_clauses": len(cl), "solve_seconds": round(time.time() - t0, 3)}316317318def beta(n, k, p, start=11, solver_name="cadical195"):319 """Exact beta_p(n) = min |dX| over |X| = p, if <= start; else returns start+1 meaning '>= start+1'."""320 best = None321 bX = None322 bnd = start323 while True:324 r = boundary_sat(n, k, p, bnd, solver_name)325 if not r["sat"]:326 break327 best = len(r["boundary"])328 bX = r["X"]329 bnd = best - 1330 if bnd < 0:331 break332 return {"n": n, "k": k, "p": p, "beta": best if best is not None else start + 1, "exact": best is not None,333 "X": bX, "boundary": (boundary_of(n, k, bX) if bX else None)}334335336if __name__ == "__main__":337 # smoke tests on known cases338 for (n, k, expect) in [(5, 2, 5), (11, 5, 6)]:339 r1 = round_sat(n, k, expect)340 r0 = round_sat(n, k, expect - 1)341 f1 = fort_cover(n, k, expect)342 f0 = fort_cover(n, k, expect - 1)343 print("P(%d,%d): round_sat(%d)=%s round_sat(%d)=%s fort(%d)=%s fort(%d)=%s iters=%d/%d" % (344 n, k, expect, r1["sat"], expect - 1, r0["sat"], expect, f1["sat"], expect - 1, f0["sat"], f1["iterations"], f0["iterations"]))345 assert r1["sat"] and not r0["sat"] and f1["sat"] and not f0["sat"]346 print("boundary P(19,4) p=16 b=9:", boundary_sat(19, 4, 16, 9)["sat"], " (Krishnan-Ma-Zhang example: expected True)")347 print("boundary P(20,4) p=16 b=9:", boundary_sat(20, 4, 16, 9)["sat"], " (expected False)")348 print("SAT smoke tests passed")349