31:T3023,viewof apollonian_viewer = { const core = apollonian_core; // ---- canvas size ---- let canvasW, canvasH; if (width <= 720) { canvasW = width; canvasH = width; } else { canvasW = 0.7 * width; canvasH = 0.7 * width; } // ---- UI ---- const root = document.createElement("div"); root.style.cssText = ` font-family: ui-sans-serif, system-ui, -apple-system, Segoe UI, Roboto, Helvetica, Arial; width: 100%; `; const controls = document.createElement("div"); controls.style.cssText = `display:flex; gap:12px; align-items:center; flex-wrap:wrap; padding: 8px 0 12px 0;`; const label = document.createElement("div"); label.textContent = "Example:"; label.style.cssText = `font-weight:600;`; const sel = document.createElement("select"); sel.style.cssText = `padding:6px 10px; border:1px solid #ccc; border-radius:8px; background:#fff;`; core.exampleNames.forEach((name, i) => { const opt = document.createElement("option"); opt.value = String(i); opt.textContent = `${i + 1}. ${name}`; sel.appendChild(opt); }); const stats = document.createElement("div"); stats.style.cssText = `opacity:0.75; font-size:12px; white-space:nowrap;`; controls.append(label, sel, stats); root.appendChild(controls); const canvas = document.createElement("canvas"); canvas.width = canvasW; canvas.height = canvasH; canvas.style.cssText = `border:1px solid #ddd; border-radius:12px; background:#fff; display:block;`; root.appendChild(canvas); const ctx = canvas.getContext("2d"); ctx.textAlign = "center"; ctx.textBaseline = "middle"; // ---- projector ---- function makeProjector(view) { const { Rat, rAdd, rSub, rMul, rInv, rToNumber, rAbs } = core; const { xminR, yminR, ddR } = view; const invDD = rInv(ddR); const ymaxR = rAdd(yminR, ddR); const pxN = (xR) => rToNumber(rMul(rMul(rSub(xR, xminR), invDD), Rat(BigInt(canvasW), 1n))); const pyN = (yR) => rToNumber(rMul(rMul(rSub(ymaxR, yR), invDD), Rat(BigInt(canvasH), 1n))); const rPixN = (kBig) => { const rR = rAbs(core.Rat(1n, kBig)); const rp = rMul(rMul(rR, invDD), core.Rat(BigInt(canvasW), 1n)); return rToNumber(rp); }; return { pxN, pyN, rPixN }; } // ---- label that fits inside circle (using actual measureText) ---- function drawLabelIfFits(kStr, cx, cy, rp) { // if the circle is small, don't label if (rp < 18) return; // target box: a bit smaller than diameter const maxW = 1.6 * rp; const maxH = 1.0 * rp; // start guess (depends on rp) let s = Math.min(22, 0.22 * rp); s = Math.max(7, s); // shrink to fit using measureText (few iterations) for (let iter = 0; iter < 6; iter++) { ctx.font = `${Math.floor(s)}px sans-serif`; const w = ctx.measureText(kStr).width; const h = s; // approx ascent+descent if (w <= maxW && h <= maxH) break; s *= 0.85; if (s < 7) return; } ctx.font = `${Math.floor(s)}px sans-serif`; ctx.fillStyle = "#000"; ctx.fillText(kStr, cx, cy); } // ---- numerically stable drawing of huge-radius circles ---- // For very large rp, ctx.arc(cx,cy,rp) can drift due to floating precision. // We draw a local polyline near the closest point to viewport center. function strokeCircleStable(cx, cy, rp) { // threshold: above this, arc can become unstable and slow const BIG = 20000; // px if (rp <= BIG && Number.isFinite(cx) && Number.isFinite(cy) && Number.isFinite(rp)) { ctx.beginPath(); ctx.moveTo(cx + rp, cy); ctx.arc(cx, cy, rp, 0, 2 * Math.PI); ctx.stroke(); ctx.closePath(); return; } // viewport center const vx = canvasW / 2; const vy = canvasH / 2; // vector from center to viewport center let dx = vx - cx; let dy = vy - cy; let d = Math.hypot(dx, dy); if (!(d > 0)) d = 1; // unit radial towards viewport center const ux = dx / d; const uy = dy / d; // unit tangent (perp) const tx = -uy; const ty = ux; // closest point on circle to viewport center const px = cx + rp * ux; const py = cy + rp * uy; // choose angular span to cover canvas diagonal, thetaMax ~ diag/(2r) const diag = Math.hypot(canvasW, canvasH); const thetaMax = Math.min(Math.PI / 2, diag / (2 * rp)); // small for big rp // number of segments const segs = Math.max(8, Math.min(80, Math.ceil(thetaMax * 60))); ctx.beginPath(); for (let i = -segs; i <= segs; i++) { const t = (i / segs) * thetaMax; // local expansion around closest point: // point(t) = p + r*(cos t - 1)*u + r*sin t * tangent const ct1 = Math.cos(t) - 1; const st = Math.sin(t); const x = px + rp * (ct1 * ux + st * tx); const y = py + rp * (ct1 * uy + st * ty); if (i === -segs) ctx.moveTo(x, y); else ctx.lineTo(x, y); } ctx.stroke(); ctx.closePath(); } // ---- generation knobs ---- // Deep zoom goal: keep minPixRadius small (0.2–0.5) // Speed goal: increase minPixRadius and/or lower maxCircles. const GEN = { maxCircles: 60000, minPixRadius: 0.35 }; // ---- renderer (keeps your zoom flow: zoom redraw; end regenerate) ---- function drawQuad(quad, opts = {}) { const { bounds = "auto", xmin, ymin, dd, interactive = true, show_curvatures = true, show_more = true } = opts; ctx.clearRect(0, 0, canvasW, canvasH); const [c0, c1, c2, c3] = quad; // view window let view; if (bounds === "auto") { const r0 = core.rDiv(core.Rat(-1n, 1n), core.Rat(c0[2], 1n)); const R = core.rMul(core.rAbs(r0), core.ratApprox(1.1, 10_000_000_000n)); view = { xminR: core.rNeg(R), yminR: core.rNeg(R), ddR: core.rMul(R, core.Rat(2n, 1n)) }; } else { view = { xminR: xmin, yminR: ymin, ddR: dd }; } const proj = makeProjector(view); const rect = core.rectFromView(view); const Bnum = core.computeB(c0, view.ddR); const all_circles = show_more ? core.gasket(c0, c1, c2, c3, rect, Bnum, { ...GEN, proj }) : [c0, c1, c2, c3]; // Cache visible circles in pixel coords const pixelAll = []; const pixelRest = []; for (let idx = 0; idx < all_circles.length; idx++) { const [xR, yR, k] = all_circles[idx]; const cx = proj.pxN(xR); const cy = proj.pyN(yR); const rp = proj.rPixN(k); // cull offscreen cheaply if (cx + rp < 0 || cx - rp > canvasW || cy + rp < 0 || cy - rp > canvasH) continue; const rec = { cx, cy, rp, k }; pixelAll.push(rec); if (idx >= 4) pixelRest.push(rec); } // draw most circles (stable) for (const { cx, cy, rp, k } of pixelRest) { strokeCircleStable(cx, cy, rp); if (show_curvatures) drawLabelIfFits(String(k), cx, cy, rp); } // shade first three for (const [xR, yR, k] of [c1, c2, c3]) { const cx = proj.pxN(xR); const cy = proj.pyN(yR); const rp = proj.rPixN(k); // fill as arc when reasonable; otherwise fill a polygon is overkill, so we skip fill for huge radii if (rp <= 20000) { ctx.beginPath(); ctx.moveTo(cx + rp, cy); ctx.arc(cx, cy, rp, 0, 2 * Math.PI); ctx.fillStyle = "#dddddd"; ctx.fill(); ctx.stroke(); ctx.closePath(); } else { // huge circle: stroke only (fill not meaningful visually) ctx.save(); ctx.strokeStyle = "#000"; strokeCircleStable(cx, cy, rp); ctx.restore(); } if (show_curvatures) drawLabelIfFits(String(k), cx, cy, rp); } // exterior { const cx = proj.pxN(c0[0]); const cy = proj.pyN(c0[1]); const rp = proj.rPixN(c0[2]); ctx.save(); ctx.lineWidth = 2; strokeCircleStable(cx, cy, rp); ctx.restore(); // optional label if (show_curvatures) { ctx.font = `36px sans-serif`; ctx.fillStyle = "#000"; ctx.fillText(String(c0[2]), 0.12 * canvasW, 0.12 * canvasH); } } stats.textContent = `circles: ${all_circles.length} (drawn: ${pixelAll.length}) minPixRadius: ${GEN.minPixRadius}`; // zoom: redraw cached pixel circles; end: compute new exact window and regenerate if (interactive) { d3.select(canvas).call( d3.zoom() .scaleExtent([0, 8]) .on("zoom", function () { ctx.save(); ctx.clearRect(0, 0, canvasW, canvasH); for (const { cx, cy, rp } of pixelAll) { const [tcx, tcy] = d3.event.transform.apply([cx, cy]); const [tcxr] = d3.event.transform.apply([cx + rp, cy]); const trp = Math.abs(tcxr - tcx); // during zoom, use plain arc for speed when possible, otherwise a quick polyline if (trp <= 20000) { ctx.beginPath(); ctx.moveTo(tcx + trp, tcy); ctx.arc(tcx, tcy, trp, 0, 2 * Math.PI); ctx.stroke(); ctx.closePath(); } else { // approximate huge circles as stable polyline in zoomed space // reuse the stable stroker by temporarily drawing in transformed coords: ctx.save(); // We can’t easily reuse strokeCircleStable with transformed coords cleanly, // so do a coarse local polyline here: const vx = canvasW / 2, vy = canvasH / 2; let dx = vx - tcx, dy = vy - tcy; let d = Math.hypot(dx, dy); if (!(d > 0)) d = 1; const ux = dx / d, uy = dy / d; const tx = -uy, ty = ux; const px = tcx + trp * ux, py = tcy + trp * uy; const diag = Math.hypot(canvasW, canvasH); const thetaMax = Math.min(Math.PI / 2, diag / (2 * trp)); const segs = Math.max(8, Math.min(60, Math.ceil(thetaMax * 50))); ctx.beginPath(); for (let i = -segs; i <= segs; i++) { const t = (i / segs) * thetaMax; const ct1 = Math.cos(t) - 1; const st = Math.sin(t); const x = px + trp * (ct1 * ux + st * tx); const y = py + trp * (ct1 * uy + st * ty); if (i === -segs) ctx.moveTo(x, y); else ctx.lineTo(x, y); } ctx.stroke(); ctx.closePath(); ctx.restore(); } } ctx.restore(); }) .on("end", function () { const t = d3.event.transform; const inv = t.invert.bind(t); const [pBLx, pBLy] = inv([0, canvasH]); const [pTRx, pTRy] = inv([canvasW, 0]); const fx0 = core.rDiv(core.rFromNumber(pBLx), core.Rat(BigInt(canvasW), 1n)); const fy0 = core.rDiv(core.rFromNumber(canvasH - pBLy), core.Rat(BigInt(canvasH), 1n)); const fx1 = core.rDiv(core.rFromNumber(pTRx), core.Rat(BigInt(canvasW), 1n)); const new_xmin = core.rAdd(view.xminR, core.rMul(view.ddR, fx0)); const new_ymin = core.rAdd(view.yminR, core.rMul(view.ddR, fy0)); const new_xmax = core.rAdd(view.xminR, core.rMul(view.ddR, fx1)); const new_dd = core.rSub(new_xmax, new_xmin); drawQuad(quad, { bounds: "set", xmin: new_xmin, ymin: new_ymin, dd: new_dd, interactive: true, show_curvatures, show_more }); d3.zoom().transform(d3.select(canvas), d3.zoomIdentity); }) ); } return { view, all_circles, pixelAll }; } // ---- state ---- let idx = 0; let quad = core.examples[idx]; function redrawAuto() { drawQuad(quad, { bounds: "auto", interactive: true, show_curvatures: true, show_more: true }); } sel.addEventListener("change", () => { idx = +sel.value; quad = core.examples[idx]; redrawAuto(); root.value = idx; root.dispatchEvent(new CustomEvent("input")); }); redrawAuto(); root.value = idx; root.dispatchEvent(new CustomEvent("input")); return root; } 32:T294c,apollonian_core = { // ------------------------- // BigInt "lazy" rationals // (no gcd reduction in hot path) // ------------------------- function Rat(n, d = 1n) { if (d === 0n) throw new Error("Rat: zero denominator"); if (d < 0n) { n = -n; d = -d; } return { n, d }; } const R0 = Rat(0n, 1n); const rNeg = (a) => Rat(-a.n, a.d); const rAdd = (a, b) => Rat(a.n * b.d + b.n * a.d, a.d * b.d); const rSub = (a, b) => Rat(a.n * b.d - b.n * a.d, a.d * b.d); const rMul = (a, b) => Rat(a.n * b.n, a.d * b.d); const rDiv = (a, b) => Rat(a.n * b.d, a.d * b.n); const rAbs = (a) => Rat(a.n < 0n ? -a.n : a.n, a.d); const rInv = (a) => Rat(a.d, a.n); const rToNumber = (a) => Number(a.n) / Number(a.d); function ratApprox(x, maxDen = 1_000_000_000n) { if (!Number.isFinite(x)) throw new Error("ratApprox: non-finite"); if (x === 0) return R0; const sign = x < 0 ? -1 : 1; x = Math.abs(x); let a0 = Math.floor(x); let p0 = 1n, q0 = 0n; let p1 = BigInt(a0), q1 = 1n; let frac = x - a0; let it = 0; while (frac > 0 && it++ < 80) { const inv = 1 / frac; const ai = Math.floor(inv); const p2 = BigInt(ai) * p1 + p0; const q2 = BigInt(ai) * q1 + q0; if (q2 > maxDen) break; p0 = p1; q0 = q1; p1 = p2; q1 = q2; frac = inv - ai; } return Rat(sign < 0 ? -p1 : p1, q1); } const rFromNumber = (x) => (Number.isInteger(x) ? Rat(BigInt(x), 1n) : ratApprox(x, 10_000_000_000n)); function rCmp(a, b) { const L = a.n * b.d; const R = b.n * a.d; return L < R ? -1 : (L > R ? 1 : 0); } // ------------------------- // Circles: [xRat, yRat, kBigInt] // ------------------------- const Circle = (x, y, k) => [x, y, k]; function toBigCircle([x, y, k]) { const xR = (x?.n !== undefined) ? x : (typeof x === "bigint" ? Rat(x, 1n) : rFromNumber(x)); const yR = (y?.n !== undefined) ? y : (typeof y === "bigint" ? Rat(y, 1n) : rFromNumber(y)); const kB = (typeof k === "bigint") ? k : BigInt(Math.trunc(k)); return Circle(xR, yR, kB); } // your Descartes-based update, exact function circle(c0, c1, c2, c3) { const b0 = c0[2], b1 = c1[2], b2 = c2[2], b3 = c3[2]; const b = 2n * b0 + 2n * b1 + 2n * b2 - b3; const numX = rSub( rAdd( rAdd(rMul(Rat(2n * b0, 1n), c0[0]), rMul(Rat(2n * b1, 1n), c1[0])), rMul(Rat(2n * b2, 1n), c2[0]) ), rMul(Rat(b3, 1n), c3[0]) ); const numY = rSub( rAdd( rAdd(rMul(Rat(2n * b0, 1n), c0[1]), rMul(Rat(2n * b1, 1n), c1[1])), rMul(Rat(2n * b2, 1n), c2[1]) ), rMul(Rat(b3, 1n), c3[1]) ); return Circle(rDiv(numX, Rat(b, 1n)), rDiv(numY, Rat(b, 1n)), b); } function tangency(c1, c2) { const k1 = c1[2], k2 = c2[2]; const r1 = Rat(1n, k1); const r2 = Rat(1n, k2); const d = rAdd(r1, r2); const x = rDiv(rAdd(rMul(c1[0], r2), rMul(c2[0], r1)), d); const y = rDiv(rAdd(rMul(c1[1], r2), rMul(c2[1], r1)), d); return [x, y]; } function inner_triangle(c1, c2, c3) { return [tangency(c1, c2), tangency(c2, c3), tangency(c3, c1)]; } // ------------------------- // Pruning geometry // We make pruning cheap: numeric first, exact only if needed. // ------------------------- function leftOf(A, B, P) { const [x0, y0] = A, [x1, y1] = B, [x, y] = P; // x1*(y-y0) + x*(y0-y1) + x0*(y1-y) const s = rAdd( rAdd(rMul(x1, rSub(y, y0)), rMul(x, rSub(y0, y1))), rMul(x0, rSub(y1, y)) ); return s.n > 0n; } function sameSide(A, B, P, Q) { const ab = leftOf(A, B, P); const cd = leftOf(A, B, Q); return (ab && cd) || (!ab && !cd); } function some_sameSide(A, B, P, rect) { for (let i = 0; i < rect.length; i++) { if (sameSide(A, B, P, rect[i])) return true; } return false; } // numeric bbox helpers for fast reject function bboxPolyN(poly) { let xmin = Infinity, xmax = -Infinity, ymin = Infinity, ymax = -Infinity; for (const p of poly) { const x = rToNumber(p[0]), y = rToNumber(p[1]); if (x < xmin) xmin = x; if (x > xmax) xmax = x; if (y < ymin) ymin = y; if (y > ymax) ymax = y; } return [xmin, ymin, xmax, ymax]; } function bboxOverlapN(a, b) { return !(a[2] < b[0] || b[2] < a[0] || a[3] < b[1] || b[3] < a[1]); } // exact triangle-rect intersection (used rarely; numeric pretest first) function orient(a, b, c) { const bax = rSub(b[0], a[0]); const bay = rSub(b[1], a[1]); const cax = rSub(c[0], a[0]); const cay = rSub(c[1], a[1]); return rSub(rMul(bax, cay), rMul(bay, cax)); } function onSegment(a, b, p) { const minX = rCmp(a[0], b[0]) <= 0 ? a[0] : b[0]; const maxX = rCmp(a[0], b[0]) >= 0 ? a[0] : b[0]; const minY = rCmp(a[1], b[1]) <= 0 ? a[1] : b[1]; const maxY = rCmp(a[1], b[1]) >= 0 ? a[1] : b[1]; return rCmp(minX, p[0]) <= 0 && rCmp(p[0], maxX) <= 0 && rCmp(minY, p[1]) <= 0 && rCmp(p[1], maxY) <= 0; } function segIntersects(a, b, c, d) { const o1 = orient(a, b, c); const o2 = orient(a, b, d); const o3 = orient(c, d, a); const o4 = orient(c, d, b); const s1 = o1.n === 0n ? 0 : (o1.n > 0n ? 1 : -1); const s2 = o2.n === 0n ? 0 : (o2.n > 0n ? 1 : -1); const s3 = o3.n === 0n ? 0 : (o3.n > 0n ? 1 : -1); const s4 = o4.n === 0n ? 0 : (o4.n > 0n ? 1 : -1); if (s1 === 0 && onSegment(a, b, c)) return true; if (s2 === 0 && onSegment(a, b, d)) return true; if (s3 === 0 && onSegment(c, d, a)) return true; if (s4 === 0 && onSegment(c, d, b)) return true; return (s1 !== s2) && (s3 !== s4); } function pointInConvexPoly(pt, poly) { let sign = 0; for (let i = 0; i < poly.length; i++) { const a = poly[i]; const b = poly[(i + 1) % poly.length]; const o = orient(a, b, pt); const s = o.n === 0n ? 0 : (o.n > 0n ? 1 : -1); if (s === 0) continue; if (sign === 0) sign = s; else if (sign !== s) return false; } return true; } function intersectPolyExact(T, R) { for (let i = 0; i < T.length; i++) { const a = T[i], b = T[(i + 1) % T.length]; for (let j = 0; j < R.length; j++) { const c = R[j], d = R[(j + 1) % R.length]; if (segIntersects(a, b, c, d)) return true; } } if (pointInConvexPoly(T[0], R)) return true; if (pointInConvexPoly(R[0], T)) return true; return false; } function intersectPoly(T, R) { // numeric bbox reject first (fast) if (!bboxOverlapN(bboxPolyN(T), bboxPolyN(R))) return false; // then exact return intersectPolyExact(T, R); } // ------------------------- // View helpers // ------------------------- function rectFromView(view) { const { xminR, yminR, ddR } = view; const xmaxR = rAdd(xminR, ddR); const ymaxR = rAdd(yminR, ddR); return [ [xminR, yminR], [xmaxR, yminR], [xmaxR, ymaxR], [xminR, ymaxR] ]; } function computeB(c0, ddR) { const ddN = rToNumber(ddR); return (800 * (Number((-c0[2])) ** 0.3)) / ddN; } // ------------------------- // Gasket generator with pixel-radius cutoff (key for speed) // Exact circles; cutoff only decides "too small to see in this view". // ------------------------- function gasket(c0, c1, c2, c3, rect, Bnum, opts = {}) { const { maxCircles = 40000, minPixRadius = 0.35, // smaller => more circles; keep low for deep zoom proj = null // expects { rPixN(kBig) } } = opts; const out = [c0, c1, c2, c3]; let count = 4; function* circlesGen(a, b, c, d) { if (count >= maxCircles) return; const cc = circle(a, b, c, d); const k = cc[2]; // curvature cutoff if (!(Number(k < 0n ? -k : k) < Bnum)) return; // pixel cutoff (view dependent) — huge speed win if (proj) { const rp = proj.rPixN(k); if (rp < minPixRadius) return; } // pruning region test let yield_more = false; const k0 = a[2], k1 = b[2], k2 = c[2]; if (k0 > 0n && k1 > 0n && k2 > 0n) { const T = inner_triangle(a, b, c); if (intersectPoly(T, rect)) yield_more = true; } else { if (k0 < 0n) { const A = [b[0], b[1]]; const B = [c[0], c[1]]; const P = [cc[0], cc[1]]; if (some_sameSide(A, B, P, rect)) yield_more = true; } else if (k1 < 0n) { const A = [a[0], a[1]]; const B = [c[0], c[1]]; const P = [cc[0], cc[1]]; if (some_sameSide(A, B, P, rect)) yield_more = true; } else if (k2 < 0n) { const A = [a[0], a[1]]; const B = [b[0], b[1]]; const P = [cc[0], cc[1]]; if (some_sameSide(A, B, P, rect)) yield_more = true; } } if (!yield_more) return; yield cc; yield* circlesGen(a, b, cc, c); yield* circlesGen(b, cc, c, a); yield* circlesGen(c, cc, a, b); } const gens = [ circlesGen(c1, c2, c3, c0), circlesGen(c2, c3, c0, c1), circlesGen(c3, c0, c1, c2), circlesGen(c0, c1, c2, c3) ]; for (const g of gens) { for (const cc of g) { out.push(cc); count++; if (count >= maxCircles) break; } if (count >= maxCircles) break; } return out; } // ------------------------- // Examples // ------------------------- const examples = [ [[0, 0, -1], [-1 / 2, 0, 2], [1 / 2, 0, 2], [0, -2 / 3, 3]], [[0, 0, -2], [-1 / 6, 0, 3], [1 / 3, 0, 6], [3 / 14, -2 / 7, 7]], [[0, 0, -3], [-2 / 15, 0, 5], [1 / 6, -1 / 8, 8], [1 / 6, 1 / 8, 8]], [[0, 0, -6], [-5 / 66, 0, 11], [8 / 105, -2 / 35, 14], [3 / 50, 4 / 50, 15]], [ [0, 0, -13], [-10 / 299, 0, 23], [14 / 325, -2 / 300, 30], [22 / 1235, 9 / 190, 38] ], [ [0, 0, -76], [-1 / 5852, 0, 77], [1 / 77, 0, 5852], [2000 / 154053, -2 / 5853, 5853] ] ].map((quad) => quad.map(toBigCircle)); const exampleNames = [ "(-1; 2,2,3)", "(-2; 3,6,7)", "(-3; 5,8,8)", "(-6; 11,14,15)", "(-13; 23,30,38)", "(-76; 77,5852,5853)" ]; return { // rationals Rat, rNeg, rAdd, rSub, rMul, rDiv, rAbs, rInv, rToNumber, rFromNumber, ratApprox, rCmp, // circles + geometry toBigCircle, circle, tangency, inner_triangle, // view + generation rectFromView, computeB, gasket, // examples examples, exampleNames }; } 33:Taf3,## Mathematical notes An [Apollonian packing](https://en.wikipedia.org/wiki/Apollonian_gasket) is generated from a quadruple of mutually tangent circles by recursively filling in more circles — each tangent to three of the circles in the configuration. In the special cases shown here the reciprocal of the radius of each circle is an integer — a configuration sometimes called a *bowl of integers*. This notebook displays several such packings and allows you to zoom in for a closer look. Of course, zooming in will eventually produce numerical errors. ### Curvature and Descarte's circle theorem The reciprocal of the radius of a circle is known more succinctly as its *curvature*, which are shown by the labels in the diagram. If the curvature of four mutually tangent circles are ${tex`\kappa_1`}, ${tex`\kappa_2`}, ${tex`\kappa_3`}, and ${tex`\kappa_4`}, then they must satisfy [Descartes' circle theorem](https://en.wikipedia.org/wiki/Descartes%27_theorem), namely: ${tex.block`(\kappa_1+\kappa_2+\kappa_3+\kappa_4)^2 = 2(\kappa_1^2+\kappa_2^2+\kappa_3^2+\kappa_4^2).`} This is true even for the exterior circle, assuming we take its curvature to be negative. *With this formula, one can find two solutions of the fourth curvature given the first three; or for our purposes, given one of the two possible values of ${tex`\kappa_4`}, we can find the other, ${tex`\kappa_4'`}, by a *linear* relation* ${tex.block`\kappa_4+\kappa_4'=2\kappa_1+2\kappa_2+2\kappa_3`} *so that, if we start out with 4 integer curvatures, all the resulting curvatures will be integers! The map ${tex`(\kappa_1, \kappa_2, \kappa_3, \kappa_4)\mapsto (\kappa_1, \kappa_2, \kappa_3, \kappa_4')`} is an involution on ${tex`\mathbb Z^4`}:* ${tex.block` S_4 = \begin{pmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 & 1 & 0 \\ 2 & 2 & 2 & -1\end{pmatrix} \in GL(4,\mathbb Z)`} *and similarly for ${tex`S_1, S_2, S_3`} that are "reflections" through three other mirrors. The group they generate is called the Apollonian group.* *Remarkably still, the centers of the four circles, when expressed as complex numbers and weighted by the respective curvatures, satisfy the *same* quadratic relation, a result apparently only discovered in the late 1990s when trying to generate this kind of picture on a computer; see https://arxiv.org/abs/math/0101066. One of the consequences is that the centers are all *rational* (if they start out that way).* ### *Number theory* *The study of Apollonian packings has become part of number theory. Among the many deep results, it is known (Sarnak, 2007) that there are infinitely many primes (shown in red) that show up in any integer Apollonian packing; in fact, there are infinitely many primes that surround any given circle (hence in particular, infinitely many kissing primes)!*34:[] 2f:{"title":"Integer Apollonian Packings","cells":[{"id":0,"value":"# Integer Apollonian Packings\n\n*LLM takes in the work of the original author, [Mark McClure](https://observablehq.com/@mcmcclur), and realizes my suggestions:*\n\n- *whole numbers are now in `BigInt`;*\n- TODO: make rationals as pairs of `BigInt`;\n- TODO: make it go deeper;\n- *show prime numbers.*","pinned":false,"hidden":false,"mode":"md"},{"id":190,"value":"$31","pinned":false,"hidden":false,"mode":"ojs"},{"id":35,"value":"$32","pinned":false,"hidden":false,"mode":"ojs"},{"id":1089,"value":"$33","pinned":false,"hidden":false,"mode":"md"},{"id":791,"value":"## Code\n\nrefer to the original.","pinned":false,"hidden":false,"mode":"md"},{"id":1077,"value":"### Initial Appollonian packings\n\nThere are infinitely many integer Apollonian packings but they are not necessarily easy to find. Here's a list of the six that are displayed above.","pinned":false,"hidden":false,"mode":"md"},{"id":826,"value":"aps = [\n [[0, 0, -1n], [-1/2, 0, 2n], [1/2, 0, 2n], [0, -2/3, 3n]],\n [[0, 0, -2n], [-1/6, 0, 3n], [1/3, 0, 6n], [3/14, -2/7, 7n]],\n [[0, 0, -3n], [-2/15, 0, 5n], [1/6, -1/8, 8n], [1/6, 1/8, 8n]],\n [[0, 0, -6n], [-5/66, 0, 11n], [8/105, -2/35, 14n], [3/50, 4/50, 15n]],\n [\n [0, 0, -13n],\n [-10/299, 0, 23n],\n [14/325, -2/300, 30n],\n [22/1235, 9/190, 38n]\n ],\n [\n [0, 0, -76n],\n [-1/5852, 0, 77n],\n [1/77, 0, 5852n],\n [2000/154053, -2/5853, 5853n]\n ]\n].map(x=>{ })","pinned":true,"hidden":false,"mode":"ojs"},{"id":1080,"value":"### Imports","pinned":false,"hidden":false,"mode":"md"},{"id":372,"value":"Flatten = require('@flatten-js/core')","pinned":true,"hidden":false,"mode":"ojs"},{"id":4,"value":"d3 = require('d3@5')","pinned":true,"hidden":false,"mode":"ojs"},{"id":1144,"value":"html`\n`","pinned":true,"hidden":false,"mode":"ojs"}],"files":"$W34","stdlib":"1","titling":"auto"}