import numpy as np

def compute_volume(D, N = 100000):
    vol = 0.
    for i in range(N):
        x = np.random.random_sample(D) * 2. - 1. # draw a d-dimensional vector with coordinates between -1 and 1 (in the cube containing the sphere).
        if np.sum(x ** 2) <= 1.: # check whether the vector is in the sphere.
            vol += 1.
    vol /= N
    vol *= 2 ** D # multiply by the probability density of the uniform distribution between -1 and 1 in d dimensions.
    return vol

nrepet = 20
volume_3d = 0.
volume_10d = 0.
for k in range(nrepet):
    volume_3d += compute_volume(3)
    volume_10d += compute_volume(10)
volume_3d /= nrepet
volume_10d /= nrepet
print(volume_3d, 4 * np.pi / 3.)
print(volume_10d, np.pi ** 5 / 120.)