Output image from example code

Attachment: geodesic_bug.png

The current implementation of hyperbolic geodesics
relies on a conversion to the half-plane model.

Any geodesic g has its corresponding
half-plane model geodesic stored as
g._cached_geodesic.

The example geodesic has endpoints
I and 2.222222222222222e-14 + 5.666666666666666*I.

The corresponding complete geodesic is a
euclidean semicircle. Due to rounding errors
on its centre and radius, its endpoints are
computed in CDF as (-2.375, 1400000000000002.0),
whereas computations in CC give them much more
correctly as (0.0, 1.4e15).

Converting back to the disc model, the effect
is a disaster.

Aha! I suspect it's possible to find the ideal endpoints of a disk model geodesic in a more numerically stable way. I can work on a test implementation.

If the test implementation works well, do you know what we'd have to do to incorporate it into hyperbolic_geodesic.py? Is it just a matter of overriding ideal_endpoints in HyperbolicGeodesicPD?

Replying to @Vectornaut:

Aha! I suspect it's possible to find the ideal endpoints
of a disk model geodesic in a more numerically stable way.
I can work on a test implementation.

Good!

If the test implementation works well, do you know
what we'd have to do to incorporate it into
hyperbolic_geodesic.py? Is it just a matter of
overriding ideal_endpoints in HyperbolicGeodesicPD?

Yes.

Changed keywords from none to hyperbolic, geodesic, plot, precision

Finding ideal endpoints of geodesics on the Poincaré disk: example implementation

Attachment: hyperbolic_geodesic_PD.sage.gz

Attachment: drawing-geodesics.pdf.gz

Finding ideal endpoints of geodesics on the Poincaré disk: method explanation

I just attached an example implementation. I've done some basic tests, but no systematic challenges or comparison with the current implementation (via the upper half-plane).

Maybe unrelated or not (If you think it's unrelated I can open a new ticket on this) I've just discovered working on ticket [ticket:23427] that ideal endpoints on PD geodesics can be mapped to imag()<0 points which is an error, and when plotted often draw a geodesic int the lower half plane.
For example

Sage:    P = PD.get_point(-0.920704875684763 - 0.390259569889459*I)
Sage:    P
Boundary point in PD -0.920704875684763 - 0.390259569889459*I
Sage:    UHP(P)
Boundary point in UHP -0.662253938491478 - 1.66533453693773e-16*I

sage: g=HyperbolicPlane().PD().random_geodesic()
Sage:  g
Geodesic in PD from -0.920704875684763 - 0.390259569889459*I to 0.0632992206159446 + 0.997994593507106*I
sage: g.to_model('UHP')  
Geodesic in UHP from -0.662253938491478 - 5.55111512312578e-17*I to 31.5642842686728 - 8.34887714518118e-14*I

Not being aware of this discussion I reported a similar issue here.

The one line change I mentioned there seems to me to fix the problem both in my use cases and in the case reported in the top post above. Maybe avoiding testing of exact equality of floting point numbers in the UHP code would be a good thing to do to improve the situation?

I found a case with similar symptoms, but which I think reflects a separate bug of the same type:

r = exp((pi*I/2).n())
G = hyperbolic_triangle(0, r, r*(1 + I)/2, model='PD', fill=True, alpha=.3)
G._show_axes = False
show(G)

(The central angle should be pi/4.) I think this is the root cause

PD = HyperbolicPlane().PD()
UHP = HyperbolicPlane().UHP()
r = exp((pi*I/2).n())
p = PD.get_point(r)
UHP(p)

which yields Boundary point in UHP 3.26624787063907e16 - 1.00000000000000*I A fix might be to change this line to

    if abs(x - I) < EPSILON:

(and add from sage.geometry.hyperbolic_space.hyperbolic_constants import EPSILON to that file's preamble).