Inaccurate volumes of Dirichlet domains

[NathanDunfield]

The volume of a DirichletDomain often does not match that of the associated Manifold. Here is a quite simple manifold where they only agree to three digits:

>>> M = snappy.Manifold('m004(1,2)')
>>> M.volume()
1.39850888415
>>> M.dirichlet_domain().volume()
1.393551248521783
>>> H = M.high_precision(); H.volume()
1.398508884150806640509594326663004920276433947325827118548221874
>>> H.dirichlet_domain().volume()
1.393551248521774626126463081235516657883615640326408683944746082

The problem is not with the location of the vertices of the domain. I did an alternate computation of the volume from the domain which gave the correct answer to 1e-62 when working at quad-double precision: I put a vertex at the barycenter of each face, coned out from that to divide each face into triangles, and then coned each of those triangles to the center, and used the general volume formula for a hyperbolic tetrahedron involving dilogarithms.

In contrast, the SnapPy kernel uses a decomposition of the Dirichlet domain into birectangular tets (possibly to simplify the volume formula) which can involve tetrahedra with both positive and negative volume, possibly creating some sort of accumulating numerical error.

[NathanDunfield]

Root problem is two-fold:

1. The algorithm used for computing the volume of birectangular tets involves inverting a matrix. If the matrix is close to singular, the algorithm decides that the tetrahedra is degenerate and returns the volume of 0. However, the threshold for calling the matrix singular is too stringent, resulting in some tets with volume as much as 0.0004 being viewed as degenerate.

2. For quad-double precision, the implementation of the Lobachevsky function used is the same as for double precision and lacks sufficient terms in the Taylor expansion.

The problem was largely fixed in 20debb0c5ac9b3059efaf2ca0ad64ff6c7728ecf:

>>> M = snappy.Manifold('m004(1,2)')
>>> M.dirichlet_domain().volume() - M.volume()
-2.15494289079743e-12
>>> H = M.high_precision()
>>> H.dirichlet_domain().volume() - H.volume()
-1.41506748774410824854345245597920110295062521238937867865487793e-60

though there are still some cases where the Dirichlet volume is not quite as accurate as one might expect.