Hyperbolic Tessellations of a Pentagon in Mathematica

I want to implement this in Mathematica:

I realize I can use ParametricPlot to get the lines, and then use the Mesh option to fill in the colors. Can anyone tell me what a general formula for the equations of these lines are? Does it generalize to a regular n-gon?

I wrote some code in 1999 that defines a PTiling function that does something close to what you want. I've included that code at the bottom of this response. Once you execute that to define PTiling, you should be able to generate a picture like so:

Show[PTiling[5, 2, 3], ImageSize -> 600, PlotRange -> 1.1]

Obviously, we have just a black and white picture illustrating the boundaries of the pentagons that you want but, hopefully, this is a good start. I think that one could use portions of disks, rather than circles, to get closer to what you want.

A few of comments on the code are in order. The ideas are all described in Saul Stahl's excellent book, The Poincare Half-Plane - specifically, the chapter on the Poincare disk. I wrote the code to illustrate some ideas for a geometry class that I was teaching in 1999 so it must have been for version 3 or 4. I have done nothing to try to optimize the code for any subsequent version.

Also, it's not the case that all combinations of n and k work so well. I think they actually need to be chosen somewhat carefully; I honestly don't recall exactly how they should be chosen.

(* Generates a groovy image of the Poincare disk *)
(* Not all combinations of n and k work great *)

PTiling[n_Integer, k_Integer, depth_ : 2] := Module[
    {aH, a, q, r, idealPoints, z1Ideal, z2Ideal, init, initCircs},
    
    aH = ArcCosh[(Cos[Pi/n] Cos[Pi/2] + 
            Cos[(n - 2) Pi/(2 k)])/(Sin[Pi/n] Sin[Pi/2])];
    a = (Exp[aH] - 1)/(Exp[aH] + 1);
    q = (a + 1/a)/2;
    r = q - a;
    
    (* The Action *)
    idealPoints = {x, y} /. NSolve[{x^2 + y^2 == 1, 
            (x - q)^2 + y^2 == r^2}, {x, y}]; 
    {z1Ideal, z2Ideal} = toC /@ idealPoints;
    
    init = N@IdealPLine[{z1Ideal, z2Ideal}];
    initCircs = NestList[RotateCircle[#, 2 Pi/n] &, init, n - 1];
    
    Show[Graphics[{Nest[Iter, initCircs, depth], PGamma}],
        AspectRatio -> Automatic]
];


(* Ancillary code *)

(* One step in the iteration *)
Iter[PLineList_List] := Union[PLineList, Select[Flatten[
        Outer[Reflect, PLineList, PLineList]], (# =!= Null &)],
  SameTest -> sameCircleQ
];


(* Generate the ideal Poincare line through the points z1 and z2 *)
(* Should be an arc, if z1 and z2 are not on the same diameter *)
IdealPLine[{z1_, z2_}] := Module[
    {center},
    center = Exp[I (Arg[z2] + Arg[z1])/2] / 
            Cos[(Arg[z2] - Arg[z1])/2];
    arc[{z1, z2}, center]
];
arc[{z1_, z2_}, z0_] := Module[{theta1, theta2},
        theta1 = Arg[z1 - z0]; 
        theta2 = Arg[z2 - z0];
    
        If[Abs[N[theta1 - theta2]] > Pi, 
            If[N[theta1] < N[theta2], 
                theta1 = theta1 + 2 Pi, 
                theta2 = theta2 + 2 Pi]
    ];
    
        Circle[toR2[z0], Abs[z1 - z0],
            Sort[{theta1, theta2}, N[#1] < N[#2] &]
    ]
];

(* Circle operations *)
Invert[Circle[c_, r1_, ___], z_] := 
        r1^2/Conjugate[z - toC[c]] + toC[c];
Reflect[circ1_Circle, Circle[c2_, r2_, thetaList_]] := 
        IdealPLine[
      Invert[circ1, toC[c2 + r2 *{Cos[#], Sin[#]}]] & /@ thetaList
    ];
RotateCircle[Circle[c_, r_, thetaList_], theta_] := 
    Circle[RotationMatrix[theta] . c, r, theta + thetaList];

cSameCircleQ = Compile[
   {{c1, _Real, 1}, {r1, _Real}, {th1, _Real, 1},
    {c2, _Real, 1}, {r2, _Real}, {th2, _Real, 1}},
   (c1[[1]] - c2[[1]])^2 + (c1[[2]] - c2[[2]])^2 + (r1 - r2)^2 +
     (th1[[1]] - th2[[1]])^2 + (th1[[2]] - th2[[2]])^2 < 0.00001];
sameCircleQ[Circle[c1_, r1_, th1_], Circle[c2_, r2_, th2_]] := 
  cSameCircleQ[c1, r1, th1, c2, r2, th2];


(* Basics *)
toR2 = {Re[#], Im[#]} &;
toC = #[[1]] + #[[2]] I &;
PGamma = {Thickness[.008], GrayLevel[.3],Circle[{0, 0}, 1]};