Plot a hyperbolic plane tessellation

Question

Make a plot (Poincare disk) of a tessellation on a hyperbolic plane, such as:

The program takes four inputs:

1) How many edges/polygon (three in this example).

2) How many intersect at each vertex (seven in this example).

3) How many steps away from the center vertex to render (5 in this example, if you look closely). This means that a vertex is included iff it can be reached in 5 or less steps form the center. Edges are rendered iff both their vertexes are included.

4) The resolution of the image (a single number of pixels, the image is square).

The output must be an image. Edges must be rendered as circle-arcs, not lines (the Poincaré disk projection turns lines into circles). Points do not need to be rendered. When the user puts in something that is not hyperbolic (i.e. 5 triangles meeting at each vertex), the program does not have to work properly. This is code-golf, so the shortest answer wins.

Answer 1

Mathematica, 2535 bytes

Taken from here (hence why it's community wiki). Not really that golfed. View the provided link for the author's explanation of his code.

Also, I'm no Mathematica expert, but I bet Martin could do wonders on the code length. I don't even understand the math behind it.

I left it readable, but if the question doesn't get closed, I'll golf it past readability and move the 2 other parameters inside the caller function.

Currently invalid, feel free to help improve it:

• I think this uses lines rather than arcs.

• Centered on a face, rather than a vertex.

HyperbolicLine[{{Px_, Py_}, {Qx_, Qy_}}] := 
 If[N[Chop[Px Qy - Py Qx]] =!= 0., 
  Circle[OrthoCentre[{{Px, Py}, {Qx, Qy}}], 
   OrthoRadius[{{Px, Py}, {Qx, Qy}}], 
   OrthoAngles[{{Px, Py}, {Qx, Qy}}]], Line[{{Px, Py}, {Qx, Qy}}]]

OrthoCentre[{{Px_, Py_}, {Qx_, Qy_}}] := 
 With[{d = 2 Px Qy - 2 Py Qx, p = 1 + Px^2, q = 1 + Qx^2 + Qy^2}, 
  If[N[d] =!= 0., {p Qy + Py^2 Qy - Py q, -p Qx - Py^2 Qx + Px q}/d, 
   ComplexInfinity]]

OrthoRadius[{{Px_, Py_}, {Qx_, Qy_}}] := 
 If[N[Chop[Px Qy - Py Qx]] =!= 0., 
  Sqrt[Total[OrthoCentre[{{Px, Py}, {Qx, Qy}}]^2] - 1], Infinity]

OrthoAngles[{{Px_, Py_}, {Qx_, Qy_}}] := 
 Block[{a, b, c = OrthoCentre[{{Px, Py}, {Qx, Qy}}]}, 
  If[(a = N[Apply[ArcTan, {Px, Py} - c]]) < 0., a = a + 2 \[Pi]];
  If[(b = N[Apply[ArcTan, {Qx, Qy} - c]]) < 0., 
   b = b + 2 \[Pi]]; {a, b} = Sort[{a, b}];
  If[b - a > \[Pi], {b, a + 2 \[Pi]}, {a, b}]]

Inversion[Circle[{Cx_, Cy_}, r_], {Px_, Py_}] := {Cx, Cy} + 
  r^2 {Px - Cx, Py - Cy}/((Cx - Px)^2 + (Cy - Py)^2)
Inversion[Circle[{Cx_, Cy_}, r_, {a_, b_}], {Px_, Py_}] := {Cx, Cy} + 
  r^2 {Px - Cx, Py - Cy}/((Cx - Px)^2 + (Cy - Py)^2)

Inversion[Circle[{Cx_, Cy_}, r_, {a_, b_}], p_Line] := 
 Map[Inversion[Circle[{Cx, Cy}, r], #] &, p, {2}]

Inversion[Circle[{Cx_, Cy_}, r_, {a_, b_}], p_Polygon] := 
 Map[Inversion[Circle[{Cx, Cy}, r], #] &, p, {2}]

Inversion[Line[{{Px_, Py_}, {Qx_, Qy_}}], {Ux_, Uy_}] := 
 With[{u = Px - Qx, 
   v = Qy - Py}, {-Ux (v^2 - u^2) - 2 u v Uy, 
    Uy (v^2 - u^2) - 2 u v Ux}/(u^2 + v^2)]
Inversion[Line[{{Px_, Py_}, {Qx_, Qy_}}], p_Polygon] := 
 Map[Inversion[Line[{{Px, Py}, {Qx, Qy}}], #] &, p, {2}]

Inversion[Circle[{Cx_, Cy_}, r_], c_List] := 
 Map[Inversion[Circle[{Cx, Cy}, r], #] &, c]


PolygonInvert[p_Polygon] := 
 Map[Inversion[HyperbolicLine[#], p] &, 
  Partition[Join[p[[1]], {p[[1, 1]]}], 2, 1]]
PolygonInvert[p_List] := Flatten[Map[PolygonInvert[#] &, p]]

LineRule = Polygon[x_] :> Line[Join[x, {x[[1]]}]];
HyperbolicLineRule = 
  Polygon[x_] :> 
   Map[HyperbolicLine, Partition[Join[x, {x[[1]]}], 2, 1]];

CentralPolygon[p_Integer, q_Integer, \[Phi]_: 0] := 
 With[{r = (Cot[\[Pi]/p] Cot[\[Pi]/q] - 1)/
     Sqrt[Cot[\[Pi]/p]^2 Cot[\[Pi]/q]^2 - 1], \[Theta] = \[Pi] Range[
       1, 2 p - 1, 2]/p}, 
  r Map[{{Cos[\[Phi]], -Sin[\[Phi]]}, {Sin[\[Phi]], Cos[\[Phi]]}}.# &,
     Transpose[{Cos[\[Theta]], Sin[\[Theta]]}]]]

PolygonUnion[p_Polygon, tol_: 10.^-10] := p
PolygonUnion[p_List, tol_: 10.^-10] := 
 With[{q = p /. Polygon[x_] :> N[Polygon[Round[x, 10.^-10]]]}, 
  DeleteDuplicates[q]]
HyperbolicTessellation[p_Integer, q_Integer, \[Phi]_, k_Integer, 
  t_: 10.^-10] := 
 Map[PolygonUnion[#, t] &, 
   NestList[PolygonInvert, Polygon[CentralPolygon[p, q, \[Phi]]], 
     k][[{-2, -1}]]] /; k > 0

HyperbolicTessellation[p_Integer, q_Integer, \[Phi]_, k_Integer, 
  t_: 10.^-10] := Polygon[CentralPolygon[p, q, \[Phi]]] /; k == 0
HyperbolicTessellationGraphics[p_Integer, q_Integer, \[Phi]_, 
  k_Integer, rule_RuleDelayed, opts___] := 
 Graphics[{Circle[{0, 0}, 1], 
   HyperbolicTessellation[p, q, \[Phi], k, 10.^-10] /. rule}, opts]

Called like:

HyperbolicTessellationGraphics[7, 3, 0., 5, HyperbolicLineRule, ImageSize -> 100, PlotLabel -> "{7,3}"]
HyperbolicTessellationGraphics[3, 7, 0., 7, HyperbolicLineRule, ImageSize -> 300, PlotLabel -> "{7,7}"]

Answer 2

JavaScript/HTML/SVG, 5086 bytes

Not golfed, but the logic/math is short and hopefully relatively clear; also does the planar and spherical tessellations ((p-2)*(q-2) <= 4) without too much extra logic (it was easier to implement those than to leave gaps in the table!).

The calling code in the snippet below calls the function several times: first to display the 512x512 {3,7} picture (not shown here; it looks just like the above), and then to display a 128x128 {p,q} (meaning q p-gons surround each vertex) for each p,q from 3 to 9:

<!DOCTYPE html>
<html>
 <!--
 https://codegolf.stackexchange.com/questions/55848/plot-a-hyperbolic-plane-tessellation
 -->
 <style> svg { border: 1px solid black; } </style>
 <body>
  <script>
   "use strict";
   console.log("    in script");
   // Do an immediately invoked async function,
   // so we can sleep in it (to flush drawing).
   (async () => {
 console.log("        in script async function");
 const t0 = performance.now();

 // ==========================
 // IMPLEMENTATION BEGINS HERE

 // some complex math
 const neg = ([x,y]) => [-x,-y];
 const conj = ([x,y]) => [x,-y];
 const plus = (a,b) => [a[0]+b[0],a[1]+b[1]];
 const minus = (a,b) => [a[0]-b[0],a[1]-b[1]];
 const cross = (a,b) => a[0]*b[1]-a[1]*b[0];  // area of spanned parallelogram
 const times = (a,b) => [a[0]*b[0]-a[1]*b[1], a[0]*b[1]+a[1]*b[0]];
 const timesreal = ([x,y],r) => [x*r, y*r];
 const abs2 = ([x,y]) => x**2+y**2;
 const abs = z => Math.sqrt(abs2(z));
 const dist = (a,b) => abs(minus(a,b));
 const inverse = z => timesreal(conj(z), 1/abs2(z));
 const dividedby = (a,b) => times(a,inverse(b));

 // Transform z by the translation that takes the origin to t.
 // Works with any curvature; in particular:
 //   -1 = hyperbolic (poincare disk)
 //    0 = planar,
 //    1 = spherical (stereographic projection)
 const Plus = (z,t, curvature) => {
   // (z + t) / (1 - curvature*conj(t)*z)
   return dividedby(plus(z,t),
                    minus([1,0],
                          timesreal(times(conj(t), z), curvature)));
 };
 // Transform z by the translation that takes t to the origin.
 const Minus = (z,t, curvature) => Plus(z, neg(t), curvature);

 // Simple monotonic function of actual distance. Specifically:
 // - hyperbolic distance is 2*atanh(chord distance)
 // - spherical distance is 2*atan(chord distance)
 // - planar distance is chord distance
 // (I might be off by a factor of 2 in some of these claims, doesn't matter)
 const ChordDistance = (a,b, curvature) => abs(Minus(a,b, curvature));

 // signed radius of the circle passing through a,b,c
 const circumradius = (a,b,c) => dist(a,b)*dist(b,c)*dist(c,a) /
                                 (2 * cross(minus(b,a), minus(c,a)));

 // signed radius of circle needed to render segment a,b as circular arc
 // in the picture; may be infinite
 const ArcRadius = (a,b, curvature) => {
   if (curvature === 0) return Infinity;  // planar picture: all straight lines
   if (abs(a) > 1e3 || abs(b) > 1e3) return Infinity;  // one of them is infinite
   if (Math.abs(cross(a,b)) < 1e-3) return Infinity;  // collinear with origin
   // Let c be any third point on the arc
   const c = Plus(timesreal(Minus(b,a, curvature), .5), a, curvature);
   return circumradius(a,b,c);
 };


 const doit = (svgsize,p,q,steps,maxverts) => {
  const curvature = Math.sign(4-(p-2)*(q-2));  // -1:hyperbolic, 0:planar, 1:spherical
  const Dist = (a,b) => ChordDistance(a,b, curvature);

  const euclidean_first_edge_length = curvature===0 ? 1 :
    Math.sqrt(Math.abs((Math.sin(Math.PI/q)/Math.cos(Math.PI/p))**2 - 1))

  // a,b are generators of the symmetry group (not including reflections):
  //  - a is rotation by 2pi/q about origin
  //  - b is rotation by pi about first edge midpoint
  const primitive_qth_root_of_unity = [Math.cos(2*Math.PI/q),
                                       Math.sin(2*Math.PI/q)];
  const a = z => times(z, primitive_qth_root_of_unity);
  const b = ([x,y]) => Plus([-x, -y], [0,-euclidean_first_edge_length], curvature);

  // generate verts
  const vs = [[0,0]];
  for (let i = 0; i < steps; ++i) {
    for (const g of [b,a]) {
      for (let j = 0; j < vs.length; ++j) {  // while vs is growing
        if (vs.length >= maxverts) break;
        const v = g(vs[j]);
        if (vs.findIndex(w=>Dist(v,w)<=1e-6)<0) vs.push(v);
      }
    }
  }

  // generate edges (pairs of points, not pairs of indices)
  const es = [];
  if (vs.length >= 2) {
    const d = Dist(vs[0], vs[1]);
    for (let i = 0; i < vs.length; ++i) {
      for (let j = 0; j < i; ++j) {
        if (Math.abs(Dist(vs[i],vs[j]) - d) <= 1e-6) es.push([vs[i], vs[j]]);
      }
    }
  }

  const scale =
      curvature>0 ? .35 :  // spherical: get all verts in the picture
      curvature===0 ? 1/5.5 :  // planar: get a good sampling in the picture
      1;  // hyperbolic: get the poincare disk in the picture

  let svg_html = '<svg width="'+svgsize+'" height="'+svgsize+'"><path d="';
  for (let [a,b] of es) {
    const r = ArcRadius(a,b,curvature)*scale/2*svgsize;

    // if a or b is at infinity, then draw straight outward
    // from the other point
    if (abs(a) > 1e2) a = timesreal(b, 1e2);
    if (abs(b) > 1e2) b = timesreal(a, 1e2);

    const x0 = (a[0]*scale+1)/2*svgsize;
    const y0 = (a[1]*scale+1)/2*svgsize;
    const x1 = (b[0]*scale+1)/2*svgsize;
    const y1 = (b[1]*scale+1)/2*svgsize;

    if (!Number.isFinite(r)) {
      svg_html += 'M'+x0+','+y0+'L'+x1+','+y1;  // straight line
    } else {
      // circular arc of signed radius r
      svg_html += 'M'+x0+','+y0+'A'+r+','+r+',0,0,'+(r<0?1:0)+','+x1+','+y1;
    }
  }
  svg_html += '" stroke="black" fill="none"/></svg>';
  return svg_html;
 }; // doit

 // IMPLEMENTATION ENDS HERE
 // ========================


 document.body.innerHTML += doit(512, 3, 7, 5, 1000)
 document.body.innerHTML += "<br>";

 const sleep = ms => new Promise(resolve => setTimeout(resolve, ms));

 const svgsize = 128;
 const max = 9;
 const steps = 5;
 const maxverts = 1000;  // keep from going off the rails
 for (let p = 3; p <= max; ++p) {
   for (let q = 3; q <= max; ++q) {
     await sleep(0);  // flush previously drawn stuff before computing next
     document.body.innerHTML += doit(svgsize, p, q, steps, maxverts);
   }
   document.body.innerHTML += "<br>";
 }

 const t1 = performance.now();
 console.log("        out script async function in "+(t1-t0)/1000+"s")
   })();
   console.log("    out script")
  </script>
 </body>
</html>