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

Opakování

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

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

  • 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í.)

Příklady pro dnešek vyzobneme z poznámek k přednášce, str. 2-4, kde si o nich můžete přečíst. Tam jsou ve Fortranu, zde je převedeme do Pascalu a Pythonu, pokud se vám už zalíbil jiný z jazyků vyjmenovaných na předchozí přednášce, zkuste ten. Určitě si ukázky zkuste, v některém jazyce, s různým vstupním parametrem, s úpravami vstupu či výstupu dat.

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 první variantu (druhá je 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

Z poznámek k přednášce přejímáme součtový vzorec ve Fortranu:

! Fortran: vzorec
read *,n     ! prikaz vstupu
s=n*(n+1)/2  ! prirazovaci prikaz
print *,n,s  ! prikaz vystupu
end

Připomínáme znalost od minula, že k překladu tohoto 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). Pro n máme v paměti 4 byty, rozsah do zhruba 2 miliard. Fortran podrží datový typ operandů i pro výsledek, takže pozor na efektivně druhou mocninu n*(n+1) – ta zde přeteče už při \(n\ge46341\), takže doplňme ještě popis s přidělením 8bytového datového typu proměnné n:

! Fortran: vzorec s vyšší přesností
integer(8) n
! real(8) s
read *,n
s=n*(n+1)/2
print *,n,s
end

a vydržíme pak s n až k miliardám. Nepřesný výsledek pro miliardu potlačíme, když reálné proměnné s přidělíme také 8bytový prostor: real(8) s.

V Pascalu takový program se zkušeností od minula už dokážeme také. V Pascalu se navíc nemusíme starat o přetečení n až do miliard, neboť Pascal povyšuje automaticky vyčíslení do 8bytového typu Int64, a reálná proměnná s je implicitně 8bytová.

// Pascal: vzorec
program a;
var
  n : integer;
  s : real;
begin
  readln(n);
  s:=n*(n+1)/2;
  writeln(n,s);
end.

Překládáme buď příkazem fpc -Mobjfpc a.pas nebo pomocí Code Runneru nebo ve vývojovém prostředí Lazarus. Ať už zadáme jakkoliv velké n, běh programu stojí na 3 elementárních operacích a nebude trvat skoro nic.

Python bude stručný podobně jako Fortran:

# Python: vzorec
n=int(input('Zadejte n: '))
s=n*(n+1)/2
print(n,s)

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).

Ukázka č. 2: spaghetti code

Nyní se pokusíme náš algoritmus zapsat ve Fortranu prakticky doslovně. Ať to máme vedle sebe:

  • 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\).

! Fortran: spaghetti code
real(8) s
read *,n
i=0; s=0
111 if (i>=n) goto 999  ! řádek s návěštím, podmíněný příkaz, příkaz skoku
i=i+1; s=s+i
goto 111
999 print *,n,s
end

Jak vidíme, v programu jsme negovali podmínku algoritmu, zápis programu je pak poněkud přirozenější. (Následující ukázka č. 3 bude naopak přirozenější při zápisu ve shodě s algoritmem.) Problémem 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).

Ekvivalent v Pascalu:

// Pascal: spaghetti code
program b;
var
  n,i : integer;
  s : real;
label
  111,999;
begin
  readln(n);
  i:=0; s:=0;
111: if i>=n then goto 999;
  i:=i+1; s:=s+i;
  goto 111;
999: writeln(n,s);
end.

V Pythonu skoky nejsou a špagety nebudou. Ostatně takto programovat nebudeme ani ve Fortranu, Pascalu a nikde jinde.

Ukázka č. 3: cykly

A nyní strukturované programování. Program je zde sekvencí příkazů a příkazových konstrukcí, což jsou především a) rozvětvení a b) cykly, obojí s jasně daným vstupním a výstupním bodem. Fortran:

! Fortran: cyklus s podmínkou
real(8) s
read *,n
i=0; s=0
do while (i<n)  ! do-while cyklus
   i=i+1; s=s+i
end do
print *,n,s
end

Vidíme zde cyklus do while s prováděcí podmínkou: dokud je podmínka splněna, tělo (blok) cyklu se opakovaně provádí, je-li při testu podmínka nesplněna, program pokračuje za tělem cyklu. V Pascalu takto:

// Pascal: cyklus s podmínkou
program c1;
var
  n,i : integer;
  s : real;
begin
  readln(n);
  i:=0; s:=0;
  while i<n do begin  // while cyklus
    i:=i+1; s:=s+i;
  end;
  writeln(n,s);
end.

Python:

# Python: cyklus s podmínkou
n=int(input('Zadejte n: '))
i=0; s=0.
while i<n:  # while cyklus
  i=i+1; s=s+i
print(n,s)

Takto zapsáno už je to čisté strukturované programování. Ale 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:

! Fortran: indexovaný cyklus
real(8) s
read *,n
s=0
do i=1,n  ! indexovaný cyklus do
  s=s+i
end do
print *,n,s
end

Indexovaný cyklus tedy 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. V Pascalu se indexovaný cyklus uvozuje slůvkem for, což mu dává další jeho běžný název: for-cyklus.

// Pascal: indexovaný cyklus
program c2;
var
  n,i : integer;
  s : real;
begin
  readln(n);
  s:=0;
  for i:=1 to n do begin  // for cyklus
    s:=s+i;
  end;
  writeln(n,s);
end.

A konečně poslední verze poslední ukázky:

# Python: indexovaný cyklus
n=int(input('Zadejte n: '))
s=0.
for i in range(1,n+1):  # for cyklus s funkcí range
  s=s+i
print(n,s)

Python pro indexovaný cyklus vytěžuje svůj vcelku obecný for-cyklus s funkcí range(start,stop,step).

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, Pascalu, 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 v GNU Fortranu volbou gfortran -O a.f90 (default -O je druhá hladina, -Ofast je sada nejúčinnějších nastavení) a ve Free Pascalu podobně, fpc -Mobjfpc -O4 a.pas (pro naši ukázku je účinná až čtvrtá hladina).

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 Pascalu 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

fpc

1000000000

1.9

1.9

1.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, Pascal, C) by pro n rovno miliardě měly být v řádu sekund, což dovoluje k měření použít jen hodinky. Pro vyšší přesnost měření 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 v podobě tří příkazů timer && echo 1000000000 | a && timer /s, kde prostředním příkazem posílám na vstup programu a příkazem echo 1000000000 žádanou miliardu. V Matlabu a Octavu lze měřit ve vývojovém prostředí pomocí dvojice příkazů tic; a; toc.

Vidíme, že bez optimalizace je běh programů přeložených Fortranem, C a Pascalem stejně rychlý. Optimalizace zde zrychluje na dvojnásobek. (U Free Pascalu jen ve Windows, v Linuxu s mírně starší verzí překladače nikoliv.) 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

fpc

1000000000

1.8

1.8

1.8

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í.