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:
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}"]
