25:["$","title",null,{"children":"A formula for Hilbert's curve | Mark McClure | Observable"}] 49:T42e,## Comments If you've seen approximations to Hilbert's curve before, they probably don't look like the approximations above. The ${tex`k^{\text{th}}`} approximation in the interactive image above is obtained by computing the points ${tex.block`h(p/2^k) \text{ for } p\in\{0,1,2,\ldots,2^k\}`} and then connecting the dots. That was *not* Hilbert's approach and that distinction is rather the point of this notebook. ### Hilbert's construction Hilbert constructs his curve in a particular way to help him *prove* that the curve has the specific properties he wants: namely, he wants the curve to be > *continuous* and *onto*. There are zillions of computer generated images of Hilbert's construction, including a few right here on Observable. You could hardly find a more authoritative reference, though, than Hilbert's own rendition, which appears in his [paper of 1891](https://www.digizeitschriften.de/id/235181684_0038|log40?tify=%7B%22pages%22%3A%5B476%5D%2C%22pan%22%3A%7B%22x%22%3A0.474%2C%22y%22%3A0.71%7D%2C%22view%22%3A%22info%22%2C%22zoom%22%3A0.485%7D):4a:T5cf,### A formula A reasonable question, as asked in [this Math.StackExchange](https://math.stackexchange.com/questions/921985/what-is-the-hilbert-curves-equation/) post, that's not addressed in Hilbert's paper is: > Is there a formula for Hilbert's curve? It turns out that there is such a formula discovered by Hans Sagan in the early 1990s. It was published in his fun little book called [Space Filling Curves](https://link.springer.com/book/10.1007/978-1-4612-0871-6), where it appears as formula 2.4.3 on page 18. Sagan's formula, is expressed in terms of an infinite sum involving the base four expansion of the input ${tex`t\in[0,1)`}. That is, we write ${tex.block`t=0_{\dot 4}q_1q_2q_3\ldots,`} where each ${tex`q_i`} is ${tex`0`}, ${tex`1`}, ${tex`2`}, or ${tex`3`}. Then, ${tex.block`h(t) = \sum_{j=1}^{\infty} \frac{(-1)^{e_{0j}}}{2^j}\text{sign}(q_j) \left( \begin{array}{c} (1-d_j)q_j - 1 \\ 1-d_jq_j \end{array} \right),`}` where ${tex`e_{kj}`} denotes the number of ${tex`k`}s preceding ${tex`q_j`} and ${tex`d_j=e_{0j}+e_{3j} \mod 2`}. The notation ${tex`\text{sign}(x)`} denotes [the signum function](https://en.wikipedia.org/wiki/Sign_function) that's ${tex`1`} if ${tex`x`} is positive, ${tex`0`} if ${tex`x`} is zero, or ${tex`-1`} if ${tex`x`} is negative (though ${tex`q_j`} never is negative). I suppose this isn't the simplest thing to grasp but it's not too bad to program. For me, programming a mathematical concept is often a step towards understanding.4b:[{"id":"b1ede53a825fd1befb33d31998ed453519ca116118cea1fff50347bdc6f96ae56d6cc3b089c2d77a5668f919a95406eb84a1f76598c552741863276ead0c16fd","href":"https://static.observableusercontent.com/files/b1ede53a825fd1befb33d31998ed453519ca116118cea1fff50347bdc6f96ae56d6cc3b089c2d77a5668f919a95406eb84a1f76598c552741863276ead0c16fd","name":"Screen Shot 2023-08-10 at 5.37.34 PM.png","type":"image/png","lastModified":1691703726735.249,"size":543379}] 29:{"id":"14837cb8287ef850","slug":null,"title":"A formula for Hilbert's curve","thumbnail":"86ccaa2592953d63725963b92d4e52fb7004dd4912d9c4fd24308af15b3c5c90","access_level":"public","create_time":"$D2023-08-10T19:57:44.762Z","update_time":"$D2026-05-09T16:57:40.023Z","view_time":null,"first_public_time":"$D2023-08-11T15:03:57.323Z","has_unpublished_edits":null,"fork":false,"stars":3,"starred":false,"owner":{"type":"team","id":"d771c8ae812d5130","login":"mcmcclur","name":"Mark McClure","avatar_url":"https://avatars-next.observableusercontent.com/e615e819ce83d43faef607e48c53c5ff63943aa0ba9f5848499bd02fef7396f4","tier":"starter_2024","bio":"I'm a professor of mathematics at the University of North Carolina Asheville. I've also done a fair amount of consulting work over the years focusing on scientific and data visualization.\n\nI'm happy to be part of the first Observable Ambassador's cohort.","followed_by_you":false,"follows_you":false},"version":408,"publish_version":null,"role":null,"authored":false,"authors":[{"type":"individual","id":"c3301e515dd08a41","login":"mcmcclur","name":"Mark McClure","avatar_url":"https://avatars-next.observableusercontent.com/e615e819ce83d43faef607e48c53c5ff63943aa0ba9f5848499bd02fef7396f4","tier":null,"bio":"I'm a professor of mathematics at the University of North Carolina Asheville. I've also done a fair amount of consulting work over the years focusing on scientific and data visualization.\n\nI'm happy to be part of the first Observable Ambassador's cohort.","followed_by_you":false,"follows_you":false}],"members":[{"type":"individual","id":"c3301e515dd08a41","login":"mcmcclur","name":"Mark McClure","avatar_url":"https://avatars-next.observableusercontent.com/e615e819ce83d43faef607e48c53c5ff63943aa0ba9f5848499bd02fef7396f4","tier":null,"role":"editor"}],"comments":[],"custom_thumbnail":true,"trashed":false,"lapsed":null,"subscribed":null,"forks":1,"collection":null,"type":"notebook","body":{"titling":"auto","title":"A formula for Hilbert's curve","cells":[{"id":0,"value":"# A formula for Hilbert's curve\n\nHilbert's space filling curve is a continuous function ${tex`h`} that maps the unit interval onto the unit square. That is, it accepts a number ${tex`t`} such that ${tex`0\\leq t \\leq 1`} and it returns a point ${tex`(x,y)`} where ${tex`x`} and ${tex`y`} are also between zero and one.\n\nYou can use the slider below to generate an ${tex`n^{\\text{th}}`} level approximation of Hilbert's curve.","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"","database":"$undefined"},{"id":87,"value":"viewof n = Inputs.range([0, 14], { value: 6, step: 1 })","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":106,"value":"pic = curve(n)","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":132,"value":"$49","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"","database":"$undefined"},{"id":134,"value":"FileAttachment(\"Screen Shot 2023-08-10 at 5.37.34 PM.png\").image()","pinned":false,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":139,"value":"The numbering scheme allows him to show that the nesting structure of the squares corresponds to exactly the same nesting structure in the interval; this allows him to derive its properties.","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"","database":"$undefined"},{"id":142,"value":"$4a","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"","database":"$undefined"},{"id":285,"value":"## Implementation\n\nHere's where we actually define the Javascript function `h(t)`, some auxiliary functions, and a function to produce the image.\n\nNote that, if ${tex`t`} has the form ${tex`p/2^n`}, then the sum in ${tex`h(t)`} will terminate and the computation is exact. For any other rational number, the base four expansion is periodic and the value can still be computed exactly using a geometric series. In [this Math.StackExchange post](https://math.stackexchange.com/questions/921985/what-is-the-hilbert-curves-equation/1739118#1739118), for example, it's shown that\n${tex.block`h\\left(\\frac{4}{7}\\right) = \\left(\n\\begin{array}{c}\n 38/63 \\\\\n 52/63 \\\\\n\\end{array}\n\\right).`}\n\nFor irrational numbers, the sum can be truncated to obtain a good decimal approximation.","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"","database":"$undefined"},{"id":334,"value":"## Code","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"","database":"$undefined"},{"id":12,"value":"// Accepts a real number between 0 and 1 and computes the image of\n// that number in the unit square.\n// Should be exact if the input is a rational number whose \n// denominator is a power of 2.\n// Yields an approximation for all other numbers.\nfunction h(t) {\n if (t == 0) {\n return [0, 0];\n } else if (t == 1) {\n return [1, 0];\n } else if (typeof t == \"number\" && 0 < t && t < 1) {\n let qs = t\n .toString(4)\n .split(\".\")[1]\n .split(\"\")\n .map((q) => parseInt(q));\n return h0(qs);\n }\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":312,"value":"// Example\nh(3 / 8)","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":38,"value":"// Accepts an array of numbers, denoted qs, that are assumed to represent the\n// base four expansion of a real number in the unit interval.\n// Thus, each entry in qs should be 0, 1, 2, or 3.\n// Returns the image of the number\nfunction h0(qs) {\n let cs = qs.map((q, j) => ((-1) ** e(0, j, qs) / 2 ** j) * Math.sign(q));\n let x = d3.sum(qs.map((q, j) => cs[j] * ((1 - d(j, qs)) * q - 1))) / 2;\n let y = d3.sum(qs.map((q, j) => cs[j] * (1 - d(j, qs) * q))) / 2;\n return [x, y];\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":317,"value":"// Example\nh0([1, 0, 3])","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":301,"value":"Sagan's ${tex`d_j`}:","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"","database":"$undefined"},{"id":10,"value":"d = (j, qs) => (e(0, j, qs) + e(3, j, qs)) % 2","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":304,"value":"Sagan's ${tex`e_{kj}`}:","pinned":false,"hidden":false,"mode":"md","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"","database":"$undefined"},{"id":4,"value":"function e(k, j, qs) {\n let count = 0;\n for (let i = 0; i < j; i++) {\n if (qs[i] === k) {\n count++;\n }\n }\n return count;\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":283,"value":"// Generate the nth level approximation to the curve.\nfunction curve(n, size = 800) {\n let pts = d3\n .range(0, 1, 1 / 2 ** n)\n .concat([1])\n .map(function (t) {\n let [x, y] = h(t);\n return { t, x, y };\n });\n\n return Plot.plot({\n width: size,\n height: size,\n inset: 10,\n x: { domain: [0, 1], ticks: 8 },\n y: { domain: [0, 1], ticks: 8 },\n marks: [\n Plot.line(pts, { x: \"x\", y: \"y\", strokeWidth: 3 }),\n Plot.dot(pts, {\n fill: \"black\",\n r: 5,\n x: \"x\",\n y: \"y\",\n channels: { t: \"t\" },\n tip: true\n }),\n Plot.frame()\n ]\n });\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":102,"value":"// // This code is used to produce the image in the StackExchange post\n// {\n// let container = d3.create(\"div\");\n// d3.range(12).forEach(function (n) {\n// let div = container\n// .append(\"div\")\n// .style(\"width\", `${width / 4}px`)\n// .style(\"display\", \"inline-block\");\n// div.append(() => curve(n, width / 3));\n// container.append(() => div.node());\n// });\n\n// return container.node();\n// }\n","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":365,"value":"curve(5, width / 2)","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":395,"value":"hilbertExact(4, 7)","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":372,"value":"function hilbertExact(a, b) {\n // assumes 0 <= a < b; handle t = 1 separately as 0.3333...\n let r = a;\n let e0 = 0;\n let e3 = 0;\n let j = 1;\n\n const seen = new Map();\n const terms = [];\n\n while (true) {\n if (r === 0) {\n // terminating base-4 expansion; remaining q_j are 0\n // and contribute zero because sign(0) = 0\n return finiteSum(terms);\n }\n\n const state = `${r},${e0},${e3}`;\n\n if (seen.has(state)) {\n const m = seen.get(state);\n const L = j - m;\n return prefixPlusGeometricTail(terms, m, L);\n }\n\n seen.set(state, j);\n\n const x = 4 * r;\n const q = Math.floor(x / b);\n r = x % b;\n\n const d = (e0 + e3) % 2;\n const sgn = e0 === 0 ? 1 : -1;\n const sigq = q === 0 ? 0 : 1;\n\n const vx = sgn * sigq * ((1 - d) * q - 1);\n const vy = sgn * sigq * (1 - d * q);\n\n terms.push([vx, vy]); // this is v_j\n\n if (q === 0) e0 = 1 - e0;\n if (q === 3) e3 = 1 - e3;\n\n j += 1;\n }\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":380,"value":"function gcd(a, b) {\n a = a < 0n ? -a : a;\n b = b < 0n ? -b : b;\n\n while (b !== 0n) {\n [a, b] = [b, a % b];\n }\n\n return a;\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":382,"value":"\nfunction reduce(num, den) {\n if (den < 0n) {\n num = -num;\n den = -den;\n }\n\n const g = gcd(num, den);\n return [num / g, den / g];\n}\n\n","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":385,"value":"function addFrac([a, b], [c, d]) {\n return reduce(a * d + b * c, b * d);\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":387,"value":"function scaleFrac([a, b], k) {\n return reduce(a * k, b);\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":389,"value":"function subFrac([a, b], [c, d]) {\n return reduce(a * d - b * c, b * d);\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":391,"value":"function finiteSum(terms) {\n let x = [0n, 1n];\n let y = [0n, 1n];\n\n for (let j = 0; j < terms.length; j++) {\n const denom = 1n << BigInt(j + 1);\n\n const [vx, vy] = terms[j];\n\n x = addFrac(x, [BigInt(vx), denom]);\n y = addFrac(y, [BigInt(vy), denom]);\n }\n\n return { x, y };\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":393,"value":"function prefixPlusGeometricTail(terms, m, L) {\n // prefix\n let px = [0n, 1n];\n let py = [0n, 1n];\n\n for (let j = 1; j < m; j++) {\n const denom = 1n << BigInt(j);\n\n const [vx, vy] = terms[j - 1];\n\n px = addFrac(px, [BigInt(vx), denom]);\n py = addFrac(py, [BigInt(vy), denom]);\n }\n\n // one repeating block\n let bx = [0n, 1n];\n let by = [0n, 1n];\n\n for (let j = m; j < m + L; j++) {\n const denom = 1n << BigInt(j);\n\n const [vx, vy] = terms[j - 1];\n\n bx = addFrac(bx, [BigInt(vx), denom]);\n by = addFrac(by, [BigInt(vy), denom]);\n }\n\n // divide by (1 - 2^{-L})\n const geomDen = (1n << BigInt(L)) - 1n;\n\n bx = scaleFrac(bx, 1n << BigInt(L));\n by = scaleFrac(by, 1n << BigInt(L));\n\n bx = [bx[0], bx[1] * geomDen];\n by = [by[0], by[1] * geomDen];\n\n bx = reduce(...bx);\n by = reduce(...by);\n\n return {\n x: addFrac(px, bx),\n y: addFrac(py, by)\n };\n}","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"},{"id":398,"value":"hilbertExact(1, 12)","pinned":true,"hidden":false,"mode":"ojs","author":"$undefined","summary":"$undefined","tool":"$undefined","providerMetadata":"$undefined","attachments":"$undefined","output":"$undefined","database":"$undefined"}],"files":"$W4b","stdlib":"1"}} 24:null