# Standalone ordinary Mahler measure N(t), using mpmath in Sage or Python.
# Pass an integer, a decimal string, or a callable such as
# lambda: 18-6*mp.sqrt(3), evaluated afresh at each working precision.
# Numerical convergence estimates and precision agreement are not certified bounds.
import mpmath as mp

def hessian_measure(argument, working=140):
    with mp.workdps(working):
        t = mp.mpf(argument() if callable(argument) else argument)
        if not mp.isfinite(t):
            raise ValueError('The argument must be finite and real.')
        k = mp.sign(t) * mp.root(abs(t), 3)
        def integrand(theta):
            z = (k*mp.exp(mp.j*theta/2) - mp.exp(3*mp.j*theta/2))/2
            return abs(mp.re(mp.acosh(z)))
        points = [mp.mpf(0), mp.pi]
        if -1 < t < 27:
            points.insert(1, mp.acos((k-1)/2))
        tolerance = mp.power(10, -working+20)
        def integrate(a, b, depth=0):
            value, error = mp.quad(integrand, [a, b],
                                   error=True, maxdegree=8)
            if mp.isfinite(value) and error <= tolerance*(b-a)/mp.pi:
                return value
            if depth == 24:
                raise ArithmeticError('Quadrature did not converge.')
            c = (a+b)/2
            return integrate(a, c, depth+1) + integrate(c, b, depth+1)
        pieces = []
        for a, b in zip(points, points[1:]):
            c = (a+b)/2
            pieces.extend([integrate(a, c), integrate(c, b)])
        return mp.fsum(pieces)/mp.pi

low = hessian_measure('24', 140)
high = hessian_measure('24', 190)
with mp.workdps(200):
    assert abs(low-high) < mp.mpf('1e-115')*max(1, abs(high))
    print(mp.nstr(high, 100))