о том, как (не) стоит считать π, и о магии им. Эйлера


1.
одна из первых идей, приходящих в голову человеку, изучавшему базовый анализ — разложить подходящую функцию в степенной ряд и что-то в него подставить

если взять arctg(x) = x - x³/3 + x⁵/5 - …, подставить¹ x=1, то получается замечательно простая формула Лейбница π/4 = 1 - 1/3 + 1/5 - 1/7 + … — чего еще желать?..

¹ и проговорить нужные слова, потому что мы оказывается, вообще-то, не внутри, а на границе круга сходимости

в прошлый раз для сотни правильных знаков хватило 4 итераций, сейчас попробуем взять… ну, скажем, 1000 членов:

from mpmath import *

n = 1000
mypi = mpf(0)
s = mpf(4)
for k in range(n):
    mypi += s / mpf(2*k+1)
    s = -s

print(n,nstr(mypi,10),
      "diff:",nstr(pi-mypi,3,min_fixed=1))

1000 3.140592654 diff: 1.0e-3

видно, что погрешность ~1/1000… и вообще, можно понять, что для N слагаемых такой способ дает погрешность ~1/N — то есть чтобы получить 100 знаков нам нужно сложить ~10^(100) слагаемых (!!!)


2.
кажется, что это совсем тупиковая ветвь… и все же, не будем унывать сразу

посмотрим, скажем, на сумму 5 млн членов

погрешность получается ожидаемая, порядка 1/N, т.е. шесть цифр после запятой верны, а седьмая уже нет (получается 3,1415924… вместо 3,1415926…)

но то, что видно дальше, иначе как магией не назовешь: например, следующая дюжина цифр снова правильная!

from mpmath import *
mp.dps = 100

n = 2500000
mypi = mpf(0)
for k in range(1,n+1):
    mypi += mpf(8) / ((4*k-3)*(4*k-1))

for (s1,s2) in zip(nstr(mypi,60),nstr(pi,60)):
    print(s1 if s1==s2 else s1+'\u0332',end='')
print()

вот что получается в результате суммирования (неправильные цифры подчеркнуты):

3.14159245358979323846464338327950278419716939938730582097494

и, на самом деле, такого количества слагаемых уже достаточно, чтобы найти не только первые 6, но и первые 60 знаков π


продолжение ниже