Programování pro fyziky 2023/24 – přednáška č. 2

Opakování

To podstatné z první přednášky bylo toto:

  • zkusili jsme použít různé programovací jazyky (Python, Fortran, C ad.) užitečné pro fyzika, vše z příkazového řádku (psychologický backup) a prostředí VS Code,

  • přiřazovacími příkazy jsme ukládali do jednoduchých (skalárních) proměnných, vyčíslovali jsme výrazy s užitím aritmetických operátorů a sdělovali výsledky pomocí příkazů pro výpis na obrazovku (samé příkazy, proto – imperativní programování),

  • pozorovali jsme vlastnosti dvou numerických datových typů (integer, real) v dílčích variantách.

Dnešní cíl

Dnes vezmeme numerický problém, popíšeme algoritmus (etymologie slova) pro jeho řešení a ten naprogramujeme pomocí postupů imperativního a strukturovaného programování. (A pro odstrašení i pomocí nestrukturovaného, špagetového programování.)

Budeme řešit matematickou úlohu ve dvou variantách:

\[S_n=\sum_{i=1}^n i\,,\qquad H_n=\sum_{i=1}^n 1/i\,,\]

tedy budeme počítat součet aritmetické posloupnosti a harmonická čísla.

Algoritmus

Pro výpočet součtu řady (harmonická čísla jsou analogická) můžeme algoritmus popsat takto:

  • zjisti parametr \(n\),

  • inicializuj nulami index \(i\) a proměnnou \(s\) pro kumulaci součtu \(S_n\),

  • otestuj, zda \(i<n\),

  • pokud ano, zvětši \(i\) o 1, přičti aktuální hodnotu \(i\) k mezisoučtu \(s\) a vrať se na předchozí bod (na test),

  • pokud ne, pokračuj na následující bod,

  • vypiš slovem, co vypisuješ, pak hodnotu parametru \(n\) a výsledek \(s\).

Tento slovní popis se někdy, zvláště v počáteční fázi výkladu algoritmů, graficky znázorňuje pomocí vývojových diagramů. Pěkný diagram zachycující podmínky, vnořené podmínky a cyklus s podmínkou je zde. Jak je však patrné, už i jen minimalistické algoritmy o dvou podmínkách a jednom cyklu zaberou celou obrazovku, záhy i papír, a v praxi si s těmito diagramy nikdo nehraje. Vlastně si je člověk může jen tiše představovat a ještě lépe udělá, když si je v mysli zkondenzuje do konstrukcí svého programovacího jazyka a algoritmy si představuje rovnou tak. Což je samozřejmě třeba trénovat, a to tady hodláme dělat.

Pojďme tedy k ukázkám a těšme se na zápis algoritmu pro naši úlohu v programovacích jazycích. Ovšem – a to pokud možno před programováním vždy, ještě trochu přibrzděme. Je tak běžné, že problém rozumně formulovaný a řešený na papíře má na počítači jiné řešení, matematicky ekvivalentní, ale na počítači efektivnější. Zde jde o to, že při hodnotě parametru \(n\) bude výše popsaný algoritmus po počítači chtít přibližně \(3n\) elementárních (matematických či logických) operací: inkrementaci \(i\), update \(s\) a test \(i<n\). Je zvykem to shrnout jako závislost časové náročnosti algoritmu na vstupním parametru, zde: \(O(n)\), čímž rozumíme lineární časovou složitost algoritmu, počet operací úměrný parametru \(n\). Algoritmy mívají časovou složitost \(O(n)\), \(O(n\log n)\), \(O(n^2)\), \(O(n^3)\) apod. Lineární časová složitost problému je výborná, s \(n\) třeba u miliardy čekáme zhruba miliardu operací a miliardu elementárních operací provede procesor o GHz frekvenci řádově za sekundy.

Ale může to být ještě lepší. Pro součet aritmetické posloupnosti přece existuje součtový vzorec, který v našem případě nabude tvaru \(S_n=n(n+1)/2\), což pro počítač znamená konstantní počet elementárních operací (tři). To zapisujeme jako \(O(1)\) (konstantní časová složitost) a je to malý zázrak, když takový algoritmus potkáme. Ale stalo se a my místo doslovného programování sumace nejprve sáhneme po součtovém vzorci.

Ukázka č. 1: vzorec

Skript se součtovým vzorcem může být takto stručný:

# Python: vzorec
n=int(input('Zadejte n: '))  # vstup dat z klávesnice, konverze řetezce na celé číslo
# n=1000                     # vstup dat přiřazením v zdrojovém souboru
s=n*(n+1)/2
print('n,s =',n,s)

K přiřazovacímu příkazu a příkazu pro výpis zde přibyla výzva ke vstupu dat z klávesnice. Na výstupu funkce input je řetězec, pro použití v součtovém vzorci jej převedeme konverzní funkcí int na celé číslo. To celé uloženo jako skript spustíme buď příkazem python a.py nebo ve VS Code šipkou vpravo nahoře pro Run Python Code nebo Run Code (Ctrl+Alt+N).

Na obrazovce vidíme 1000 500500.0, tedy vstup jako celé a výsledek jako reálné číslo. To tam vniklo pomocí operátoru / pro reálné dělení, a – jak jsme viděli minule – povede to při větších hodnotách n k nepřesnému výsledku (typ float má jen 15 platných desítkových číslic). Přesvědčme se o tom vložením vstupní hodnoty n=1_000_000_001 (vidíme schopnost Pythonu zapisovat velká čísla čitelně pomocí podtržítek) a porovnejme s výpočtem, ve kterém užitím operátoru // pro celočíselné dělení udržíme výsledek v celočíselném typu (v Pythonu neomezeného rozsahu):

# Python: vzorec s typy float a int
# n=int(input('Zadejte n: '))        # vstup dat z klávesnice
n=1_000_000_001
s=n*(n+1)/2                          # reálné dělení
print('n,s =',n,s,'%.1f'%s,type(s))  # type(s) float
s=n*(n+1)//2                         # celočíselné dělení
print('n,s =',n,s,'%.1f'%s,type(s))  # type(s) int

Výpis

1000000001 5.000000015e+17 500000001500000000.0 <class 'float'>
1000000001 500000001500000001 500000001500000000.0 <class 'int'>

dokládá, že s operátorem / jsme získali reálné s, které po vyžádání všech číslic před desetinnou tečkou postrádá jedničku, zatímco s operátorem // jsme zůstali u přesného celočíselného s. Ovšem, pokud toto celočíselné s postoupíme formátové specifikaci pro výpis reálných čísel '%.1f'%s, Python hodnotu zkonvertuje do typu float a je po jedničce.

To samé s užitím NumPy vypadá, že bude počítat stejně, ale nebude:

# Python & NumPy: vzorec s typy np.float64 a np.int64
import numpy as np
# n=np.int64(input('Zadejte n: '))
n=np.int64(1_000_000_001)
s=n*(n+1)/2                          # reálné dělení
print('n,s =',n,s,'%.1f'%s,type(s))  # type(s) np.float64
s=n*(n+1)//2                         # celočíselné dělení
print('n,s =',n,s,'%.1f'%s,type(s))  # type(s) np.int64

Zvětšíme-li hodnotu n na 10 miliard, její druhá mocnina v součtovém vzorci přeteče horní mez typu np.int64 a máme tento pozoruhodný výsledek:

10000000001 3.883139830726121e+18 3883139830726120960.0 <class 'numpy.float64'>
10000000001 3883139830726120961 3883139830726120960.0 <class 'numpy.int64'>

Zpátky k realitě bychom se v rámci NumPy mohli vrátit použitím reálného n, např. příkazem n=np.float64(10_000_000_001).

Pro srovnání, fortranský ekvivalent má reakcí nejblíže ke skriptu Python & NumPy. Ovšem – operátor / se ve Fortranu (i C) rozhoduje pro celočíselné nebo reálné dělení podle typu operandů, a typ cílové proměnné v přiřazovacích příkazech není dán dynamicky typem výrazu na pravé straně (jako v Pythonu), ale statickými deklaracemi zkraje programu. Také podtržítka v literálech mají jiný význam než v Pythonu:

! Fortran: vzorec s typy real(8) a integer(8)
real(8) s       ! deklarace proměnné typu real(8)
integer(8) n,si ! deklarace proměnných typu integer(8)
! print '(a$)','Zadejte n: '; read *,n  ! vstup dat z klávesnice
n=1000000001_8  ! 8bytový literál
s=n*(n+1)/2     ! celočíselné dělení, přiřazení s implicitní konverzí do real(8)
print *,'n,s =',n,s; print '(f0.1)',s
si=n*(n+1)/2    ! celočíselné dělení, výsledek typu integer(8)
print *,'n,s =',n,si; print '(f0.1)',real(si,kind=8)
end program

Připomínáme znalost od minula, že k překladu fortranského zdrojového textu do proveditelného souboru by měl stačit příkaz pro volání překladače ve tvaru gfortran a.f90 z příkazového řádku Windows nebo terminálu VS Code nebo ve VS Code pomocí extenze Code Runner (Ctrl+Alt+N).

Ukázka č. 2: spaghetti code

Nyní se pokusíme náš algoritmus (zde s otočenou podmínkou, špagety budou propletenější) zapsat ve Fortranu prakticky doslovně:

  • zjisti parametr \(n\),

  • inicializuj nulami index \(i\) a proměnnou \(s\) pro kumulaci součtu \(S_n\),

  • otestuj, zda \(i\ge n\),

  • pokud ano, pokračuj na závěrečný výpis,

  • pokud ne, zvětši \(i\) o 1, přičti aktuální hodnotu \(i\) k mezisoučtu \(s\) a vrať se na předchozí test,

  • vypiš slovem, co vypisuješ, pak hodnotu parametru \(n\) a výsledek \(s\).

! Fortran: spaghetti code
real(8) s              ! real(8) promenna
integer(8) n,i         ! integer(8) promenne

! print '(a$)','Zadejte n: '; read *,n  ! vstup dat z klávesnice
n=1000000001_8         ! 8bytovy literal
i=0; s=0               ! inicializace
1 if (i>=n) goto 2     ! podmineny prikaz, prikaz skoku, navesti
i=i+1; s=s+i           ! update promennych i a s
goto 1                 ! prikaz skoku na radek s navestim 1
2 print *,'n,s =',n,s  ! vypis
end program

Kámen úrazu jsou příkazy skoku goto. Zde jsou sice použity k přesné realizaci rozumného algoritmu, ale v praxi by příkazy skoku mohly mířit celkem kamkoliv, úmyslně i omylem. A odchytit případně mylné příkazy skoku ve větším programu, to je jako probírat na talíři špagety carbonara. Jazyky tedy příkazy skoku ještě tolerují (Python však nikoliv), ale pedagogové svým studentům používání goto promíjet nemohou. Komu by však příkaz goto přišel jako čajíček, může si diskrétně rozšířit svůj arzenál o příkaz comefrom (wikipedia, originální článek).

V Pythonu skoky nejsou a špagety nebudou. (Pravda, 1. dubna 2004 byl navržen modul goto.) Ostatně takto programovat nebudeme ani ve Fortranu, C a nikde jinde.

Ukázka č. 3: cykly

A nyní strukturované programování. Viděli jsme, jak vyjádřit cyklus neboli opakování části programu pomocí explicitního příkazu skoku. Nyní místo toho použijeme příkazové konstrukce pro cykly, které v sobě skok sice také obsahují, ale implicitně a bezpečně, právě jen ten skok, který je potřeba, bez možnosti zmýlení. Příkazové konstrukce mají jasně dán svůj vstupní a výstupní bod (na začátku a konci); to s příkazy goto nebylo možné zajistit.

Cyklus while s prováděcí podmínkou realizuje náš algoritmus prakticky doslovně, přičemž zmínky o skocích můžeme učinit implicitními i ve slovním popisu:

  • zjisti parametr \(n\),

  • inicializuj nulami index \(i\) a proměnnou \(s\) pro kumulaci součtu \(S_n\),

  • dokud \(i<n\), zvětšuj \(i\) o 1 a přičítej aktuální hodnotu \(i\) k mezisoučtu \(s\),

  • vypiš slovem, co vypisuješ, pak hodnotu parametru \(n\) a výsledek \(s\).

Dokud je podmínka splněna, blok (tělo) cyklu se opakovaně provádí, je-li při testu podmínka nesplněna, program pokračuje (implicitním skokem) za blok cyklu. V Pythonu se blok cyklu vyznačuje odsazením.

# Python: cyklus s podmínkou (cyklus while)
# n=int(input('Zadejte n: '))  # vstup dat z klávesnice
n=100_000_001
i=0; s=0.
while i<n:                     # cyklus s podmínkou
  i+=1; s+=i                   # blok cyklu: update proměnných i a s
print('n,s =',n,s)

V ukázce máme přiřazovací příkazy s běžnou zkratkou pro update: místo i=i+1 lze psát i+=1, místo s=s+i lze s+=i a existují další varianty. Můžeme poukázat i na relační operátory, které jsme dnes začali používat – je jich šest a zapisují se v Pythonu takto: ==, != (Fortran: /=), <, <=, > a >=. Platí zde to samé, co jsme probrali u první ukázky, totiž že přiřazením s=0. budeme sčítat v reálné aritmetice, tedy při větších n nepřesně, zatímco s=0 zavede celočíselné s a sčítat se bude v neomezeném typu int. Ovšem – rozdílu bychom se měli dočkat až pro n=1_000_000_001, což v Pythonu bude trvat trochu dlouho, přes minutu.

Fortran to bude mít s miliardou hotové v sekundách:

! Fortran: cyklus s podmínkou
real(8) s
integer(8) n,i
! print '(a$)','Zadejte n: '; read *,n  ! vstup dat z klávesnice
n=1000000001_8
i=0; s=0
do while (i<n)  ! cyklus s podminkou
  i=i+1; s=s+i  ! blok cyklu: update proměnných i a s
end do
print *,'n,s =',n,s
end program

Takto zapsáno už je to čisté strukturované programování. Přesto lze ještě lépe. Cyklus o přesném počtu opakování, řízený indexem, je v praxi natolik frekventovaný, že si vysloužil vlastní konstrukci, zvanou indexovaný cyklus neboli cyklus s řídicí proměnnou. Indexovaný cyklus implicitně zajišťuje růst/klesání indexu od počáteční meze ke koncové mezi a testování, zda už index dosáhl koncové meze. Krok bývá implicitně +1, většinou lze v jazycích nastavit i jinak (neměl by být nulový). V těle cyklu zde tedy zůstává explicitně jen update proměnné s. Přizpůsobíme slovní popis algoritmu:

  • zjisti parametr \(n\),

  • inicializuj proměnnou \(s\) pro kumulaci součtu \(S_n\),

  • pro všechny hodnoty indexu \(i\) od 1 do \(n\) přičítej aktuální hodnotu \(i\) k mezisoučtu \(s\),

  • vypiš slovem, co vypisuješ, pak hodnotu parametru \(n\) a výsledek \(s\).

Python pro indexovaný cyklus vytěžuje svůj vcelku obecný cyklus for s funkcí range(start,stop,step); argumenty start a step jsou nepovinné s default hodnotami po řadě 0 a 1, index projde hodnoty od start s krokem step až před stop, tedy bez koncové hodnoty stop.

# Python: indexovaný cyklus (cyklus for)
# n=int(input('Zadejte n: '))
n=100_000_001
s=0.
for i in range(1,n+1):  # indexovaný cyklus od 1 do n (bez koncove hodnoty)
  s+=i   # blok cyklu: update proměnné s explicitně, update indexu i implicitně
print('n,s =',n,s)

V našich ukázkách je pythonský indexovaný cyklus bezmála dvakrát rychlejší než cyklus s podmínkou.

Fortran má indexovaný cyklus v podobě konstrukce do index=start,stop[,step]end do, kde hodnota stop je také zahrnuta.

! Fortran: indexovaný cyklus
real(8) s
integer(8) n,i
! print '(a$)','Zadejte n: '; read *,n
n=1000000001_8
s=0
do i=1,n  ! indexovany cyklus od 1 do n vcetne, krok implicitne 1
  s=s+i   ! blok cyklu: update promenne s explicitne, update indexu i implicitne
end do
print *,'n,s =',n,s
end program

Není rozdílu v rychlosti běhu fortranských ukázek s příkazy goto, cyklem do while a cyklem for.

Benchmark

Máme konečně v ruce něco, co bude počítat tak dlouho, že už můžeme měřit na hodinách dobu běhu (wall-clock time), tedy čas, který výpočtáře zajímá nejvíc. Ze soutěže diskvalifikujeme příliš rychlou ukázku č. 1 a poměříme si (provedeme benchmark) řešení ukázek č. 2 a 3 ve Fortranu, C, Pythonu a dalších jazycích. Překladače mívají různé volby, kterými lze zvýšit numerický výkon exe-souboru (a obvykle snížit kvalitu diagnostiky případných chyb). Úroveň této optimalizace nastavujeme např. v GNU Fortranu volbou gfortran -O a.f90 (default -O je druhá hladina, -Ofast je sada nejúčinnějších nastavení).

Podobný benchmark provedeme i pro výpočet harmonického čísla \(H_{1000000000}\). Porovnáme tím náročnost operace dělení proti všem ostatním operacím potřebným pro chod cyklu, které jsme změřili výpočtem čísla \(S_{1000000000}\). Ve Fortranu je třeba kumulativní příkaz zapsat výslovně v real(8) přesnosti: s=s+1._8/i, v Pythonu postačí s=s+1/i.

Tabulka dob běhu ukázek pro výpočet \(S_n\) (typ procesoru: Intel i9-9900K, Linux, zdrojové soubory pro S_n)

jazyk

n

ukázka 2

ukázka 3a

ukázka 3b

gfortran

1000000000

1.8

1.8

1.8

gfortran -O

1000000000

0.8

0.8

0.8

gcc

1000000000

1.8

1.8

1.8

gcc -O

1000000000

0.8

0.8

0.8

Python

100000000

8.8

4.9

Matlab

1000000000

0.8

0.8

Octave

10000000

14.1

5.8

Časy překládaných programů (Fortran, C) by pro n rovno miliardě měly být v řádu sekund, což dovoluje k měření použít jen hodinky; softwarové možnosti zmiňujeme níže. Vidíme, že bez optimalizace je běh programů přeložených Fortranem a C stejně rychlý. Optimalizace zde zrychluje na dvojnásobek. Není rozdíl mezi rychlostí špagetového kódu a strukturovaných cyklů. Optimalizovaného překladu cyklů je zjevně schopen i Matlab. Smutný pohled je na rychlost cyklů v Octavu, o tři řády níže než v Matlabu. I pythonské cykly jsou téměř o dva řády pomalejší než cykly přeložené.

Tabulka dob běhu ukázek pro výpočet \(H_n\) (typ procesoru: Intel i9-9900K, Linux, zdrojové soubory pro H_n)

jazyk

n

ukázka 2

ukázka 3a

ukázka 3b

gfortran

1000000000

4.4

4.4

4.4

gfortran -O

1000000000

0.9

0.9

0.9

gcc

1000000000

4.4

4.4

4.4

gcc -O

1000000000

0.9

0.9

0.9

Python

100000000

9.1

5.5

Matlab

1000000000

0.9

0.8

Octave

10000000

15.8

7.4

Dělení se ukazuje jako násobně pomalejší operace než všechny ostatní operace cyklu dohromady. Optimalizace překladů to ovšem zahladí. V interpretovaných (za chodu překládaných) skriptech Pythonu a Octavu se časová náročnost dělení také ztrácí.

Pro softwarové měření doby běhu programů lze použít vhodný timer. V Linuxu je to běžná součást, na program a se nasadí příkazem time a. Do Windows je třeba si něco sehnat – já zrovna teď mám timer odtud, který používám k měření doby běhu programu a v sekvenci tří příkazů timer /q && a && timer /s /q. Ve vývojovém prostředí Matlabu a Octavu lze měřit pomocí dvojice příkazů tic; a; toc. V inspiraci tímto názvoslovím můžeme v Pythonu použít modul time:

# Python: měření doby běhu pomocí modulu time
import time
tic=time.time()  # start
time.sleep(2)    # dvousekundový mikrospánek
toc=time.time()  # stop
print('sleep duration:',toc-tic)

Všechny ukázkové programy jsou sbaleny zde: ukazky-02-sumace.zip.