Plot a hyperbolic plane tessellation 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.
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]]) &lt; 0., a = a + 2 \[Pi]];
  If[(b = N[Apply[ArcTan, {Qx, Qy} - c]]) &lt; 0., 
   b = b + 2 \[Pi]]; {a, b} = Sort[{a, b}];
  If[b - a &gt; \[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], #] &amp;, p, {2}]

Inversion[Circle[{Cx_, Cy_}, r_, {a_, b_}], p_Polygon] := 
 Map[Inversion[Circle[{Cx, Cy}, r], #] &amp;, 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}}], #] &amp;, p, {2}]

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


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

LineRule = Polygon[x_] :&gt; Line[Join[x, {x[[1]]}]];
HyperbolicLineRule = 
  Polygon[x_] :&gt; 
   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]]}}.# &amp;,
     Transpose[{Cos[\[Theta]], Sin[\[Theta]]}]]]

PolygonUnion[p_Polygon, tol_: 10.^-10] := p
PolygonUnion[p_List, tol_: 10.^-10] := 
 With[{q = p /. Polygon[x_] :&gt; N[Polygon[Round[x, 10.^-10]]]}, 
  DeleteDuplicates[q]]
HyperbolicTessellation[p_Integer, q_Integer, \[Phi]_, k_Integer, 
  t_: 10.^-10] := 
 Map[PolygonUnion[#, t] &amp;, 
   NestList[PolygonInvert, Polygon[CentralPolygon[p, q, \[Phi]]], 
     k][[{-2, -1}]]] /; k &gt; 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[3, 7, 0., 7, HyperbolicLineRule, ImageSize -&gt; 300, PlotLabel -&gt; &quot;{7,7}&quot;]