Diagonal pieces: domicule tilings of the 2n × 2n board by their count (A401673)

const { drawTilings } = require('./site/notes.js');
const counts = [0, 8, 16, 32];
const svgs = counts.flatMap(k => drawTilings({ rowCount: 8, colCount: 8, pieces: 'domicule', diagonals: k, ranks: [0] }).svgs);

figure(svgs, 'Tilings of the 8 × 8 board with k = 0, 8, 16 and 32 diagonal pieces.', { highlight: true });

const reference = require('./js/solvers/domicule-rectangle-polynomial-bigint.js');
const data = require('./js/core/cases');
const stored = side => data.testCasesOf('domicule-rectangle-polynomial-modular').find(t => t.rowCount === side && t.colCount === side);
const domino = side => data.testCasesOf('domino-rectangle-formula').find(t => t.rowCount === side && t.colCount === side).expectedResult;

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

for (let n = 0; n <= 10; n++) {
  const side = 2 * n;
  const board = { rowCount: side, colCount: side };
  let coefficients;
  let source;

  if (n <= 6) {
    coefficients = reference.calculateDiagonalPolynomial(board);
    source = 'BigInt';
  } else if (stored(side)) {
    coefficients = stored(side).expectedResult.split(',').map(BigInt);
    source = 'spec/cases/';
  } else {
    continue;
  }

  const total = coefficients.reduce((sum, c) => sum + c, 0n);
  const row = coefficients.filter((_, k) => k % 2 === 0);
  const problems = [
    total !== (await nb.context.countDomicule(board)) && 'the sum is not A239273',
    String(coefficients[0]) !== domino(side) && 'T(n, 0) is not A004003',
    coefficients.some((c, k) => k % 2 === 1 && c !== 0n) && 'an odd number of diagonal pieces',
    row.at(-1) !== 1n && `T(n, n²) = ${row.at(-1)}`,
    n > 0 && row.at(-2) !== BigInt(n * n * (n + 1)) && `T(n, n² − 1) = ${row.at(-2)}`,
    source === 'BigInt' && stored(side) && stored(side).expectedResult !== coefficients.join(',') && 'the BigInt row differs from spec/cases/',
  ].filter(Boolean);

  if (problems.length) {
    throw new Error(`n = ${n}: ${problems.join(', ')}`);
  }

  const peak = row.indexOf(row.reduce((max, c) => (c > max ? c : max), 0n));

  nb.rows[n] = row;
  rows.push([n, `${side}×${side}`, row.length, peak, String(row[peak]).length, source]);
}

table(['n', 'board', 'terms T(n, 0..n²)', 'j of the largest term', 'digits of the largest term', 'solver'], rows);
print('The first rows:');

for (let n = 0; n <= 3; n++) {
  print(`n = ${n}: ${nb.rows[n].join(', ')}`);
}

The first rows:
n = 0: 1
n = 1: 2, 1
n = 2: 36, 126, 105, 12, 1
n = 3: 6728, 84132, 413932, 934213, 1000638, 491689, 99554, 6638, 36, 1

const rows = Object.entries(nb.rows)
  .filter(([n]) => Number(n) > 0)
  .map(([n, row]) => {
    const total = row.reduce((sum, c) => sum + c, 0n);
    const scale = 10n ** 15n;
    const mean = Number((row.reduce((sum, c, j) => sum + 2n * BigInt(j) * c, 0n) * scale) / total) / Number(scale);
    const square = Number((row.reduce((sum, c, j) => sum + 4n * BigInt(j * j) * c, 0n) * scale) / total) / Number(scale);
    const pieces = 2 * Number(n) * Number(n);

    return [n, `${2 * n}×${2 * n}`, (mean / pieces).toFixed(6), ((square - mean * mean) / (2 * pieces)).toFixed(6)];
  });

table(['n', 'board', 'share of diagonal pieces', 'variance per site'], rows);

const weights = [[0, 1], [1, 2], [1, 1], [3, 2], [2, 1], [3, 1]];
const { logCount } = nb.context;
const evaluate = (row, [p, q]) => {
  const degree = row.length - 1;

  return row.reduce((sum, c, j) => sum + c * BigInt(p) ** BigInt(2 * j) * BigInt(q) ** BigInt(2 * (degree - j)), 0n);
};

table(
  ['n', 'board', ...weights.map(([p, q]) => `x = ${q === 1 ? p : `${p}/${q}`}`)],
  Object.entries(nb.rows)
    .filter(([n]) => Number(n) >= 3)
    .map(([n, row]) => {
      const cells = 4 * Number(n) * Number(n);

      return [
        n,
        `${2 * n}×${2 * n}`,
        ...weights.map(([p, q]) => ((logCount(evaluate(row, [p, q])) - 2 * (row.length - 1) * Math.log(q)) / cells).toFixed(6)),
      ];
    }),
);

const freeEnergy = (row, x) => {
  const terms = row.map((c, j) => (c === 0n ? -Infinity : logCount(c) + (j === 0 ? 0 : 2 * j * Math.log(x))));
  const top = Math.max(...terms);

  return top + Math.log(terms.reduce((sum, t) => sum + Math.exp(t - top), 0));
};
const shown = [4, 6, 8, 10].filter(n => nb.rows[n]);
const xs = Array.from({ length: 61 }, (_, i) => i / 20);

chart({
  title: 'log P(x) / (4n²) against the weight x of a diagonal piece, for boards 2n × 2n',
  x: { label: 'weight x of a diagonal piece', min: 0, max: 3, ticks: [0, 0.5, 1, 1.5, 2, 2.5, 3] },
  y: { label: 'log P(x) / (4n²)', min: 0.2, max: 0.9, ticks: [0.2, 0.4, 0.6, 0.8] },
  series: shown.map((n, i) => ({
    name: `${2 * n}×${2 * n}`,
    color: `var(--ordinal-${i + 1})`,
    points: xs.map(x => [x, freeEnergy(nb.rows[n], x) / (4 * n * n)]),
    markers: weights.map(([p, q]) => p / q).map(x => [x, freeEnergy(nb.rows[n], x) / (4 * n * n)]),
  })),
  describe: ([x, value]) => `x = ${x}, log P(x)/(4n²) = ${value.toFixed(6)}`,
});

const { formatTime } = nb.context;

table(
  ['board', 'modular, Rust', 'modular, Python', 'modular, JavaScript', 'BigInt, Rust', 'BigInt, Python', 'BigInt, JavaScript'],
  nb.context.diagonalBenchmark.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/A401673/A401673.gp'), 'utf-8').trim());

\\ Domicule tilings of an h X w board by the number of diagonal pieces (OEIS A401673): P(h, w) is a polynomial
\\ in x, whose coefficient of x^k counts the tilings with k pieces joining cells that share only a corner.
\\ The dynamic programming of A239264.gp, with a piece to the down-left or down-right neighbour weighted by x.
\\ The number of diagonal pieces is even, so row(n) lists the tilings of the 2n X 2n board with 0, 2, ..., 2n^2 of them.
P(h, w) = {
  my(S = 2^(w+2), v = vector(S), u);
  if (h*w == 0, return(1));
  v[1] = 1;
  for (p = 0, h*w-1,
    my(r = p\w, c = p%w, d = [], e = []);
    if (c+1 < w, d = concat(d, 1); e = concat(e, 1));
    if (r+1 < h,
      if (c > 0, d = concat(d, w-1); e = concat(e, x));
      d = concat(d, w); e = concat(e, 1);
      if (c+1 < w, d = concat(d, w+1); e = concat(e, x)));
    u = vector(S);
    for (s = 0, S-1,
      my(y = v[s+1]);
      if (y,
        if (bitand(s, 1),
          u[(s>>1) + 1] += y,
          for (i = 1, #d, if (!bittest(s, d[i]), my(t = (s + 2^d[i])>>1); u[t+1] += e[i]*y)))));
    v = u);
  v[1]
};
row(n) = my(q = P(2*n, 2*n)); vector(n^2+1, k, polcoef(q, 2*(k-1)));