Programování pro fyziky 2023/24

Úvod a cíl

Toto je první z průvodních textů k přednášce Programování pro fyziky. Jeho části jsou nabídnuty v interaktivním prostředí Google Colab (notebook 01), kde je také shrnuto několik otázek. Odpovědi na ně lze vkládat v prostředí Moodle. Níže použité skripty jsou sbaleny v tomto zipu. Pro další praktické zkušenosti lze zajít na stránky cvičení (přednostně ovšem na stránky svého cvičícího).

Cílem přednášky je poskytnout povědomí o postupech stylů programování nazývaných

  • imperativní (program s podrobným zápisem příkazů realizujících algoritmy),

  • strukturované (omezený repertoár příkazů a příkazových konstrukcí pro lepší srozumitelnost programu),

  • procedurální (členění programu do procedur).

Postupně více a více nám bude ulevovat deklarativní styl (naznačíme, co chceme, nepopisujeme podrobně, jak to udělat). Stylům z vyšší ligy programování, jako objektově orientované programování, funkcionální programování nebo paralelní programování, se budeme věnovat jen trochu, to je úkolem jiných předmětů.

Pozornost zaměříme také k numerickým metodám používaným v našem oboru, konkrétně

  • řešení soustav lineárních algebraických rovnic,

  • polynomická aproximace,

  • numerické integrování,

  • hledání kořenů,

  • řešení soustav obyčejných diferenciálních rovnic.

Programy se zapisují v programovacích jazycích. Pro fyzika jsou relevantní především

  • C/C++ a Fortran, řekněme klasické programovací jazyky, díky realizaci pomocí překladačů produkující rychlé programy a nabízející rozvinutou diagnostiku,

  • Python, Matlab/Octave, R, Julia ad., řekněme skriptovací jazyky, obvykle interpretované, tj. překládané po částech během provádění programu, proto běžící pomaleji, volené často spíše pro menší programy (skripty).

Tento výčet lze konfrontovat např. s webem Národního superpočítačového centra IT4Innovations v Ostravě: tamní software.

Fortran

Zahájíme přípravou našeho nultého programu, a to hned v několika z uvedených jazyků. Začneme tím, čím to všechno začalo a dosud neskončilo – s Fortranem (také: fortran-lang.org). Vznikl jako první programovací jazyk světa roku 1957, v roce 1962 už prý dokázal srazit družici, prochází sedmou oficiální (ANSI/ISO) revizí, na IT4Innovations je uveden hned na druhé pozici, jeho překladače vyrábějí jak společenství GNU, tak firmy jako Intel a Nvidia. Pro následující akci nám stačí volně dostupný překladač z GNU Compiler Collection (GCC), který si pro Windows můžete obstarat např. equation.com. (Instalujme do C:\gcc, jak je doporučeno v závěrečných poznámkách zde dole.)

Dále pokračovat lze různými cestami. Tou nejobyčejnější je jednoduchá editace a překlad z příkazového interpretu. Command Prompt nalezneme buď v nabídce Windows System nebo pomocí Search řádku nebo jej spustíme v okně Run (klávesová kombinace Win+R) zadáním cmd. Vhodný je přechod do pracovního adresáře, pro který můžeme zvolit pracovní plochu: cd Desktop (change directory). Otevřeme editor pro editaci souboru a.f90 (přípony často indikují, v jakém programovacím jazyce je program psán; starodávnou příponu for nahradila po reformě Fortranu v devadesátých letech přípona f90): notepad a.f90. Napíšeme program:

print *,'ahoj'
print *,1,1.
end

Soubor uložíme a anglicky psaný zdrojový soubor (source file) programu přeložíme do spustitelného (executable) souboru a.exe příkazem gfortran a.f90. Instrukcím v něm zapsaným už bude rozumět procesor počítače, takže exe-soubor můžeme spustit: a. (Nemusíme psát a.exe.)

Efektivnější a uspokojivější než Notepad jsou programátorské editory, mezi nimiž na špičce stojí Visual Studio Code. Ten budeme používat i my na přednášce. Postup pro jeho konfiguraci pro námi užívané jazyky shrneme v jiném textu.

Tedy ať už v Command Prompt, terminálu VS Code nebo v jiném prostředí uvidíme po spuštění našeho programu toto:

ahoj
          1   1.00000000

Příkaz výstupu print zajistil výpis literálu (doslovné hodnoty) 'ahoj', říkat literálům v apostrofech můžeme také řetězec (znaků), a výpis dvou jedniček. To už nás zajímá hodně: procesor zjevně rozumí jednak celým číslům (datový typ integer), jednak číslům s desetinnou tečkou (datový typ: floating point alias real). Už dlouhá léta nás obklopují 32bitové a 64bitové procesory a operační systémy a to se projevuje i tím, kolik prostoru v paměti je pro tyto jedničky vyčleněno. V případě integer dat to jsou nejčastěji 4 byty (32 bitů), což dovoluje vyjádřit 2^32 různých hodnot. Obvyklé je pokrýt hodnoty kladné i záporné, takže 4bytová celočíselná znaménková aritmetika umí pracovat s hodnotami v rozsahu -2^31 až +2^31-1, což jsou přibližně ±2 miliardy. Fortran nijak nezakrývá, co se stane, když vybočíme neboli přetečeme mimo tento rozsah (tzv. integer overflow):

print *,'ahoj'
print *,1,1.
MAX=huge(1)  ! funkce huge pro nejvetsi reprezentovatelne cislo datoveho typu daneho vyrazem v argumentu
print *,MAX,MAX+1  ! vypis promenne a vyrazu, v kterem k promenne MAX pricitame 1
end

Na výstupu přibyl řádek 2147483647 -2147483648 a vidíme, že vpravo za největším reprezentovatelným celým číslem leží nejmenší reprezentovatelné celé číslo; celá čísla tedy v sobě dostupném rozsahu rotují, zobrazené celé číslo je výsledkem operace modulo. Ta dlouhá mezera, kterou GNU Fortran vynechal před 1, je teď téměř zaplněná, tedy nám tím už prve naznačil, že je schopen zapsat celé číslo až o 10 cifrách.

Nyní reálná 1. Ve Fortranu je rovněž 4bytová. Vnitřní uspořádání oněch 4 bytů v reálném datovém typu probereme později, nyní se seznámíme s výsledným projevem v desítkové soustavě. Výpis naznačuje, že ona jednička není přesnější než jen na 7-8 desetinných míst. To je pro praxi málo, odedávna proto hardware i software nabízí i dvojnásobnou, 8bytovou reálnou přesnost. Ve Fortranu se na ni můžeme podívat takto:

print *,'ahoj'
print *,1,1.
MAX=huge(1)                  ! funkce huge pro nejvetsi reprezentovatelny integer
print *,MAX,MAX+1
print *,1./3.,1._8/3._8      ! za podtrzitkem velikost pameti vyhrazene pro literal v bytech (mozne hodnoty: 4, 8, nekdy 16)
print *,huge(1.),huge(1._8)  ! funkce huge pro nejvetsi reprezentovatelny 4bytovy a 8bytovy real
end

Na výstupu přibylo:

0.333333343      0.33333333333333331
 3.40282347E+38   1.7976931348623157E+308

tedy pokyn spočítat třetinu ve 4bytové, resp. 8bytové přesnosti naznačil to, co později doložíme, že výsledek bude matematicky správný jen na prvních 7, resp. 15 platných míst. Poslední řádek navíc sděluje hodnotu desítkového exponentu, ke kterému mohou došplhat reálná čísla v obou přesnostech. Výstupy takto velkých čísel se dějí v semilogaritmickém tvaru, složeném z mantisy v rozsahu <1,10) a desítkového exponentu za písmenem E.

A konečně, zkusíme vypsat číslo větší než huge s reálným argumentem (real overflow). Budeme opět potřebovat proměnné (výše jsme užili proměnnou MAX), které je tentokrát třeba explicitně popsat (deklarovat) na začátku programu: real(4) MAX4; real(8) MAX8. Správná odpověď má velmi univerzální platnost, produkuje ji totiž procesor (hardware) a programovací jazyky (software) ji sice mohou poněkud upravit, ale vlastně ji nijak nezakrývají. (Na rozdíl od integer overflow, jak dále uvidíme.) Cílová verze fortranské ukázky

! Fortran: jedničky a nekonečna
program a
real(4) MAX4; real(8) MAX8
print *,'ahoj'
print *,1,1.,1._8
MAX=huge(1)              ! funkce huge pro nejvetsi reprezentovatelne cislo datoveho typu daneho vyrazem v argumentu
MAX4=huge(1.)
MAX8=huge(1._8)
print *,MAX,MAX+1        ! vypis promenne a vyrazu, v kterem k promenne MAX pricitame 1
print *,1./3.,1._8/3._8  ! za podtrzitkem velikost pameti vyhrazene pro literal v bytech (mozne hodnoty: 4, 8, nekdy 16)
print *,MAX4,MAX8
print *,MAX4*2,MAX8*2
end

vypíše:

ahoj
          1   1.00000000       1.0000000000000000
  2147483647 -2147483648
  0.333333343      0.33333333333333331
  3.40282347E+38   1.7976931348623157E+308
        Infinity                  Infinity

Počítačová reálná čísla mají svá nekonečna. Pro 4bytový reálný typ, často nazývaný single precision, přichází po 10^38, pro 8bytovou double precision po 10^308. Počítačové nekonečno může mít znaménko a utéct z něj zpět do reprezentovatelného rozsahu lze jen pomocí dělení. Nikoliv snad jako Infinity/2, ale 1/Infinity už nula bude. Na později si necháme téma 0./0., 0.*Infinity a Infinity/Infinity.

C

Uvedeme ekvivalentní program i v jazyce C:

// C: jedničky a nekonečna
#include <stdio.h>
#include <limits.h>
#include <float.h>
int main() {
  int MAX; float MAX4; double MAX8;
  printf("%s\n","ahoj");
  printf("%d %f %.17f\n",1,1.f,1.);
  MAX=INT_MAX; MAX4=FLT_MAX; MAX8=DBL_MAX;
  printf("%d %d\n",MAX,MAX+1);
  printf("%.9f %.17f\n",1.f/3.f,1./3.);
  printf("%.8e %.16e\n",MAX4,MAX8);
  printf("%e %e\n",MAX4*2,MAX8*2);
  return 0;
}

Rozdíly jsou spíše drobné, ale různorodé:

  • příkazy a popisy v C se oddělují středníky,

  • některé z použitých jmen nejsou k dispozici v holém C, proto bývá třeba poskytnout knihovny nebo jejich části (zde stdio, limits, float),

  • proměnné (zde: MAX ad.) je v C nutné popisovat explicitně; ve Fortranu platí jistá implicitní pravidla, ale explicitní popisy jsou vhodnější i tam,

  • funkce pro formátovaný výpis printf vyžaduje explicitní formátovou specifikaci, zde %s, %d, %f a %e pro výpis řetězců (strings), celých čísel v desítkové soustavě (decimal) a reálných čísel (float), včetně semilogaritmického tvaru (obsahujícího desítkový exponent), a je třeba i explicitně odřádkovat pomocí \n (newline),

  • Pascal stejně jako Fortran zpřístupňuje 4- i 8bytový reálný typ, ale default (implicitní hodnota) je 8 bytů. O 4bytové protějšky 8bytových literálů si v C musíme říci příponou literálu f.

Překládat můžeme překladačem gcc z balíčku GCC; příkazem gcc -o b b.c se output překladače nazve b.exe. Na jeho výstupu pak bude:

ahoj
1 1.000000 1.00000000000000000
2147483647 -2147483648
0.333333343 0.33333333333333331
3.40282347e+38 1.7976931348623157e+308
inf inf

což je velmi blízké výpisu z fortranského ekvivalentu.

Python

Do třetice použijeme skriptovací jazyk Python (home: python.org). Skriptem je zvykem nazývat program pro skriptovací jazyk, obvykle je to spíše menší program. Python a skriptovací jazyk obecně není překladač, ale interpret – předložený program se nepřekládá před spuštěním, ale přímo vykonává/interpretuje. (Tedy skript se překládá po menších částech během provádění.) Python si lze opatřit z jeho domácí stránky, pohodlné je to i pomocí Microsoft Store z Windows. Možné je získat Python jako distribuci, ve které je k němu přibalena i řada populárních balíčků rozšiřujících jeho schopnosti. Distribucí oblíbenou v akademických kruzích je Anaconda, případně její menší sestra Miniconda. Snadno přidávat balíčky lze i k holému Pythonu z webu; minimální sadu vhodnou pro naše potřeby zajistí příkaz: pip install numpy scipy matplotlib numba sympy jupyter.

Po instalaci můžeme buď po spuštění příkazu python v příkazovém interpretu zadávat následující příkazy jeden po druhém (říká se tomu REPL, read-eval-print loop), nebo je uložit jako skript c.py a vše spustit příkazem python c.py:

# Python bez NumPy: jedničky a nekonečna
import sys
print('ahoj')       # řetězec typu str ('string')
print(1,1.)         # jedničky typu int (`integer`) a float ('floating-point numbers')
MAX=sys.maxsize     # 2**31-1 nebo 2**63-1
print(MAX,MAX+1,2**128,float(2**128))  # integer overflow?
print(1/3)          # kolik platných číslic?
MAX8=sys.float_info.max  # největší 64bitový float
print(MAX8,MAX8*2)  # real overflow?

Podobně jako v C, k holému Pythonu bývá potřeba připojovat moduly, zde sys. Python stejně jako jazyky v ukázkách výše zpřístupňuje celočíselné i reálné počítání. Povědomí Pythonu o 8bytových celých číslech naznačuje hodnota sys.maxsize čili 2**63-1, ovšem celočíselná aritmetika Pythonu na hardwarových 8 bytech ani jinde nebrzdí a můžeme se klidně podívat třeba na 2**128 nebo na 1000 nul umocněním celočíselné desítky, 10**1000 – vidíme ukázku neomezené celočíselné aritmetiky, emulované v Pythonu softwarově. O 4bytovou reálnou aritmetiku se holý Python nezajímá, podpoří jen 8bytovou; reálné nekonečno přiznává. Výstup skriptu:

ahoj
1 1.0
9223372036854775807 9223372036854775808 340282366920938463463374607431768211456 3.402823669209385e+38
0.3333333333333333
1.7976931348623157e+308 inf

Ovšem, bez hardwarového int a 4bytového float by to při počítání nebylo ono. To a samozřejmě mnohem více nabízí balíček pro numerické počítání v Pythonu NumPy, obvykle do Pythonu importovaný pod zkratkou np příkazem import numpy as np. O celočíselných typech informuje struktura np.iinfo, o floating-point typech struktura np.finfo. V NumPy se nalezne více variant celočíselných i reálných typů, můžeme tedy zkoumat, jak se projevuje přetečení v 32- i 64bitové přesnosti. Můžeme se i podívat na přesnost 1/3 počítané v 32- a 64bitové aritmetice. Přetečení jsou svého druhu chybový stav, i když ne příliš kritický, tak pro jejich pozorování bude vhodné vypnout pythonská varování, import warnings; warnings.filterwarnings('ignore').

# Python s NumPy: jedničky a nekonečna
import numpy as np
import warnings
warnings.filterwarnings('ignore')     # potlačíme zbytná varování
print('ahoj')
print(1,1.)
MAX=np.iinfo(np.int32).max            # integer-info o 32bitové variantě
MAX=np.int32(MAX)
print(MAX,MAX+1)                      # integer overflow? někdy ano, někdy ne
MAX=np.iinfo(np.int64).max            # integer-info o 64bitové variantě
MAX=np.int64(MAX)
print(MAX,MAX+1)                      # integer overflow? ano
MAX4=np.finfo(np.single).max          # float-info o 32bitové variantě
print(MAX4,MAX4+MAX4,MAX4*2)          # real overflow? ano, ale...
MAX8=np.finfo(np.double).max          # float-info o 64bitové variantě
print(MAX8,MAX8+MAX8,MAX8*2)          # real overflow? ano
print(np.single(1)/np.single(3),1/3)  # kolik platných číslic?
print("%.33f"%(np.single(1)/np.single(3)),"%.66f"%(1/3))  # kolik nenulových číslic?

Na výpisu tohoto skriptu, tak podobného programům ve Fortranu a C, nás může – překvapivě – leccos překvapit. Tak například obsahuje-li proměnná MAX horní mez 32bitového typu np.int32, výraz MAX+1 ve Windows přeteče; v Linuxu (a tedy i v Google Colabu) nepřeteče. Pro MAX typu np.int64 výraz MAX+1 přetéká tam i tam. Reálným nekonečnům NumPy rozumí v 32bitovém typu np.single (neboli np.float32) i 64bitovém np.double (np.float64). Ovšem do nekonečna nás nedovede (na rozdíl od Fortranu a C) výraz MAX4*2, ale MAX4+MAX4. Ostatně, mohli bychom se ptát, proč vlastně nestačí MAX4+1 (to ani ve Fortranu a C).

ahoj
1 1.0
2147483647 -2147483648
9223372036854775807 -9223372036854775808
3.4028235e+38 inf 6.805646932770577e+38
1.7976931348623157e+308 inf inf
0.33333334 0.3333333333333333
0.333333343267440795898437500000000 0.333333333333333314829616256247390992939472198486328125000000000000

Nápovědu může přinést pythonská funkce type(), sdělující datový typ argumentu (např.: print(type(1),type(1.))). Něco naznačit může i poslední řádek výpisu, kde na formátovaných výpisech floating-point čísel, ovšemže v desítkové soustavě, vidíme několik platných číslic, pak několik nesprávných nenulových číslic a pak nuly. To vše bude srozumitelnější, až se propracujeme k vyjádření floating-point čísel v dvojkové soustavě.

Konkrétněji otázky k těmto dvěma pythonským skriptům shrnuje shora odkázaný notebook v Google Colabu. Odpovědi systému ChatGPT na ně jsou nabídnuty ve shora odkázaném prostředí Moodle, kam lze vložit i vlastní názory na to, co to tu vlastně pozorujeme.

Matlab a Octave

Systém pro numerickou analýzu Matlab (home: mathworks.com) vyniká úsporností své syntaxe. V lineární algebře (při počítání mj. s maticemi, jak slibuje jeho název) posloužil jako vzor pro svět Pythonu. Na skripty aplikuje transparentně překladač, takže nebývají pomalé. Je komerční, ale k dispozici je univerzitní licence. Inspiroval vznik volně dostupného klonu GNU Octave (home: octave.org). Pro účely naší přednášky jsou oba systémy zástupné; instalace Octavu je časově i prostorově podstatně skromnější.

Skript d.m

% Matlab/Octave: jedničky a nekonečna
format compact  % mene prazdnych radku na vypisu
format long     % vice desetinnych mist na vypisu
'ahoj'
1,1.
intmax,2^31,2^63,2^128
1/3
realmax,realmax*2

se spustí v IDE (integrated development environment) obou systémů buď po řádcích interaktivně nebo uvedením jména d, nebo zvnějšku voláním matlab -batch d, resp. octave d.m (ve Windows při existenci spustitelného skriptu octave.bat s obsahem @C:\Octave\octave-launch --no-gui %1). Z výstupu Octavu

ans = ahoj
ans = 1
ans = 1
ans = 2147483647
ans = 2147483648
ans = 9.223372036854776e+18
ans = 3.402823669209385e+38
ans = 0.333333333333333
ans = 1.797693134862316e+308
ans = Inf

je cítit, že celá čísla jsou v defaultu reprezentována implicitním 8bytovým reálným typem; v oblasti používání obou systémů to nečiní zvláštní problémy.

Gnuplot

Gnuplot (home: gnuplot.info) budeme vedle Pythonu používat pro grafickou vizualizaci našich numerických výsledků, především pro stručnost jeho syntaxe, mocnou matematickou výbavu a svižnost při skriptování. Je praktický i jako kalkulačka, pro nás teď i proto, že i on rozlišuje integer a real datové typy. V následujícím skriptu nás může zaujmout výsledek snahy pomáhat s integer přetečením:

# Gnuplot: jedničky a nekonečna
print 'ahoj'
print 1,1.
pr 2**30-1+2**30,2**30-1+2**30+1,2**31,2**63
pr 1/3,1./3.
pr 1e300,1e400

K horní mezi 4bytového integer se blížíme zleva, opatrně (+1) ji překročíme, a překročíme ji také zhurta (2**31). Kombinovaný výstup gnuplotu verze 4 a soudobé verze 5 ukazuje, že svého času se gnuplot snažil bránit integer overflow přechodem od 4bytového integeru k reálnému typu, zatímco nyní obranu posílil povoláním 8bytových integerů:

ahoj
1 1.0
2147483647 -2147483648 2147483648.0 9.22337203685478e+018  # verze 4.6.7
2147483647 2147483648 2147483648 9.22337203685478e+18      # verze 5.4.2
0 0.333333333333333
1e+300 inf.0

Legračního je toho na výstupu více, mj. výsledek podílu 1/3 (to udělá i Fortran a C – brzy si objasníme, proč) a nekonečno na jedno desetinné místo. Pokud se vám podařilo nainstalovat gnuplot a spustit ho v adresáři Desktop (polohu napoví příkaz gnuplotu pwd, print working directory, a je-li čerstvý gnuplot v adresáři Documents, k přesunu poslouží příkaz cd '..\Desktop'), můžete uvedené příkazy uložit do skriptu e.gp a v gnuplotu spustit příkazem load 'e.gp'.

Shrnutí

Shrneme-li dosavadní pozorování, vidíme, že procesory počítačů podporují celočíselné a (tzv.) reálné datové typy o několika přesnostech. Programovací jazyky zpřístupňují tuto nabídku a v situacích, kdy daná přesnost nestačí na potřebné matematické hodnoty, se mohou snažit pomoci, ovšem každý jinak. Je proto zřejmě optimální zvyknout si na hardwarová omezení datových typů a nevybočovat z nich.

Během přednášky jsme si také užili editoru Visual Studio Code s jeho schopností porozumět mnoha programovacím jazykům (tj. barvit syntaxi, vyznačovat chyby a leccos poradit) a pohodlným spouštěním programů z příkazového řádku (command-line interface, CLI) nebo pomocí extenze Code Runner. Všechny ukázkové programy jsou sbaleny zde: ukazky-01-overflow.zip.