Hi,
There is a bug in hyperbolic_polygon which happens quite rarely, but for example
    l = [-1, 1, -0.0571909584179366 + 0.666666666666667*I]
    hyperbolic_polygon(l, model='PD')
results in an incorrectly drawn hyperbolic triangle. It seems to be a subtle numerical issue; if the last element of l is replaced with -0.0571909584179366 + 0.667*I everything looks fine.
Digging a bit further, there seems to be an underlying issue here:
    PD = HyperbolicPlane().PD()
    z0 = CC(-0.0571909584179366 + 0.666666666666667*I)
    z1 = CC(-1)
    PD.get_geodesic(z0, z1).ideal_endpoints()
yields the erroneous result
    [Boundary point in PD -9.52420782539595e-17 + 1.00000000000000*I,
     Boundary point in PD -0.800000000000000 + 0.600000000000000*I]
(the second endpoint should be -1). Because of this I am suspicious of lines 1267 - 1269 in src/sage/geometry/hyperbolic_space/hyperbolic_geodesic.py. There we have the code
    # We could also have a vertical line with two interior points
    if x1 == x2:
        return [M.get_point(x1), M.get_point(infinity)]
Maybe the if statement should instead be something like
    if abs(x1 - x2) < EPSILON:
Best regards,
Håkan Granath