Domicules: tilings of the 2n × 2n torus (A401670)

const { drawTilings } = require('./site/notes.js');
const { svgs, total } = drawTilings({ rowCount: 8, colCount: 8, pieces: 'domicule', surface: 'torus', spread: 4 });

figure(svgs, `4 of the ${total} domicule tilings of the 8 × 8 torus.`);

const reference = require('./js/solvers/all-lattices-cylinder-torus-bigint.js');
const modular = require('./js/solvers/all-lattices-cylinder-torus-modular.js');
const known = require('./js/core/cases').testCasesOf('all-lattices-cylinder-torus-modular');
const stored = (side, twisted) =>
  known.find(t => t.lattice === 'king' && t.torus && !t.twisted === !twisted && t.rowCount === side && t.colCount === side);

nb.torus = {};
nb.twisted = {};
const rows = [];

for (let n = 1; n <= 8; n++) {
  const side = 2 * n;
  const board = { rowCount: side, colCount: side, lattice: 'king', torus: true };
  let count;
  let twisted;
  let method;

  if (n <= 6) {
    count = await modular.calculateTotalTilingCombinations(board);
    twisted = await modular.calculateTotalTilingCombinations({ ...board, twisted: true });
    method = 'modular';

    if (n <= 4) {
      const check = reference.calculateTotalTilingCombinations(board);
      const checkTwisted = reference.calculateTotalTilingCombinations({ ...board, twisted: true });

      if (check !== count || checkTwisted !== twisted) {
        throw new Error(`the two solvers differ on the ${side}×${side} torus`);
      }

      method = 'BigInt and modular';
    }
  } else if (stored(side, false)) {
    count = BigInt(stored(side, false).expectedResult);
    twisted = stored(side, true) ? BigInt(stored(side, true).expectedResult) : null;
    method = 'spec/cases/';
  } else {
    continue;
  }

  nb.torus[side] = count;
  nb.twisted[side] = twisted;
  rows.push([n, `${side}×${side}`, String(count), String(count).length, method]);
}

table(['n', 'board', 'a(n)', 'digits', 'solver'], rows);

const fs = require('fs');
const path = require('path');
const modular = require('./js/solvers/all-lattices-cylinder-torus-modular.js');
const root = nb.context.root;
const validate = fs.readFileSync(path.join(root, 'research/king-entropy/out/out_validate.txt'), 'utf-8');
const [, triangular] = validate.match(/Kasteleyn Bloch average, 64x64 grid: s_tri = ([0-9.]+)/);
const king = JSON.parse(fs.readFileSync(path.join(root, 'research/king-entropy/data/errbar.json'), 'utf-8')).recommended.s;
const entropy = { square: 0.29156090403081878, triangular: Number(triangular), king: Number(king) };
const { logCount } = nb.context;

nb.entropy = entropy;
nb.ratio = {};
const rows = [];

for (const side of Object.keys(nb.torus).map(Number)) {
  const ratios = {};

  for (const lattice of ['square', 'triangular', 'king']) {
    const count = lattice === 'king' ? nb.torus[side] : side <= 12 ? await modular.calculateTotalTilingCombinations({ rowCount: side, colCount: side, lattice, torus: true }) : null;

    ratios[lattice] = count === null ? null : Math.exp(logCount(count) - entropy[lattice] * side * side);
  }

  nb.ratio[side] = ratios;
  rows.push([
    `${side}×${side}`,
    ...['square', 'triangular', 'king'].map(lattice => (ratios[lattice] === null ? '' : ratios[lattice].toFixed(10))),
  ]);
}

const kingSides = Object.keys(nb.ratio).map(Number).filter(side => side >= 4);
const gaps = kingSides.map(side => nb.ratio[side].king - 2);

if (gaps.some((gap, i) => !(gap > 0) || (i > 0 && !(gap < gaps[i - 1])))) {
  throw new Error(`the king's graph ratio does not approach 2 from above: ${gaps.join(', ')}`);
}

table(['torus', 'square', 'triangular', 'king'], rows);

const lattices = [
  ['square', 'square lattice', 'var(--series-1)'],
  ['triangular', 'triangular', 'var(--series-2)'],
  ['king', "king's graph", 'var(--series-3)'],
];
const sides = Object.keys(nb.ratio).map(Number);

chart({
  title: 'The torus count divided by e^(sL²) against the side L, for three lattices',
  x: { label: 'side L of the L × L torus', min: 2, max: Math.max(...sides), ticks: sides },
  y: { label: 'Z / e^(sL²)', min: 1.9, max: 2.65, ticks: [2, 2.2, 2.4, 2.6] },
  reference: { y: 2 },
  series: lattices.map(([lattice, name, color]) => ({
    name,
    color,
    points: sides.filter(side => nb.ratio[side][lattice] !== null).map(side => [side, nb.ratio[side][lattice]]),
  })),
  describe: ([side, ratio]) => `${side}×${side} torus, Z/e^(sL²) = ${ratio.toFixed(10)}`,
});

const s = nb.entropy.king;
const { logCount } = nb.context;
const rows = Object.keys(nb.torus)
  .map(Number)
  .filter(side => side >= 4)
  .map((side, i, sides) => {
    const difference = (logCount(nb.torus[side]) - Math.LN2) / (side * side) - s;
    const before = i > 0 ? (logCount(nb.torus[sides[i - 1]]) - Math.LN2) / (sides[i - 1] * sides[i - 1]) - s : null;

    return [`${side}×${side}`, difference.toExponential(2), before === null ? '' : (before / difference).toFixed(2)];
  });

table(['torus', 'log(Z/2)/L² − s', 'factor per step'], rows);

const fs = require('fs');
const path = require('path');
const levels = fs.readFileSync(path.join(nb.context.root, 'research/king-entropy/log-convexity/out_levels.txt'), 'utf-8');
const [, xi] = levels.match(/RESULT xi = 1\/m_v = ([0-9.]+)/);
const sides = Object.keys(nb.twisted).map(Number).filter(side => nb.twisted[side] !== null);
const share = side => Number(nb.twisted[side] * 10n ** 18n / nb.torus[side]) / 1e18;

table(
  ['torus', 'signed sum', 'share of Z', 'factor per step'],
  sides.map((side, i) => [
    `${side}×${side}`,
    String(nb.twisted[side]),
    share(side).toExponential(2),
    i > 0 ? (share(sides[i - 1]) / share(side)).toFixed(2) : '',
  ]),
);
print(`e^(2/ξ) = ${Math.exp(2 / Number(xi)).toFixed(2)} with ξ = ${nb.context.roundDecimal(xi, 6)}`);

e^(2/ξ) = 5.04 with ξ = 1.236091

const { formatTime } = nb.context;

table(
  ['board', 'modular, Rust', 'modular, Python', 'modular, JavaScript', 'BigInt, Rust', 'BigInt, Python', 'BigInt, JavaScript'],
  nb.context.periodicBenchmark
    .filter(r => r.surface === 'torus')
    .map(r => [
      `${r.rowCount}×${r.colCount}`,
      ...[r.languages.rs, r.languages.py, r.modular, r.bigintLanguages.rs, r.bigintLanguages.py, r.bigint].map(formatTime),
    ]),
);
print(`Run: ${nb.context.benchmarkRun}`);

Run: benchmarks/results/ee071e57836d-node26.11.1-python3.14.5-rust1.99.1

print(require('fs').readFileSync(require('path').join(nb.context.root, 'oeis/A401670/A401670.gp'), 'utf-8').trim());

\\ Domicule tilings of the h X w torus = perfect matchings of the king's graph C_h X C_w (strong product);
\\ for h or w = 2 its multiple edges are counted as separate pieces; a(n) is OEIS A401670.
\\ R(w, c, f, t) lists the ways to finish a row whose filled cells are f from column c on, as the sets t
\\ of next-row cells they cover: a free cell is matched with its right neighbour (or, in column 0,
\\ with the last column across the seam), or with the cell below, below-left or below-right.
\\ M(w)[f+1, t+1] counts the transitions from f to t, and the torus count is the trace of M(w)^h.
R(w, c, f, t) = {
  my(v = []);
  if (c == w, return([t]));
  if (bittest(f, c), return(R(w, c+1, f, t)));
  if (c+1 < w && !bittest(f, c+1), v = concat(v, R(w, c+2, f + 2^c + 2^(c+1), t)));
  if (c == 0 && !bittest(f, w-1), v = concat(v, R(w, 1, f + 1 + 2^(w-1), t)));
  for (d = -1, 1, my(j = (c+d) % w); if (!bittest(t, j), v = concat(v, R(w, c+1, f + 2^c, t + 2^j))));
  v
};
M(w) = {
  my(A = matrix(2^w, 2^w));
  for (f = 0, 2^w-1, my(v = R(w, 0, f, 0)); for (i = 1, #v, A[f+1, v[i]+1]++));
  A
};
T(h, w) = trace(M(w)^h);
a(n) = trace(M(2*n)^(2*n));