Here is the same code as posted on MO. However, here i is the place to say some words about the sage part of the implementation.

```
UHP = HyperbolicPlane().UHP() # UHP is the upper half plane IH
HM = HyperbolicPlane().HM() # HM is the hyperboloid model
a1, a2, a3 = pi/4, pi/4, pi/4 # given angles, we draw a hyperbolic triangle with these angles
def c(a1, a2, a3):
return (cos(a1) + cos(a2)*cos(a3)) / sin(a2) / sin(a3)
c1, c2, c3 = c(a1, a2, a3), c(a2, a3, a1), c(a3, a1, a2) # algebraic in the given example
def s(v, w):
"""v, w are vectors with three entries, we return the Minkowski product with signature ++-"""
return v*diagonal_matrix([1, 1, -1])*w
# a, b, p, q, r are used in the following coordinates in Hyperboloid model
# as a parametrization of points
myvars = var("a b p q r");
a, b, p, q, r = myvars
V1 = vector([0, 0, 1])
V2 = vector([0, a, b])
V3 = vector([p, q, r])
sols = solve([ s(V2, V2) == -1, s(V3, V3) == -1,
s(V1, V2) == -c3, s(V2, V3) == -c1, s(V3, V1) == -c2 ]
, myvars, solution_dict=True)
sols = [sol for sol in sols if sol[a] > 0 and sol[q] > 0] # so V2, V3 maps to IH
sol = sols[0] # first solution
a0, b0, p0, q0, r0 = [sol[v].simplify_full() for v in myvars]
S1, S2, S3 = vector([0, 0, 1]), vector([0, a0, b0]), vector([p0, q0, r0])
M1, M2, M3 = HM.get_point(S1), HM.get_point(S2), HM.get_point(S3)
H1, H2, H3 = UHP(M1), UHP(M2), UHP(M3) # using the coercion from HM to UHP
Q1, Q2, Q3 = H1.coordinates(), H2.coordinates(), H3.coordinates()
p = hyperbolic_polygon(pts=[Q1, Q2, Q3], model="UHP", fill=True, alpha=0.3)
g = Graphics()
g += p.plot()
g.show(axes=True, aspect_ratio=1)
```

It was possible to perform "exact arithmetics". Sage comes with some models of the `HyperbolicPlane`

- and here they are:

```
sage: HP = HyperbolicPlane()
sage: HP.HM()
Hyperbolic plane in the Hyperboloid Model
sage: HP.Hyperboloid()
Hyperbolic plane in the Hyperboloid Model
sage: # same
sage: HP.KM()
Hyperbolic plane in the Klein Disk Model
sage: HP.KleinDisk()
Hyperbolic plane in the Klein Disk Model
sage: # same
sage: HP.PoincareDisk()
Hyperbolic plane in the Poincare Disk Model
sage: HP.PD()
Hyperbolic plane in the Poincare Disk Model
sage: # same
sage: HP.UHP()
Hyperbolic plane in the Upper Half Plane Model
sage: HP.realizations()
[Hyperbolic plane in the Upper Half Plane Model,
Hyperbolic plane in the Hyperboloid Model,
Hyperbolic plane in the Klein Disk Model,
Hyperbolic plane in the Poincare Disk Model]
```

In the implementation, we start with the given angles, use the hyperbolic laws of sine and cosine, see also the short pdf, solve a system in the `HM`

, then move it in the world of the upper half plane $\Bbb H$ by simple coercion, and there is for $\Bbb H=$`UHP`

a coercion from `HM`

:

```
UHP = HyperbolicPlane().UHP() # UHP is the upper half plane IH
HM = HyperbolicPlane().HM() # HM is the hyperboloid model
```

and now:

```
sage: UHP.has_coerce_map_from(HM)
True
```

After solving the system, your idea to use this framework, we obtain the vectors `V1, V2, V3`

, each with three coordinates, that should be considered now in `HM`

. It turns out that we have the right order of components.

For instance the solution `S2`

(instead of `V2`

) is:

```
sage: S2
(0, sqrt(2)*sqrt(sqrt(2) + 1), sqrt(2) + 1)
sage: s(S2, S2)
-(sqrt(2) + 1)^2 + 2*sqrt(2) + 2
sage: _.simplify_full()
-1
```

And the line is accepted:

```
sage: HM(S2)
Point in HM (0, sqrt(2)*sqrt(sqrt(2) + 1), sqrt(2) + 1)
```

so the bigger coordinate is at the last place. The point `HM(S2)`

becomes now `M2`

. Similarly, we have after solving the system and mapping the vectors to `HM`

also the other points `M1, M3`

. Sage gives us the chance to compute the distance between these points in the given model:

```
sage: # HM.dist(M1, M2) # arccosh( a long expression )
sage: HM.dist(M1, M2).n()
1.52857091948100
sage: arccosh(c3).n()
1.52857091948100
```

So the distance corresponds to the one stated in the linked pdf, second cosine rule in hyperbolic geometry. The distance between each two points is the same. For instance, we have symbolically the `True`

in the lines:

```
sage: bool( HM.dist(M1, M2) == HM.dist(M2, M3) )
True
sage: bool( HM.dist(M2, M3) == HM.dist(M3, M1) )
True
```

From here, we can pass to the upper half-plane $\Bbb H=$`UHP`

, and the points are also exact, to work with them as complex number we take the "coordinates", and we can convert them now to numbers in some exact field we prefer, for instance in $\bar{\Bbb Q}$ (or with some work in some cyclotomic field). Here are some test lines, that illustrate the structure, using

```
S1, S2, S3 = vector([0, 0, 1]), vector([0, a0, b0]), vector([p0, q0, r0])
M1, M2, M3 = HM.get_point(S1), HM.get_point(S2), HM.get_point(S3)
H1, H2, H3 = UHP(M1), UHP(M2), UHP(M3) # using the coercion from HM to UHP
Q1, Q2, Q3 = H1.coordinates(), H2.coordinates(), H3.coordinates()
x1, x2, x3 = QQbar(Q1), QQbar(Q2), QQbar(Q3)
```

we have:

```
sage: UHP.dist(H1, H2).n()
1.52857091948100
sage: arccosh(c3).n() # same distance in HM and UHP
1.52857091948100
sage: x1, x2, x3
(I, 4.611581789308715?*I, -1.805790894654358? + 1.162196641748775?*I)
sage: x1.minpoly(), x2.minpoly(), x3.minpoly()
(x^2 + 1, x^8 + 20*x^6 - 26*x^4 + 20*x^2 + 1, x^8 - 4*x^6 + 22*x^4 - 4*x^2 + 1)
sage: x2.minpoly().galois_group()
Transitive group number 4 of degree 8
sage: x2.minpoly().galois_group().structure_description()
'D4'
sage: x3.minpoly().galois_group()
Transitive group number 4 of degree 8
sage: x3.minpoly().galois_group().structure_description()
'D4'
```

This may seem less important for drawing one triangle, but may become important when trying to draw a full tesselation starting from one triangle and applying Möbius transformations, details in a deep level can be made visible (also in an animation with zoom-in features)...

FYI: Although I've implemented the following algorithm, it didn't work well so I'm asking for help as well. https://mathoverflow.net/questions/46...