# Python 3 with SciPy: a double-precision example for rational 0 <= n < 5.
# Eight significant digits are displayed; no rigorous error bound is claimed.
# Tighten tolerances and check convergence, especially near n = 5.
from fractions import Fraction
from math import isfinite, pi
from scipy.integrate import solve_ivp


def pressure(index, rtol=2e-13, atol=2e-15, start=1e-4, limit=1e5):
    n_exact = Fraction(index)
    if not 0 <= n_exact < 5:
        raise ValueError('Require 0 <= n < 5.')
    n = float(n_exact)
    if not 0 <= n < 5:
        raise ValueError('Index too close to an endpoint for double precision.')
    if not 0 < start < limit:
        raise ValueError('Require 0 < start < limit.')
    # Regular central series through xi^4 and its derivative.
    theta = 1 - start**2 / 6 + n * start**4 / 120
    slope = -start / 3 + n * start**3 / 30

    def rhs(xi, state):
        y, v = state
        # A zero-density extension only locates the surface event;
        # fractional powers are evaluated only on nonnegative arguments.
        density = 1.0 if n == 0 else max(y, 0.0)**n
        return [v, -2 * v / xi - density]

    def surface(xi, state):
        return state[0]

    surface.terminal = True
    surface.direction = -1
    sol = solve_ivp(rhs, (start, limit), [theta, slope], method='DOP853',
                    events=surface, rtol=rtol, atol=atol)
    if sol.status != 1:
        raise RuntimeError('Surface not reached: ' + sol.message)
    v_surface = sol.y_events[0][0, 1]
    result = 1 / (4 * pi * (n + 1) * v_surface**2)
    if not isfinite(result) or result <= 0:
        raise ArithmeticError('Invalid pressure constant.')
    return result


if __name__ == '__main__':
    print(f"{pressure('13/4'):.8g}")