1// Package sat answers the generator's questions about one Takuzu shape2// with a CDCL SAT solver (gophersat) rather than a hand-rolled backtracker:3// complete a partial shape (randomised), decide whether a partial shape has4// any completion, and enumerate solutions up to a limit. Binary puzzles are5// flat for SAT at every size in play (the literature's 0.3 s at 14x14 is6// the slow end; a shape with fixed edges is milliseconds), where MRV7// backtracking without propagation can spend seconds failing to refute an8// infeasible prefill or to complete a tile with four fixed sides.9//10// Encoding, per shape with rule mask m, over the shape's lines (every11// maximal horizontal or vertical run of its cells):12//13// x(cell) one variable per cell, true = 114// run for every three consecutive cells of a line,15// (a | b | c) and (!a | !b | !c)16// balance every line has exactly len/2 true literals17// (a pair of cardinality constraints)18// unique for every pair of lines of one orientation and equal19// length an auxiliary d(k) <-> x(i,k) xor x(j,k) per20// position, and a clause d(1) | ... | d(len)21// fixed cells unit clauses22//23// The rules package keeps the propagation solver: that one models the24// human technique ladder (pairs, sandwiches, balance, uniqueness,25// bifurcation) and is what forced/guess attribution and difficulty rating26// mean. This package is for questions where only the answer matters.27package sat2829import (30 "math/rand/v2"31 "strings"32 "sync"33 "sync/atomic"34 "time"3536 "mosaic-engine/third_party/gophersat/solver"3738 "mosaic-engine/rules"39)4041// problem is the encoding of one shape under a rule mask, without the42// fixed cells, so it can be built once per (shape, mask) and reused.43// Variables: local index + 1 for cells (absent cells are simply unused),44// auxiliaries above W*H.45type problem struct {46 W, H int47 cells []int // local indexes of the shape's cells48 nbVars int49 constrs []solver.PBConstr50}5152// lineConstraints appends the rules for one set of lines, with x mapping a53// line position to a variable id; next is the first free auxiliary id.54func lineConstraints(lines [][]int, x func(local int) int, mask rules.Mask, next int, cs []solver.PBConstr) (int, []solver.PBConstr) {55 if mask&rules.Run != 0 {56 for _, ln := range lines {57 for k := 0; k+2 < len(ln); k++ {58 cs = append(cs, solver.PropClause(x(ln[k]), x(ln[k+1]), x(ln[k+2])))59 cs = append(cs, solver.PropClause(-x(ln[k]), -x(ln[k+1]), -x(ln[k+2])))60 }61 }62 }63 if mask&rules.Balance != 0 {64 for _, ln := range lines {65 lits := make([]int, len(ln))66 for k, c := range ln {67 lits[k] = x(c)68 }69 cs = append(cs, solver.AtLeast(append([]int(nil), lits...), len(ln)/2), solver.AtMost(lits, len(ln)/2))70 }71 }72 if mask&rules.Unique != 0 {73 // orientation: consecutive local indexes differ by 1 in a row74 row := func(ln []int) bool { return len(ln) == 1 || ln[1] == ln[0]+1 }75 for i := 0; i < len(lines); i++ {76 for j := i + 1; j < len(lines); j++ {77 a, b := lines[i], lines[j]78 if len(a) != len(b) || row(a) != row(b) {79 continue80 }81 ds := make([]int, len(a))82 for k := range a {83 d := next84 next++85 ds[k] = d86 u, v := x(a[k]), x(b[k])87 cs = append(cs,88 solver.PropClause(-d, u, v),89 solver.PropClause(-d, -u, -v),90 solver.PropClause(d, -u, v),91 solver.PropClause(d, u, -v),92 )93 }94 cs = append(cs, solver.PropClause(ds...))95 }96 }97 }98 return next, cs99}100101func encode(sh rules.Shape, mask rules.Mask) *problem {102 p := &problem{W: sh.W, H: sh.H}103 for i, ok := range sh.Mask() {104 if ok {105 p.cells = append(p.cells, i)106 }107 }108 x := func(local int) int { return local + 1 }109 next := sh.W*sh.H + 1110 next, p.constrs = lineConstraints(sh.Lines(), x, mask, next, nil)111 p.nbVars = next - 1112 return p113}114115// solverFor builds a solver for the shape with the given cells fixed and116// extra constraints (blocking clauses), optionally relabelling the shape's117// cell variables (perm maps a local index to the local index whose118// variable it uses) and shuffling the constraint order.119func (p *problem) solverFor(cells []rules.Cell, extra []solver.PBConstr, perm map[int]int, shuffle *rand.Rand) *solver.Solver {120 cs := make([]solver.PBConstr, 0, len(p.constrs)+len(cells)+len(extra))121 // every constraint is copied: gophersat takes ownership of the slices122 // it is given and may reorder them, and the cached encoding is shared123 // by concurrent callers124 relabel := func(c solver.PBConstr) solver.PBConstr {125 lits := make([]int, len(c.Lits))126 for i, l := range c.Lits {127 v := l128 if v < 0 {129 v = -v130 }131 if perm != nil {132 if r, ok := perm[v-1]; ok {133 v = r + 1134 }135 }136 if l < 0 {137 v = -v138 }139 lits[i] = v140 }141 var weights []int142 if c.Weights != nil {143 weights = append([]int(nil), c.Weights...)144 }145 return solver.PBConstr{Lits: lits, Weights: weights, AtLeast: c.AtLeast}146 }147 for _, i := range p.cells {148 v := cells[i]149 if v == rules.Empty {150 continue151 }152 lit := i + 1153 if v == 0 {154 lit = -lit155 }156 cs = append(cs, relabel(solver.PropClause(lit)))157 }158 for _, c := range p.constrs {159 cs = append(cs, relabel(c))160 }161 for _, c := range extra {162 cs = append(cs, relabel(c))163 }164 if shuffle != nil {165 // shuffling the constraint order changes the solver's first166 // decisions too167 shuffle.Shuffle(len(cs), func(a, b int) { cs[a], cs[b] = cs[b], cs[a] })168 }169 return solver.New(solver.ParsePBConstrs(cs))170}171172// model reads the shape's cells back (Empty outside the shape).173func (p *problem) model(s *solver.Solver, perm map[int]int) []rules.Cell {174 m := s.Model()175 out := make([]rules.Cell, p.W*p.H)176 for i := range out {177 out[i] = rules.Empty178 }179 for _, i := range p.cells {180 v := i181 if perm != nil {182 if r, ok := perm[i]; ok {183 v = r184 }185 }186 out[i] = 0187 if v < len(m) && m[v] {188 out[i] = 1189 }190 }191 return out192}193194type key struct {195 W, H int196 mask string197 rule rules.Mask198}199200func keyFor(sh rules.Shape, mask rules.Mask) key {201 var sb strings.Builder202 for _, b := range sh.Mask() {203 if b {204 sb.WriteByte('1')205 } else {206 sb.WriteByte('0')207 }208 }209 return key{sh.W, sh.H, sb.String(), mask}210}211212var (213 cacheMu sync.Mutex214 cache = map[key]*problem{}215)216217// call counters for the generator's timing table218var (219 calls atomic.Int64220 nanos atomic.Int64221)222223// Counters reports SAT solves and the time inside them since the start.224func Counters() (n int64, d time.Duration) {225 return calls.Load(), time.Duration(nanos.Load())226}227228// solve runs one solver. The vendored gophersat keeps its clause-learning229// buffer per solver (upstream: package-level), so solvers run concurrently.230func solve(s *solver.Solver) solver.Status {231 t0 := time.Now()232 defer func() {233 calls.Add(1)234 nanos.Add(int64(time.Since(t0)))235 }()236 return s.Solve()237}238239// problemFor returns the encoding for (shape, mask), built once. Safe for240// concurrent use: the generator builds givens in parallel waves.241func problemFor(sh rules.Shape, mask rules.Mask) *problem {242 k := keyFor(sh, mask)243 cacheMu.Lock()244 defer cacheMu.Unlock()245 if p, ok := cache[k]; ok {246 return p247 }248 p := encode(sh, mask)249 cache[k] = p250 return p251}252253// Feasible reports whether the partial shape (local W*H cells, Empty254// marks blanks) has a completion under the mask. The answer is definitive.255func Feasible(cells []rules.Cell, sh rules.Shape, mask rules.Mask) bool {256 p := problemFor(sh, mask)257 return solve(p.solverFor(cells, nil, nil, nil)) == solver.Sat258}259260// Complete returns one completion of the partial shape under all rules,261// or ok=false when none exists. rng varies the result: the cell variables262// are relabelled by a random permutation of the shape's cells and the263// constraints shuffled, so the solver's first decisions differ per draw;264// on top of that a few blank cells are fixed at random values first and265// dropped again if that makes the shape infeasible.266func Complete(cells []rules.Cell, sh rules.Shape, rng *rand.Rand) (sol []rules.Cell, ok bool) {267 p := problemFor(sh, rules.All)268 order := rng.Perm(len(p.cells))269 perm := make(map[int]int, len(p.cells))270 for i, c := range p.cells {271 perm[c] = p.cells[order[i]]272 }273 var blanks []int274 for _, i := range p.cells {275 if cells[i] == rules.Empty {276 blanks = append(blanks, i)277 }278 }279 rng.Shuffle(len(blanks), func(a, b int) { blanks[a], blanks[b] = blanks[b], blanks[a] })280 seeds := len(blanks) / 8281 for ; seeds >= 0; seeds /= 2 {282 trial := append([]rules.Cell(nil), cells...)283 for _, i := range blanks[:seeds] {284 trial[i] = rules.Cell(rng.IntN(2))285 }286 shuffle := rand.New(rand.NewPCG(rng.Uint64(), rng.Uint64()))287 s := p.solverFor(trial, nil, perm, shuffle)288 if solve(s) == solver.Sat {289 return p.model(s, perm), true290 }291 if seeds == 0 {292 return nil, false293 }294 }295 return nil, false296}297298// Enumerate returns up to limit solutions of the partial shape under the299// mask, by blocking each model found. exhausted is true when fewer than300// limit exist (every solution was returned).301func Enumerate(cells []rules.Cell, sh rules.Shape, mask rules.Mask, limit int) (sols [][]rules.Cell, exhausted bool) {302 p := problemFor(sh, mask)303 var blocks []solver.PBConstr304 for len(sols) < limit {305 s := p.solverFor(cells, blocks, nil, nil)306 if solve(s) != solver.Sat {307 return sols, true308 }309 m := p.model(s, nil)310 sols = append(sols, m)311 lits := make([]int, 0, len(p.cells))312 for _, i := range p.cells {313 if cells[i] != rules.Empty {314 continue315 }316 if m[i] == 1 {317 lits = append(lits, -(i + 1))318 } else {319 lits = append(lits, i+1)320 }321 }322 if len(lits) == 0 {323 return sols, true324 }325 blocks = append(blocks, solver.PropClause(lits...))326 }327 return sols, false328}329330// Count is Enumerate without keeping the models.331func Count(cells []rules.Cell, sh rules.Shape, mask rules.Mask, limit int) (count int, exhausted bool) {332 sols, ex := Enumerate(cells, sh, mask, limit)333 return len(sols), ex334}335336// CompleteRegion completes the target shape jointly with the other shapes337// over one board (row-major, width W): every shape's rules hold at once on338// the shared cells, so the target's completion is one that leaves its339// neighbours a completion too (and its neighbours' completions compatible340// with each other, which per-shape forward checking cannot see). Only the341// target's local W*H cells are returned (Empty outside it); the neighbours342// stay unassigned on the board. rng varies the target's completion as in343// Complete.344func CompleteRegion(board []rules.Cell, W int, others []rules.Shape, target rules.Shape, rng *rand.Rand) (sol []rules.Cell, ok bool) {345 all := append([]rules.Shape{target}, others...)346 var cs []solver.PBConstr347 next := len(board) + 1348 for _, sh := range all {349 x := func(local int) int { return (sh.Y0+local/sh.W)*W + sh.X0 + local%sh.W + 1 }350 next, cs = lineConstraints(sh.Lines(), x, rules.All, next, cs)351 for i, ok := range sh.Mask() {352 if !ok {353 continue354 }355 if v := board[x(i)-1]; v != rules.Empty {356 lit := x(i)357 if v == 0 {358 lit = -lit359 }360 cs = append(cs, solver.PropClause(lit))361 }362 }363 }364 // randomise the target: relabel its cell variables by a permutation and365 // shuffle the constraints, then seed a few of its blanks366 tx := func(local int) int { return (target.Y0+local/target.W)*W + target.X0 + local%target.W + 1 }367 var tcells []int // local indexes of the target's cells368 for i, ok := range target.Mask() {369 if ok {370 tcells = append(tcells, i)371 }372 }373 order := rng.Perm(len(tcells))374 relabel := map[int]int{}375 for i, c := range tcells {376 relabel[tx(c)] = tx(tcells[order[i]])377 }378 for ci := range cs {379 for li, l := range cs[ci].Lits {380 v := l381 if v < 0 {382 v = -v383 }384 if r, ok := relabel[v]; ok {385 if l < 0 {386 r = -r387 }388 cs[ci].Lits[li] = r389 }390 }391 }392 rng.Shuffle(len(cs), func(a, b int) { cs[a], cs[b] = cs[b], cs[a] })393 var blanks []int394 for _, c := range tcells {395 if board[tx(c)-1] == rules.Empty {396 blanks = append(blanks, c)397 }398 }399 rng.Shuffle(len(blanks), func(a, b int) { blanks[a], blanks[b] = blanks[b], blanks[a] })400 extract := func(s *solver.Solver) []rules.Cell {401 m := s.Model()402 out := make([]rules.Cell, target.W*target.H)403 for i := range out {404 out[i] = rules.Empty405 }406 for _, c := range tcells {407 out[c] = 0408 if v := relabel[tx(c)] - 1; v < len(m) && m[v] {409 out[c] = 1410 }411 }412 return out413 }414 seeds := len(blanks) / 8415 for ; seeds >= 0; seeds /= 2 {416 trial := make([]solver.PBConstr, 0, len(cs)+seeds)417 trial = append(trial, cs...)418 for _, c := range blanks[:seeds] {419 lit := relabel[tx(c)]420 if rng.IntN(2) == 0 {421 lit = -lit422 }423 trial = append(trial, solver.PropClause(lit))424 }425 s := solver.New(solver.ParsePBConstrs(copyConstrs(trial)))426 if solve(s) == solver.Sat {427 return extract(s), true428 }429 if seeds == 0 {430 return nil, false431 }432 }433 return nil, false434}435436func copyConstrs(cs []solver.PBConstr) []solver.PBConstr {437 out := make([]solver.PBConstr, len(cs))438 for i, c := range cs {439 out[i] = solver.PBConstr{Lits: append([]int(nil), c.Lits...), AtLeast: c.AtLeast}440 if c.Weights != nil {441 out[i].Weights = append([]int(nil), c.Weights...)442 }443 }444 return out445}446447// CompleteAll completes every shape at once over the board's fixed cells:448// one SAT instance over all of them, randomised like CompleteRegion but449// over every blank cell of every shape. It returns the values of the blank450// cells (by global id) or ok=false when no joint completion exists, which451// is definitive: the layout cannot be filled against its context. For a452// chunk of a few dozen tiles this replaces tile-by-tile filling with453// backjumping, which can thrash when a late tile's only redrawable454// neighbours cannot help.455func CompleteAll(board []rules.Cell, W int, shapes []rules.Shape, rng *rand.Rand) (assign map[int]rules.Cell, ok bool) {456 var cs []solver.PBConstr457 next := len(board) + 1458 var blanks []int // global ids of blank shape cells, deduplicated459 seen := map[int]bool{}460 for _, sh := range shapes {461 x := func(local int) int { return (sh.Y0+local/sh.W)*W + sh.X0 + local%sh.W + 1 }462 next, cs = lineConstraints(sh.Lines(), x, rules.All, next, cs)463 for i, ok := range sh.Mask() {464 if !ok {465 continue466 }467 id := x(i) - 1468 if v := board[id]; v != rules.Empty {469 lit := x(i)470 if v == 0 {471 lit = -lit472 }473 cs = append(cs, solver.PropClause(lit))474 } else if !seen[id] {475 seen[id] = true476 blanks = append(blanks, id)477 }478 }479 }480 if len(blanks) == 0 {481 return map[int]rules.Cell{}, true482 }483 // relabel the blank cells' variables by a permutation, shuffle the484 // constraints, seed a few blanks485 order := rng.Perm(len(blanks))486 relabel := map[int]int{}487 for i, id := range blanks {488 relabel[id+1] = blanks[order[i]] + 1489 }490 for ci := range cs {491 for li, l := range cs[ci].Lits {492 v := l493 if v < 0 {494 v = -v495 }496 if r, ok := relabel[v]; ok {497 if l < 0 {498 r = -r499 }500 cs[ci].Lits[li] = r501 }502 }503 }504 rng.Shuffle(len(cs), func(a, b int) { cs[a], cs[b] = cs[b], cs[a] })505 shuffled := append([]int(nil), blanks...)506 rng.Shuffle(len(shuffled), func(a, b int) { shuffled[a], shuffled[b] = shuffled[b], shuffled[a] })507 seeds := len(blanks) / 8508 for ; seeds >= 0; seeds /= 2 {509 trial := make([]solver.PBConstr, 0, len(cs)+seeds)510 trial = append(trial, cs...)511 for _, id := range shuffled[:seeds] {512 lit := relabel[id+1]513 if rng.IntN(2) == 0 {514 lit = -lit515 }516 trial = append(trial, solver.PropClause(lit))517 }518 s := solver.New(solver.ParsePBConstrs(copyConstrs(trial)))519 if solve(s) == solver.Sat {520 m := s.Model()521 out := make(map[int]rules.Cell, len(blanks))522 for _, id := range blanks {523 out[id] = 0524 if v := relabel[id+1] - 1; v < len(m) && m[v] {525 out[id] = 1526 }527 }528 return out, true529 }530 if seeds == 0 {531 return nil, false532 }533 }534 return nil, false535}