Как приближенно вычислить значение числа Пи с использованием объема сферы?

Приблизим осьмушку шара треугольным мешем. Вершины меша будут лежать на поверхности единичной сферы. Тело ограниченное координатными плоскостями и мешем лежит целиком внутри осьмушки шара. Его объём используем для оценки pi снизу. Отыщем множитель на который нужно умножить координаты всех вершин меша так чтобы после умножения фигура содержала осьмушку шара внутри себя. Объём новой (раздутой) фигуры используем для оценки pi сверху. Так как мы не знаем заранее сколько должно быть треугольников в меше, то будем повторять процедуру со всё более подробными мешами пока оценки pi не сблизятся на расстояние 10^-6. Поехали.

Четыре поколения поверхностей:

Начинается всё с одного треугольника. Следующая поверхность строится так: каждый треугольник предыдущей разбивается на четыре меньшего размера. Новые точки на серединах сторон оказываются под поверхностью сферы. Чтобы вернуть их на сферу их координаты нормализуются.

Объём полигонального тела складывается из объёмов тетраэдров. Тетраэдр строится на треугольнике поверхности и точке (0, 0, 0). Шестикратный объём такого тетраэдра вычисляется как определитель.

Имеем последовательность оценок числа pi снизу. Как оценить pi сверху?

Каждый треугольник поверхности определяет плоскость. Вычислим минимальное расстояние от этой плоскости до начала координат: строим единичный нормальный вектор к треугольнику, скалярно умножаем его на координаты любой вершины треугольника. Расстояние получится меньше единицы.

Если все координаты всех вершин треугольника умножить на величину обратную к найденному расстоянию, то новые вершины образуют новую плоскость (параллельную старой), для которой расстояние будет единичным.

Выберем минимум из расстояний для всех треугольников поверхности. Все вершины поверхности домножим на обратную величину. Получим новую поверхность гомотетичную старой. Новая поверхность целиком проходит вне шара. То есть объём нового тела будет больше pi. Конечно, суммировать его заново не надо, он вычисляется из объёма поверхности внутри шара умножением на некоторый коэффициент.

Код:

import math


def add(p1, p2):
    x1, y1, z1 = p1
    x2, y2, z2 = p2
    return x1 + x2, y1 + y2, z1 + z2


def sub(p1, p2):
    x1, y1, z1 = p1
    x2, y2, z2 = p2
    return x1 - x2, y1 - y2, z1 - z2


def scale(f, p):
    x, y, z = p
    return f * x, f * y, f * z


def mid(p1, p2):
    return scale(1 / 2, add(p1, p2))


def dot(p1, p2):
    x1, y1, z1 = p1
    x2, y2, z2 = p2
    return x1 * x2 + y1 * y2 + z1 * z2


def norm(p):
    return math.sqrt(dot(p, p))


def normalized(p):
    return scale(1 / norm(p), p)


def cross(p1, p2):
    x1, y1, z1 = p1
    x2, y2, z2 = p2
    return \
        y1 * z2 - z1 * y2, \
        z1 * x2 - x1 * z2, \
        x1 * y2 - y1 * x2


def normal(p1, p2, p3):
    return normalized(cross(sub(p2, p1), sub(p3, p1)))


def height(p1, p2, p3):
    return dot(normal(p1, p2, p3), p1)


def volume6(p1, p2, p3):
    x1, y1, z1 = p1
    x2, y2, z2 = p2
    x3, y3, z3 = p3
    # | x1, y1, z1 |
    # | x2, y2, z2 |
    # | x3, y3, z3 |
    return \
        ((x1 * y2) - (x2 * y1)) * z3 - \
        ((x1 * y3) - (x3 * y1)) * z2 + \
        ((x2 * y3) - (x3 * y2)) * z1


def subdiv(p1, p2, p3):
    p12 = normalized(mid(p1, p2))
    p23 = normalized(mid(p2, p3))
    p31 = normalized(mid(p3, p1))
    yield p1, p12, p31
    yield p2, p23, p12
    yield p3, p31, p23
    yield p12, p23, p31


def main():
    surface = (((1, 0, 0), (0, 1, 0), (0, 0, 1)), )
    for i in range(1, 1000):
        print(f'{i}({len(surface)}): ', flush=True, end='')
        min_height = min(height(*t) for t in surface)
        v_low = sum(volume6(*t) for t in surface)
        v_high = v_low * (1 / min_height) ** 3
        print(f'{v_low} < pi < {v_high}, {v_high - v_low}')
        if v_high - v_low < 1e-6:
            break
        surface = tuple(tt for t in surface for tt in subdiv(*t))


main()

$ time python pi_approximation.py 
1(1): 1 < pi < 5.196152422706629, 4.196152422706629
2(4): 2.207106781186547 < pi < 4.054714066307733, 1.8476072851211862
3(16): 2.8632972148712703 < pi < 3.4166084802548, 0.5533112653835297
4(64): 3.0687004650644996 < pi < 3.2136638054952016, 0.14496334043070203
5(256): 3.123152697186476 < pi < 3.1598235333433777, 0.03667083615690192
6(1024): 3.136968953405614 < pi < 3.1461637781466827, 0.009194824741068697
7(4096): 3.140435869168335 < pi < 3.142736273849948, 0.0023004046816126333
8(16384): 3.1413034037353573 < pi < 3.141878611120626, 0.0005752073852685058
9(65536): 3.141520337766175 < pi < 3.1416641462518413, 0.00014380848566641902
10(262144): 3.1415745744239363 < pi < 3.141610526960324, 3.595253638755125e-05
11(1048576): 3.1415881337850373 < pi < 3.141597121945071, 8.988160033585046e-06
12(4194304): 3.141591523637882 < pi < 3.141593770679511, 2.2470416292108553e-06
13(16777216): 3.141592371100547 < pi < 3.1415929328610557, 5.617605087770983e-07

real  2m59.102s
user  2m57.320s
sys   0m1.480s