Calcudoku: the obvious way to read the Jacobson-Matthews chain is biased, and a census of 16.9 million Latin squares shows where

I built Calcudoku (the KenKen-style arithmetic puzzle) in the browser with a six-rung solver. The puzzle part went fine. The interesting part turned out to be the step before the puzzle: drawing the Latin square the board is cut from.

The standard uniform sampler is the Jacobson–Matthews Markov chain, and the standard way to read a square off it — run a while, and if the state is improper, keep stepping until it is proper — is biased. On 7×7 it gives a mean of 10.110 intercalates (2×2 Latin subsquares) against an exact 10.529, and running ten times longer gives 10.107. It is not a mixing problem. It is a when do you look problem, and it has a one-line explanation. Reading the chain at a fixed time instead lands at 10.570.

"Exact" means exact: I enumerated all 16,942,080 reduced Latin squares of order 7 and counted intercalates in every one. And intercalates turn out to matter for the puzzle as well: an intercalate that no cage notices is a second answer for free.

Demo: https://sen.ltd/portfolio/calcudoku/ · Source: https://github.com/sen-ltd/calcudoku

The rules

• Fill an n×n grid with 1 to n so that every row and every column holds each number exactly once — a Latin square.
• Heavy outlines cut the grid into cages. Each cage shows a target and an operation; its numbers must reach the target with that operation.
• Minus and divide only label two-cell cages, and the pair can be in either order.
• A number may repeat inside a cage, as long as the repeats share no row or column.
• A one-cell cage gives its number away.

I checked the two rules that change the code (minus/divide only on pairs; repeats allowed across lines) against a primary source rather than memory. KenKen is a registered trademark of KenKen Puzzle LLC, so the project uses the generic name.

Drawing a Latin square uniformly

A board is a Latin square, cut into cages, labelled. If the square is not uniform, the bank quietly leans toward some structures.

The generator people write first fills cells in reading order, tries shuffled values, and backtracks. It reaches every square, but not with equal probability. The standard fix is the Jacobson–Matthews chain (1996): view the square as an n×n×n 0/1 cube with exactly one 1 on every axis-parallel line, and move by adding and subtracting 1 on the eight corners of a sub-box:

M[at(x, y, z)]++;   M[at(x, y1, z1)]++;
M[at(x1, y, z1)]++; M[at(x1, y1, z)]++;
M[at(x, y1, z)]--;  M[at(x, y, z1)]--;
M[at(x1, y, z)]--;
const k = at(x1, y1, z1);
M[k]--;
this.bad = M[k] < 0 ? k : -1;   // one -1 left behind: improper

Every line still sums to one, but a move can leave a single -1 behind: an improper cube, which is not a square. The next move then starts from that cell. At stationarity the chain is uniform over the proper cubes.

"Step until proper" is biased

The obvious reader:

runUntilProper(steps: number, rnd: () => number): void {
  for (let i = 0; i < steps; i++) this.step(rnd);
  while (this.bad >= 0) this.step(rnd);
}

My generator used it at first. A dry run of the stats tool (500 draws) put the 6×6 intercalate mean 1.1 below the exact value. At 500 draws that could have been waved off as noise, so I reproduced it on the order with the biggest gap: 5,000 draws after 216, 1,000 and 5,000 moves gave 7.491, 7.464 and 7.448; one long thinned chain gave 7.498. More moves did not move it. Under-mixing does not behave like that.

The real run, 20,000 draws per row:

grid	sampler	mean intercalates	exact	total variation	noise floor (mean / 95th pct)	chains per square
6×6	naive backtracking	7.348	8.265	0.0942	0.0065 / 0.0095	—
6×6	J–M, 216 moves, then wait for proper	7.513	8.265	0.0769	0.0065 / 0.0095	—
6×6	J–M, 2,160 moves, then wait for proper	7.443	8.265	0.0827	0.0065 / 0.0095	—
6×6	J–M, exactly 216 moves, keep only if proper	8.263	8.265	0.0086	0.0065 / 0.0095	5.66
7×7	naive backtracking	9.909	10.529	0.0659	0.0101 / 0.0145	—
7×7	J–M, 343 moves, then wait for proper	10.110	10.529	0.0473	0.0101 / 0.0145	—
7×7	J–M, 3,430 moves, then wait for proper	10.107	10.529	0.0472	0.0101 / 0.0145	—
7×7	J–M, exactly 343 moves, keep only if proper	10.570	10.529	0.0105	0.0101 / 0.0145	6.83

The noise floor is the total variation distance you get from 20,000 perfect draws from the exact distribution, simulated 40 times. At orders 5, 6 and 7, the fixed-time reader sits inside the 95th percentile of that floor and the wait-for-proper reader sits outside it.

Why: it is about when you look

Waiting for the chain to become proper samples it at the moment it re-enters the proper squares. At stationarity, the flow into a square from improper cubes equals the flow out of it into improper cubes, and that is proportional to the share of its moves that are not proper-to-proper.

What is a proper-to-proper move? Read it off the update. From a proper cube, pick a 0-cell (x,y,z); the result is proper only if the eighth corner (x1,y1,z1) was already a 1 — and then the four cells (x,y), (x,y1), (x1,y), (x1,y1) read a b / b a. A proper-to-proper move is exactly an intercalate flip. A square with k intercalates has 4k of them out of n³ − n² moves. Squares with few intercalates spend more of their traffic on the improper route, get re-entered more often, and are over-drawn — the same direction as the table.

The fix is to condition at a fixed time: run exactly T moves, keep the square if the cube is proper at that moment, and otherwise throw the chain away and start another.

export function sampleSquare(n: number, rnd: () => number, steps = jmSteps(n)): Square {
  for (let k = 0; ; k++) {
    // e.g. a single move from the cyclic square of odd order is never proper
    if (k === 100000) throw new Error(`no proper square after ${steps} moves in ${k} chains`);
    const jm = new JacobsonMatthews(n);
    jm.run(steps, rnd);
    if (jm.proper) return jm.square();
  }
}

That costs 6.83 chains per square on 7×7; the stationary chain is proper only about 15% of the time.

The guard in that comment has a story. I added a T = 1 row to the table, and the stats tool never finished. The cyclic square (r + c) mod n of odd order has no intercalates, so its first move is always improper, and "throw it away and try again" loops forever. I kept the too-short rows in the table (exactly 7 moves on 7×7: mean 7.685, total variation 0.4626), because they show the contrast: the under-mixing bias goes away with more moves, and the reader's bias does not.

The census it is held to

"Exact" is not an estimate. The ledger enumerates every reduced Latin square (first row and first column in order) of order 1 to 7 and counts intercalates in each. Row and column permutations never change the intercalate count, and every square is a row-and-column permutation of exactly one reduced square in exactly n!(n−1)! ways, so the reduced histogram is the distribution over all squares.

order	reduced squares	all squares	mean intercalates	share with none	most
5	56	161,280	3.571	10.71%	4
6	9,408	812,851,200	8.265	0.43%	27
7	16,942,080	61,479,419,904,000	10.529	0.10%	42

Order 7 takes 41 s. The counts match OEIS A000315 and A002860, and every stats run recomputes orders 1–6 from scratch and stops if the stored histogram disagrees. The heavy order runs as its own script that prints each order as it lands; tables are copied from that output and never from memory.

Counting intercalates has a cheap trick. Fix two rows and send each column to the column where the second row holds the same symbol; an intercalate is exactly a 2-cycle of that map. That is n steps per row pair instead of a quadruple loop:

for (let r1 = 0; r1 < n; r1++)
  for (let r2 = r1 + 1; r2 < n; r2++)
    for (let c = 0; c < n; c++) {
      const c2 = pos[r2 * n + L[r1 * n + c]];
      if (c2 > c && L[r1 * n + c2] === L[r2 * n + c]) k++;
    }

An intercalate is a second answer waiting for a cage to miss it

Flip an intercalate and the square is still Latin. So if every cage touching those four cells still reaches its target, the board has a second answer, full stop. A cage misses the flip when it holds one cell from each diagonal of the 2×2 (its multiset of numbers does not change), or when it is a two-cell minus or divide cage whose other number sits at the same distance, or the same ratio, from both symbols.

That is a check that costs nothing next to a search. Of 16,000 boards cut at random (orders 4–7 × largest cage 2–5 × 1,000 each), 8,704 had more than one answer, and on 6,658 of them (76.5%) a flipped intercalate was already one of the extra answers. The converse is a theorem, so it is a test, not a sample: 0 boards with a surviving flip ever came out unique.

It also predicts which squares make good boards. On all 16 rows, squares whose raw cut came out unique had fewer intercalates than squares whose cut did not:

grid	largest cage	unique	non-unique explained by a flip	intercalates (unique)	(not unique)	exact mean
6×6	2	580 / 1,000	335 / 420 (79.8%)	7.31 ± 0.16	9.51 ± 0.25	8.27
6×6	4	397 / 1,000	465 / 603 (77.1%)	6.93 ± 0.18	9.36 ± 0.20	8.27
7×7	3	433 / 1,000	444 / 567 (78.3%)	9.63 ± 0.16	11.28 ± 0.16	10.53
7×7	5	317 / 1,000	521 / 683 (76.3%)	9.92 ± 0.19	11.08 ± 0.14	10.53

(All 16 rows are on the demo page and in the README.) So a generator that draws a uniform square, cuts cages, and keeps only the cuts that come out unique ships a bank that leans away from intercalates — biased, even with a perfect sampler underneath. That the sampler's bias and the filter's bias point the same way is probably not a coincidence: both mechanisms react to the same movable 2×2.

Repair, don't reject

So this generator never throws a square away. It cuts random cages and labels them, then repairs: while a second answer exists, split a cage holding a cell where the two answers differ. Once unique, it tightens: try merging neighbouring cages and keep a merge only if the answer stays unique, until a full pass is refused (locally maximal). Every search is node-capped; a capped search refuses the candidate instead of guessing.

Repair always terminates — split far enough and every cell is its own cage — so 40 of 40 seeds became boards at orders 5, 6 and 7. Nothing is rejected, so the bank's squares are exactly the sampler's squares:

grid	intercalates (40 generated boards)	exact mean
5×5	3.70 ± 0.17	3.57
6×6	8.18 ± 0.70	8.27
7×7	11.03 ± 0.81	10.53

All three are within two standard errors — compare 6.93 for keep-only-unique on 6×6 above.

The intercalate check runs before every search, and found 93 of 144 repairs (64.6%) without one. It did not make the generator faster: median 206 ms per 7×7 board with it and 224 ms without, worst case 26,448 ms against 12,514 ms, and search nodes about the same (20,566 against 18,622). The wall clock is set by a few slow boards that the node count does not explain. The check stays as a proof rather than a speed-up: it names the second answer, and its four cells are where the split goes.

What an operation sign is worth

Some puzzle books sell Calcudoku with the signs removed as the harder edition. That is a clean ablation: same cages, same targets, and a cage holds if any operation its size allows reaches the target. Hiding can add answers and never remove one; the test suite checks that, and checks that every rung is sound under both readings. (Last entry I measured an ablation whose own side was unsound, and the headline was a bug. The soundness test now loops over both modes.)

Cage by cage, the sign is worth little:

cage	cages	fillings with sign	without	bits the sign carries	no change at all
2 cells, ÷	350	3.51	9.23	1.455	0 (0.0%)
2 cells, ×	300	2.32	5.30	0.931	131 (43.7%)
2 cells, −	1,059	7.63	10.09	0.468	472 (44.6%)
3 cells, ×	624	7.20	10.61	0.367	458 (73.4%)
4 cells, ×	294	25.93	27.39	0.106	277 (94.2%)
5 cells, ×	118	87.54	89.97	0.044	117 (99.2%)

Of 4,041 multi-cell cages cut at random, 2,043 (50.6%) admit exactly the same fillings with the sign hidden. The most informative sign is division on a pair: hiding it multiplies the fillings by 2.74 (1.46 bits), and no division cage survived hiding untouched, because a quotient target is always reachable another way too. From four cells up, sum and product rarely collide and the sign carries at most 0.11 bits.

Board by board the loss is small too. Of the 16,000 random cuts, 7,296 were unique with signs and 6,711 without: hiding every sign costs 8.0% of the unique boards. The shipped bank was tightened with the signs shown, and still 19 of 20 boards have exactly one answer with every sign hidden (the exception, 7x7-20, opens up to 2). The page has a hide the signs switch.

The ladder, priced rung by rung

rung	cells settled, climbing	finished, climbing	finished without it	search nodes without it
cage	43 (6.8%)	0 / 20	20 / 20	—
single	122 (19.4%)	4 / 20	20 / 20	149 (with: 149)
hidden	257 (40.8%)	9 / 20	20 / 20	229 (with: 149)
must	411 (65.2%)	13 / 20	20 / 20	459 (with: 149)
subset	453 (71.9%)	14 / 20	20 / 20	—
probe	630 (100.0%)	20 / 20	14 / 20	—

must says: a value that every filling of a cage puts on one line is off-limits to the rest of that line. It shares its enumeration with the cage rung:

for (const [l, m] of s.need) {
  if (!m) continue;
  for (const c of lines(n)[l])
    if (b.p.cageOf[c] !== b.p.cageOf[k.cells[0]]) if (!restrict(b, c, ~m, ctx)) return false;
}

Two things stand out. First, single is a special case of must: a settled cell sits in some cage, every filling of that cage puts its number on its row and column, so must removes it outside the cage and cage removes it inside. Dropping single from the search leaves the node count exactly at 149. It stays on the ladder because it is the step a person takes first. Second, drop any one rung and the rest still finish every board without search — except probe, without which 14 of 20 finish. The search leans on must and hidden: 459 and 229 nodes without them.

I price rungs one at a time because a sound rung that never fires is invisible to end-to-end tests (that bit me on an earlier puzzle). This time it found not a dead rung, but one that another rung fully contains.

What I'd tell you to steal

1. How you read an MCMC chain matters as much as the chain. "Wait until proper" samples re-entry times, and ten times more steps does not fix it. Condition at a fixed time and discard the misses.
2. Hold a sampler to an exact distribution. Here I could enumerate to order 7. The bias only became visible with both the exact value and the noise floor of perfect draws.
3. Find the cheapest shape of a second answer first. In Calcudoku it is an intercalate flip, and it explains 76.5% of non-unique boards without a search.
4. Filtering bends the distribution; repairing does not. Keep-only-unique leans away from intercalates; a generator that splits cages instead never discards a square and carries the sampler's uniformity through.
5. Measure an ablation per piece and per board. The sign is worth zero bits on half of all cages, and hiding all of them costs 8.0% of unique boards.
6. Price rungs by removing them. single settles 79 cells when you climb, and removing it costs the search zero nodes.

Every number here is emitted by npm run stats, and the README and page prose are written from src/stats.json by npm run notes. 20 boards shipped, 23 tests.

Play it · Source