import numpy as np
from scipy.integrate import quad
from scipy.special import ive

def G_func(point):
    x, y, z = np.abs(point)
    # The integrand is stable using ive (exponentially scaled Bessel)
    def integrand(t):
        return ive(x, t/3) * ive(y, t/3) * ive(z, t/3)

    # Correct integration without the extra factor of 3
    res, _ = quad(integrand, 0, np.inf, epsabs=1e-13, epsrel=1e-13)
    return res

def ever_hit_probability(point):
    # Precomputed G(0,0,0) for the 3D cubic lattice
    G0 = 1.5163860591519780
    
    if point == (0, 0, 0):
        return 1 - (1 / G0)
    
    return G_func(point) / G0

n=4
# Test results
dct = {}
for p in range(n):
    for q in range(n):
        for r in range(n):
            dct[(p, q, r)] = ever_hit_probability((p, q, r))
print(dct)