HyperShpere

import numpy as np
import pylab
import time
from math import pi

RSp = 1     # радиус сферы
MaxD = 16   # максимальная размерность

# реальные объемы с графика
Vreal = [1, 2, pi, 4*pi/3, pi**2/2, 8*pi**2/15, pi**3/6, 16*pi**3/105, pi**4/24, 32*pi**4/945, pi**5/120, 64*pi**5/10395, pi**6/720, 128*pi**6/135135, pi**7/5040, 15*pi**7/118771, 0.23533063035889312]

# x - количество точек, D - число измерений
def calcVNumPy(x, D):
    Xi = np.random.uniform(0, 1, x * D).reshape(x, D)    # набор случайных точек - вектора из [0,1]
    if D > 1:
      Dist = np.sum(Xi**2, axis=1)**0.5
      PtsInSphere = np.sum(Dist < RSp)                     # число точек внутри сферы
    elif D == 1:
      PtsInSphere = x

    return 2**D * PtsInSphere / x                        # объем методом Монте Карло

data = []

for d in range(1, MaxD+1):           # d - число измерений
    startTime = time.time()          # время начало расчета
    n = m = 10**5
    Vn = calcVNumPy(n, d)
    while True:
          Vm = calcVNumPy(m, d)
          Vnm = (n * Vn + m * Vm) / (n + m)

          epsA = abs(Vn - Vreal[d])          # абсолютная погрешность
          epsR = epsA / Vreal[d]             # относительная погрешность

          s = epsA + abs(Vnm) * epsR
          if abs(Vnm - Vn) <= s:
            break
          n = n + m
          Vn = Vnm
    endTime = time.time()

    # число измерений, количество точек, абс.погр., относит.погр., время в мс
    print("Время расчета для {:2d}-d t = {:.2f}с".format(d, endTime-startTime), end='. ')
    print("Погрешности абс={:.7f}, относ={:.7f}".format(epsA, epsR))
    data.append([d, Vn, n]) # сохраняем результаты

D, V, P = np.array(data).T
# График 1: расчитанный объем шара в зависимости от размерности
fig, axs = pylab.subplots(2, 1, figsize=(10,20))
axs[0].set_title('Расчитанный (x) и реальный объем (o)', color='C0')
axs[0].plot(D, Vreal[1:MaxD+1], 'C0', marker='o')
axs[0].plot(D, V, 'C2', marker='x')
# График 2: количество точек в зависимости от размерности пространства
axs[1].set_title('Количество точек', color='C1')
axs[1].plot(D, P, marker='x')

pylab.show()