Hyperbolic Tiling of the Poincaré Disk

// Poincaré disk tiling by {p, q} regular polygons.
// We build the fundamental p-gon centered at origin, then recursively
// reflect it across each side using hyperbolic isometries represented
// as Möbius transformations on complex numbers.

const P = 7;       // sides per polygon (heptagon)
const Q = 3;       // polygons meeting at each vertex
const MAX_DEPTH = 5;

// Complex helpers (z = [re, im])
function cadd(a, b) { return [a[0] + b[0], a[1] + b[1]]; }
function csub(a, b) { return [a[0] - b[0], a[1] - b[1]]; }
function cmul(a, b) { return [a[0]*b[0] - a[1]*b[1], a[0]*b[1] + a[1]*b[0]]; }
function cdiv(a, b) { const d = b[0]*b[0] + b[1]*b[1]; return [(a[0]*b[0]+a[1]*b[1])/d, (a[1]*b[0]-a[0]*b[1])/d]; }
function cconj(a) { return [a[0], -a[1]]; }
function cscale(a, s) { return [a[0]*s, a[1]*s]; }

// A Möbius transformation M(z) = (a z + b) / (c z + d), stored as [a,b,c,d].
function mapply(M, z) {
  return cdiv(cadd(cmul(M[0], z), M[1]), cadd(cmul(M[2], z), M[3]));
}
function mcompose(A, B) {
  // (A ∘ B)(z): coefficients are matrix product
  return [
    cadd(cmul(A[0], B[0]), cmul(A[1], B[2])),
    cadd(cmul(A[0], B[1]), cmul(A[1], B[3])),
    cadd(cmul(A[2], B[0]), cmul(A[3], B[2])),
    cadd(cmul(A[2], B[1]), cmul(A[3], B[3])),
  ];
}

// Reflection across a generalized circle (geodesic in the disk) is
// anti-holomorphic. We store each side reflection as a Möbius M plus
// a "conjugate first" flag: z -> M(conj(z)). Composition of two such
// flips Q's: holomorphic = even count of conjugations.
//
// To keep things uniform we instead build the SYMMETRY GROUP from
// rotations (holomorphic) only: for each side we use the rotation
// that maps the polygon to its neighbor across that side. That
// rotation is: reflect across the side, then reflect across the line
// through origin perpendicular to the side. Both reflections compose
// to a holomorphic Möbius transformation — a hyperbolic translation.

let sideGens = [];     // P generators (Möbius), one per side
let basePoly = [];     // vertices of fundamental polygon
let center0 = [0, 0];
let rotState = 0;

function buildBasePolygon() {
  // Vertex distance from origin (Euclidean) for regular {p,q} tiling:
  // r^2 = (cos(π/p + π/q)) / (cos(π/p - π/q))
  const a = Math.PI / P, b = Math.PI / Q;
  const num = Math.cos(a + b);
  const den = Math.cos(a - b);
  const r2 = num / den;
  const r = Math.sqrt(Math.max(0, r2));
  basePoly = [];
  for (let k = 0; k < P; k++) {
    const t = (2 * Math.PI * k) / P;
    basePoly.push([r * Math.cos(t), r * Math.sin(t)]);
  }
}

function buildGenerators() {
  // Generator for side k: hyperbolic translation along the perpendicular
  // from origin to that side's midpoint. We construct it as
  //   T_k = R(θ_k) ∘ T_d ∘ R(-θ_k)
  // where T_d is a translation along the x-axis by hyperbolic distance 2d,
  // and θ_k is the angle of side k's midpoint.
  //
  // T_d as a Möbius on the disk: z -> (z + t) / (1 + t z),  t = tanh(d).
  const a = Math.PI / P, b = Math.PI / Q;
  // Distance from origin to side midpoint (Euclidean):
  // tanh(d) = cos(π/p) * sqrt((1 - tan²(π/q) tan²(π/p)) ...).
  // Simpler closed form: the midpoint sits at Euclidean radius
  //   m = r * cos(π/p)   where r is the vertex radius above? No — for
  // a regular hyperbolic polygon centered at 0, the side midpoint is at
  //   tanh(d_m) where  cosh(d_m) = cos(π/q) / sin(π/p).
  const coshDm = Math.cos(b) / Math.sin(a);
  const sinhDm = Math.sqrt(coshDm * coshDm - 1);
  const tanhDm = sinhDm / coshDm;
  // The translation we want maps the polygon to its mirror across the side,
  // i.e. translation distance 2*d_m along the perpendicular.
  // tanh(2 d_m) = 2 t / (1 + t²) with t = tanh(d_m).
  const t = tanhDm;
  const t2 = (2 * t) / (1 + t * t);

  // T_{t2} along x-axis on the disk:
  //   z -> (z + t2) / (t2 z + 1)
  // Möbius coefficients: a=1, b=t2, c=t2, d=1.
  const Tx = [[1, 0], [t2, 0], [t2, 0], [1, 0]];

  sideGens = [];
  for (let k = 0; k < P; k++) {
    // Side k lies between vertex k and vertex k+1. Its midpoint angle is
    //   θ_k = (2k+1) π / P
    const th = ((2 * k + 1) * Math.PI) / P;
    const c = Math.cos(th), s = Math.sin(th);
    // Rotation R(θ) as Möbius: z -> e^{iθ} z
    //   a = e^{iθ}, b = 0, c = 0, d = 1
    const Rp = [[c, s], [0, 0], [0, 0], [1, 0]];
    const Rm = [[c, -s], [0, 0], [0, 0], [1, 0]];
    sideGens.push(mcompose(Rp, mcompose(Tx, Rm)));
  }
}

const IDENT = [[1, 0], [0, 0], [0, 0], [1, 0]];

// Recursive tile list: each entry is { M, depth, fromSide }
function enumerateTiles() {
  const tiles = [{ M: IDENT, depth: 0, fromSide: -1 }];
  let frontier = [{ M: IDENT, depth: 0, fromSide: -1 }];
  for (let d = 0; d < MAX_DEPTH; d++) {
    const next = [];
    for (let i = 0; i < frontier.length; i++) {
      const f = frontier[i];
      for (let k = 0; k < P; k++) {
        // Don't immediately reverse the side we just came in through.
        if (k === f.fromSide) continue;
        const M2 = mcompose(f.M, sideGens[k]);
        // The "side we came in through" in the new tile is the side that
        // pairs with k under the generator. For a {p,q} tiling with even
        // p the pairing is k <-> k. For odd p (like 7) the back-side is
        // simply the same index in the new tile's local numbering — the
        // generator T_k swaps interior and exterior across side k. So:
        next.push({ M: M2, depth: d + 1, fromSide: k });
        tiles.push({ M: M2, depth: d + 1, fromSide: k });
      }
    }
    frontier = next;
    if (frontier.length > 4000) break; // safety cap
  }
  return tiles;
}

let tiles = [];
let R = 0;
let cx = 0, cy = 0;

function init(ctx) {
  buildBasePolygon();
  buildGenerators();
  tiles = enumerateTiles();
  rotState = 0;
}

function diskToScreen(z, width, height) {
  return [cx + z[0] * R, cy - z[1] * R];
}

function drawTile(g, verts, depth, width, height) {
  // Color by depth: cycle hue, fade with depth.
  const hue = (depth * 47 + rotState * 30) % 360;
  const lum = 62 - depth * 6;
  g.fillStyle = `hsl(${hue}, 70%, ${Math.max(18, lum)}%)`;
  g.strokeStyle = `hsla(${hue}, 80%, 80%, 0.55)`;
  g.lineWidth = 0.6;
  g.beginPath();
  for (let i = 0; i < verts.length; i++) {
    const [x, y] = diskToScreen(verts[i], width, height);
    if (i === 0) g.moveTo(x, y); else g.lineTo(x, y);
  }
  g.closePath();
  g.fill();
  g.stroke();
}

function tick({ dt, ctx, width, height }) {
  cx = width / 2;
  cy = height / 2;
  R = Math.min(width, height) * 0.48;
  rotState += dt * 0.08;

  // Apply a slow rotation by composing every tile transform with a
  // rotation in front. We just transform the base polygon vertices.
  const cR = Math.cos(rotState), sR = Math.sin(rotState);
  const Rrot = [[cR, sR], [0, 0], [0, 0], [1, 0]];

  // Background — Poincaré disk
  ctx.fillStyle = '#06080d';
  ctx.fillRect(0, 0, width, height);
  ctx.save();
  ctx.beginPath();
  ctx.arc(cx, cy, R, 0, Math.PI * 2);
  ctx.fillStyle = '#0b1018';
  ctx.fill();
  ctx.clip();

  // Sort tiles back-to-front: deeper tiles are near the boundary,
  // shallower near the center. Draw deepest first so center sits on top.
  // tiles is already roughly BFS order, so iterate in reverse.
  for (let i = tiles.length - 1; i >= 0; i--) {
    const t = tiles[i];
    const M = mcompose(Rrot, t.M);
    const verts = new Array(P);
    let visible = false;
    for (let v = 0; v < P; v++) {
      const z = mapply(M, basePoly[v]);
      verts[v] = z;
      if (z[0] * z[0] + z[1] * z[1] < 0.9999) visible = true;
    }
    if (!visible) continue;
    drawTile(ctx, verts, t.depth, width, height);
  }

  ctx.restore();

  // Boundary circle
  ctx.strokeStyle = 'rgba(180, 200, 255, 0.35)';
  ctx.lineWidth = 1;
  ctx.beginPath();
  ctx.arc(cx, cy, R, 0, Math.PI * 2);
  ctx.stroke();

  // HUD
  ctx.fillStyle = 'rgba(220, 230, 255, 0.85)';
  ctx.font = '12px system-ui, sans-serif';
  ctx.fillText(`{${P}, ${Q}} tiling   depth ${MAX_DEPTH}   tiles ${tiles.length}`, 10, 18);
}