Cvičení č. 2: Počítáme π

V našem programování pokročíme počítáním podle algoritmu pro součet aritmetické posloupnosti, kterému jsme se věnovali na přednášce. Analogickým postupem získáme součty několika dalších řad, tentokrát s líbivým výsledkem (souhrn vzorců v PDF). Pracujeme v imperativním strukturovaném stylu; pro inspiraci si také naznačíme, jak by se naše úlohy řešily i v dnes oblíbeném funkcionálním stylu programování.

Machinovy vzorce

Budeme počítat \(\pi\). Některé jazyky ho nabízejí jako předdefinovanou konstantu, jiné ho nevyzdvihují a člověk si ho musí připravit sám. Lépe než pamatovat si či dohledávat jeho 15 či 30 desetinných míst v 8- či 16bytové přesnosti je zapamatovat si navždy vzorec, na který si člověk vzpomene i po probuzení z hlubokého spánku:

\[\pi=4\arctan1\]

Obdobné vzorce s funkcí \(\arctan\) nesou jméno Johna Machina (viz např. odvození posledního z nich):

\[\begin{split}\frac\pi4=\arctan\frac12+\arctan\frac13=2\arctan\frac12-\arctan\frac17=2\arctan\frac13+\arctan\frac17 \\ \frac\pi4=4\arctan\frac15-\arctan\frac1{239}\end{split}\]

Program, zatím pouhou sekvenci příkazů, zapíšeme v Pythonu a Fortranu:

# Python: předdefinované pí a Machinovo pí
from math import *                                      # import všech jmen z modulu math
format='%25s%20.15f'                                    # formátová specifikace
print(format%('pi podle math =',pi))                    # printf-style string formatting
# print('{:>25}{:20.15f}'.format('pi podle math =',pi)) # string format() method
# print(f'{"pi podle math =":>25}{pi:20.15f}')          # formatted string literals
piMachin=(4*atan(1/5)-atan(1/239))*4                    # reálné dělení, (1/5) je 0.2
print(format%('pi podle Machina =',piMachin))

Pro předdefinované pi je třeba jít do modulu math, stejně jako pro funkci atan (arkus tangens). Operátor / v Pythonu stojí pro reálné dělení, jeho operandy můžeme ponechat v celočíselném typu int. Pro výpis jsme naznačili tři pythonské cesty, jak formátovat výstup. Všechny zahrnují explicitní pokyn pro zarovnání řetězce do 25 znaků a reálného čísla do 20 znaků s 15 desetinnými místy.

! Fortran: půlnoční pí a Machinovo pí
program p2a
implicit none                            ! zákaz implicitních popisů jmen
real(8) pi,piMachin                      ! explicitní popis real(8) proměnných
character(100) :: format='(a25,f20.15)'  ! formátová specifikace

! kdyby vás o půlnoci vzbudili...
pi=atan(1._8)*4                          ! argument musí být typu real, raději real(8)
print format,'pulnocni pi =',pi

! Machin
piMachin=(4*atan(1._8/5)-atan(1._8/239))*4  ! pozor na integer dělení, (1/5) by byla 0
print format,'pi podle Machina =',piMachin
end program

Fortran předdefinované pi nemá, tedy ho počítáme. Fortranské funkce typicky vyhodnocují v přesnosti argumentu, tedy pro real(8) přesnost funkce atan vložíme do argumentu 8bytovou reálnou jedničku 1._8. Fortranský operátor / se rozhoduje pro celočíselné nebo reálné dělení podle typu operandů; je-li alespoň jeden reálný, dělí se v reálné aritmetice. Výstup formátujeme obdobně jako v pythonské ukázce explicitně.

Leibnizův a Eulerův součet

Na přednášce jsme obecný sumační algoritmus zapsali ve strukturovaném stylu pomocí cyklů. Indexovaný cyklus z přednášky použijeme jako vzor pro výpočet pi součtem Leibnizovy řady,

\[\frac\pi4=\sum_{n=1}^\infty\frac{(-1)^{n-1}}{2n-1}=1-\frac13+\frac15-\frac17+\cdots\]

a Eulerovy řady,

\[\frac{\pi^2}6=\sum_{n=1}^\infty\frac1{n^2}=1+\frac1{2^2}+\frac1{3^2}+\frac1{4^2}+\cdots\]

Pythonský indexovaný cyklus vyjádříme pomocí funkce range(start,stop,step), vracející aritmetickou posloupnost celých čísel z polouzavřeného intervalu <start,stop) s nenulovým krokem step:

# Python a imperativní pí: Leibnizova a Eulerova řada

from math import *
print(f'{"pi podle math =":>25}{pi:20.15f}')

# Leibniz
nmax=1000000
piLeibniz=0.
z=1.
for n in range(1,nmax+1):
  # piLeibniz=piLeibniz+(-1)**(n-1)/(2*n-1)  # pomalé umocňování
  piLeibniz=piLeibniz+z/(2*n-1)              # alternace znaménka pomocí proměnné z
  z=-z
piLeibniz=piLeibniz*4
print(f'{"pi podle Leibnize =":>25}{piLeibniz:20.15f}')

# Euler: kumulace od největších členů
nmax=1000000
piEuler=0.
for n in range(1,nmax+1):
  piEuler=piEuler+1/(n*n)     # dělení součinem
  # piEuler=piEuler+1/n/n     # dvě pomalá dělení
  # piEuler=piEuler+1/n**2    # dělení pomalou mocninou
piEuler=sqrt(piEuler*6)
print(f'{"pi podle Eulera =":>25}{piEuler:20.15f}')

# Euler: kumulace od nejmenších členů
nmax=1000000
piEuler=0.
for n in range(nmax,0,-1):    # cyklus od nmax do 1
  piEuler=piEuler+1/(n*n)
piEuler=sqrt(piEuler*6)
print(f'{"pi podle Eulera =":>25}{piEuler:20.15f}')

Na cvičení jsme se sledováním doby běhu přesvědčili, jak pomalé může být numerické umocňování, a nahradili umocňování (-1)**(n-1), v podstatě matematický symbol pro alternaci znaménka, explicitním alternováním znaménka jedničky v proměnné z. Sledovali jsme rovněž rychlost konvergence; z absolutní hodnoty n-tého členu se zdá, že s každým řádem n přibude ve výsledku jedna platná číslice. Podobnou rychlost konvergence bychom mohli čekat i od Eulerovy metody. Shledali jsme oproti očekávání jistou nepřesnost. Pro její příčinu jsme ukázali na směr sumace řady: při sumování od největších hodnot se začnou vlivem zaokrouhlování mezisoučty kazit dříve, než při sčítání od nejmenších hodnot. (Tento kaz se projeví až při nmax řádu miliard, do té doby je pomalá konvergence zdrojem větší chyby. K miliardě se ovšem v Pythonu bez dalšího nepropočteme, ukáže to až Fortran.) I zde jsme testovali několik variant zápisu členu řady – z variant a) jedno dělení součinem, b) dvě dělení, c) dělení mocninou jsme zde vyhodnotili jako nejrychlejší variantu a).

Fortranský indexovaný cyklus má podobu konstrukce do index=start,stop,step; ...; end do, jejíž tělo se vykoná pro index z uzavřeného intervalu <start,stop> s nenulovým krokem step:

! Fortran a imperativní pí: Leibnizova a Eulerova řada
program p2b
implicit none
real(8) pi,piLeibniz,piEuler
real(8) z,rn
integer nmax,n
character(100) :: format='(a25,f20.15)'

! kdyby vás o půlnoci vzbudili...
pi=atan(1._8)*4
print format,'pulnocni pi =',pi

! Leibniz
nmax=1000000
piLeibniz=0.
z=1.
do n=1,nmax
! piLeibniz=piLeibniz+(-1._8)**(n-1)/(2*n-1)  ! pomalé umocňování
  piLeibniz=piLeibniz+z/(2*n-1)               ! alternace znaménka pomocí proměnné z
  z=-z
end do
piLeibniz=piLeibniz*4
print format,'pi podle Leibnize =',piLeibniz

! Euler: kumulace od největších členů
nmax=1000000
piEuler=0.
do n=1,nmax
  ! piEuler=piEuler+(1._8)/(n*n)    ! integer overflow
  ! piEuler=piEuler+(1._8)/n/n      ! dvě pomalá dělení
  ! piEuler=piEuler+1/real(n,8)**2  ! dělení reálnou mocninou
  rn=n; piEuler=piEuler+1/(rn*rn)   ! dělení reálným součinem
end do
piEuler=sqrt(piEuler*6)
print format,'pi podle Eulera =',piEuler

! Euler: kumulace od nejmenších členů
nmax=1000000
piEuler=0.
do n=nmax,1,-1                      ! cyklus od nmax do 1
  rn=n; piEuler=piEuler+1/(rn*rn)
end do
piEuler=sqrt(piEuler*6)
print format,'pi podle Eulera =',piEuler

end program

V Eulerově součtu se 4bytovým typem integer hrozí ve variantě s (1._8)/(n*n) integer overflow, přešli jsme proto k vyčíslení pomocí reálné proměnné rn. Ukázalo se také, že pokud překládáme s optimalizací, gfortran -O, překladač smysl výrazu (-1._8)**(n-1) rozezná, mocninu nahradí jiným řešením a rychlost běhu zrychlí na rychlost naší alternativy.

Viètovy odmocniny

Viètův nekonečný součin

\[\frac2\pi=\sqrt{\frac12}\sqrt{\frac12+\frac12\sqrt{\frac12}} \sqrt{\frac12+\frac12\sqrt{\frac12+\frac12\sqrt{\frac12}}}\cdots\]

bude potřebovat vyjádření vývoje \(n\)-tého členu součinu: \(a_0=0,\) \(a_n=\sqrt{\frac12+\frac12a_{n-1}},\) \(n=1,\dots\) Stačit k tomu bude update jedné (skalární, neindexované) proměnné: a=sqrt(0.5+0.5*a).

Python:

# Python a imperativní pí: Vietův součin

from math import *
print(f'{"pi podle math =":>25}{pi:20.15f}')

# Viete
nmax=24              # rychlá konvergence
piViete=1.
a=0.
for n in range(1,nmax+1):
  a=sqrt(0.5+0.5*a)  # update členu řady
  piViete=piViete*a  # dílčí součin
piViete=2/piViete
print(f'{"pi podle Vieta =":>25}{piViete:20.15f}')

Fortran:

! Fortran a imperativní pí: Vietův součin
program p2c
implicit none
real(8) pi,piViete
real(8) a
integer nmax,n
character(100) :: format='(a25,f20.15)'

! kdyby vás o půlnoci vzbudili...
pi=atan(1._8)*4
print format,'pulnocni pi =',pi

! Viete
nmax=24              ! rychlá konvergence
piViete=1.
a=0.
do n=1,nmax
  a=sqrt(0.5+0.5*a)  ! update členů řady
  piViete=piViete*a  ! dílčí součin
end do
piViete=2/piViete
print format,'pi podle Vieta =',piViete

end program

Povšimli jsme si také rychlé konvergence Vietova součinu, k dosažení 15 platných míst stačilo 25 členů.

DÚ: Wallisův součin

Za domácí úkol zkuste doplnit náš zdrojový text (v Pythonu, ve Fortranu nebo v obou) tak, aby počítal \(\pi\) i aproximací Wallisova nekonečného součinu:

\[\frac\pi2=\prod_{n=1}^\infty\frac{(2n)^2}{(2n-1)(2n+1)} =\frac{2\cdot2}{1\cdot3}\cdot \frac{4\cdot4}{3\cdot5}\cdot\frac{6\cdot6}{5\cdot7}\cdot \frac{8\cdot8}{7\cdot9}\cdots\]

Py a f90 zdrojové soubory s dnešními ukázkami: zip. Náš souhrn vzorců: PDF. Více návodů, jak počítat \(\pi\): MathWorld, Wikipedia, Piphilology.