Programování pro fyziky 2022/23

Úvod a cíl

Toto je první z průvodních textů k přednášce Programování pro fyziky. K samostudiu jsou určeny Poznámky k přednášce a také Poznámky k paralelní přednášce. Pro další praktické zkušenosti lze zajít na stránky cvičení (přednostně samozřejmě na stránky 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++, Fortran, (pro výuku) Pascal, ř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. 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 1 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)
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)
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

program a
real(4) MAX4; real(8) MAX8
print *,'ahoj'
print *,1,1.
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
  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.

Pascal

Dost bylo Fortranu, nyní Pascal pro totéž. Překladač Free Pascal je dobré používat v jeho vývojovém prostředí (IDE, Integrated Development Environment) Lazarus, to ale můžeme i nechat stranou a zůstat s Free Pascalem ve Visual Studio Code. Ekvivalent výše vyvinutého fortranského programu vychází v Pascalu takto:

program b;
uses Math;
var MAX : longint;  // longint neboli 4bytový integer
begin
  writeln('ahoj');
  writeln(1,1.0);
  MAX:=MaxLongint;
  writeln(MAX,' ',MAX+1);
  writeln(1.0/3.0,single(1.0)/single(3.0),double(1.0)/double(3.0));
  writeln(MaxSingle,MaxDouble);
  writeln(MaxSingle*2,MaxDouble*2);
end.

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

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

  • některé z použitých jmen nejsou k dispozici v holém Pascalu, proto je třeba připojit modul (jednotku, unit) Math,

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

  • reálný literál 1.0 před pravou závorkou píšeme ze zcela prapodivné příčiny – zkuste si sami nevýznamnou nulu umazat a program přeložit,

  • Pascal neoddělí bez pokynu mezerou dvě po sobě zapsaná celá čísla, proto tak činíme vložením literálu o jedné mezeře sami,

  • 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 Pascalu musíme říci konverzní funkcí single.

Uložíme-li zdrojový soubor jako b.pas, můžeme jej přeložit příkazem volajícím Free Pascal Compiler fpc b.pas a výsledek spustit příkazem b. Nerozumí-li systém volání překladače fpc, je třeba přidat cestu k němu (u mě: C:\lazarus\fpc\3.2.2\bin\x86_64-win64) do proměnné PATH (Windows: Properties of This PC – Advanced system settings – Environment variables nebo Win+R a systempropertiesadvanced). Výstup:

ahoj
1 1.000000000E+00
2147483647 2147483648
 3.333333433E-01 3.333333433E-01 3.3333333333333331E-001
 3.4028234664000001E+038 1.7976931348623157E+308
 6.8056469328000002E+038                    +Inf

a tedy rozdíly proti Fortranu jsou nejen kosmetické: na třetím řádku je výraz MAX+1 (na rozdíl od Fortranu) matematicky správně, ale mimo rozsah typu longint. Free Pascal zde v tichosti povýší provedení celočíselného výpočtu na 8bytový typ Int64. Pokud se pokusíme podobně jako ve Fortranu zvětšit i největší reálná čísla, uvidíme u 4bytové přesnosti single také tiché povýšení datového typu pro výpočet na 8bytový double, a v Linuxu (nikoliv ve Windows) můžeme vidět i povýšení 8bytového double na 10bytový typ extended. (Viz datové typy Pascalu.)

Python

Do třetice použijeme skriptovací jazyk Python. 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. Dnešní Anaconda obsahuje 5 GB softwaru, zatímco Miniconda je na počátku malinká a lze do ní snadno přidávat balíčky připravené ve webovém repozitáři velké Anacondy. Snadno přidávat balíčky lze i k holému Pythonu z webu nebo Microsoft Store; minimální sadu vhodnou pro naše potřeby zajistí příkaz: pip install numpy scipy matplotlib numba.

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, nebo je uložit jako skript c.py a vše spustit příkazem python c.py:

import sys
print('ahoj')
print(1,1.)
print(sys.maxsize,sys.maxsize+1)
print(1/3)
print(sys.float_info.max)

Podobně jako u Pascalu, 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 8 bytech ani jinde nebrzdí a můžeme se klidně podívat třeba na 1000 nul umocněním celočíselné desítky, print(10**1000). Samotný Python se nezajímá o 4bytovou reálnou aritmetiku, podpoří jen 8bytovou. Výstup skriptu:

ahoj
1 1.0
9223372036854775807 9223372036854775808
0.3333333333333333
1.7976931348623157e+308
inf

Matlab a Octave

Systém pro numerickou analýzu Matlab 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. Pro účely naší přednášky jsou oba systémy zástupné; instalace Octavu je časově i prostorově podstatně skromnější.

Skript d.m

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

C

Pro úplný dojem uvedeme ekvivalentní program i v jazyce C:

#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.,1.);
  MAX=INT_MAX; MAX4=FLT_MAX; MAX8=DBL_MAX;
  printf("%d %d\n",MAX,MAX+1);
  printf("%f %.17f\n",1.f/3.f,1./3.);
  printf("%e %.17e\n",MAX4,MAX8);
  printf("%e %.17e\n",MAX4*2,MAX8*2);
  return 0;
}

Z ukázky je znát technická náročnost programování v C, zde už jen při pouhém vypisování hodnot. Překládat můžeme překladačem gcc z balíčku GCC; příkaz gcc e.c vytvoří a.exe, příkazem gcc -o e e.c se output překladače nazve e.exe. Na jeho výstupu pak bude:

ahoj
1 1.000000 1.00000000000000000
2147483647 -2147483648
0.333333 0.33333333333333331
3.402823e+38 1.79769313486231571e+308
inf inf

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

Gnuplot

Gnuplot budeme používat pro grafickou vizualizaci našich numerických výsledků, především pro jeho 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:

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 nebo 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 f.gp a v gnuplotu spustit příkazem load 'f.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 extenze Code Runner. Všechny ukázkové programy jsou sbaleny zde: ukazky-01-overflow.zip.

Poznámky

GNU Compiler Collection a mnohé další vývojářské systémy vznikají v Linuxu a po jejich přenosu (portování) do Windows mívají různé zvláštnosti. Zdá se, že GCC i dalším dělá lépe, když nejsou instalovány v adresáři, v jehož názvu je mezera či diakritika; přitom domácí adresáře uživatelů ve Windows to mohou obsahovat snadno. Ukazuje se tak jako robustnější řešení instalovat rovnou do kořenového adresáře disku C:, např. tedy GCC do C:\gcc, Lazarus do C:\lazarus, Gnuplot do C:\gnuplot a GNU Octave do C:\Octave. Dále se ukazuje, že ačkoliv se GCC snaží umístit odkaz na své adresáře s překladači a dalšími nástroji do proměnné prostředí (environment variable) PATH, Windows to mohou a nemusí vzít bezprostředně na vědomí. Viděli jsme dnes na vlastní oči, jak po umístění adresářů do proměnné PATH v Properties of This PC–Advanced system settings–Environment variables, otevření nového příkazového okna a zadání příkazu path systém zadané adresáře na vědomí nevzal. Pomohlo nám univerzální řešení podobně nejasných situací – restart počítače. Pomůže ovšem i zavření okna Environment Variables; my sice zavírali okno Edit environment Variable, ale nikoliv okno Environment Variables.

Je-li adresář s gfortran.exe a gcc.exe pro překladače Fortranu a C (a podobně s fpc.exe pro Free Pascal, python.exe pro Python a wgnuplot.exe pro gnuplot) uveden v proměnné PATH (a ověřen výpisem příkazu path), lze je volat odkudkoliv. Případné hlášení o nenalezení f951.exe znamenalo opět jen problém s cestou k němu. GCC kromě C:\gcc\bin potřebuje v proměnné PATH i C:\gcc\libexec\gcc\x86_64-w64-mingw32\12.2.0. (Míněno teď, když aktuální verze GCC je 12.2.0.)

Chce-li programátor vyvíjet soubor a.f90 na disku D: v adresáři D:\MujAdresar, může si přepnout pracovní adresář na disk D: sekvencí příkazů D: a cd \MujAdresar nebo sloučit oba kroky do jednoho příkazem cd /d D:\MujAdresar. Nápovědu k příkazu ve Windows lze získat volbou /?, zde tedy příkazem cd /? – tak se člověk dozví o volbě /d neboli /D. (Anebo vhodným dotazem na Googlu.)

Mnoho voleb mají i překladače. Překladače GCC zobrazují některé své volby po gfortran --help a mnoho dalších po upřesnění, např. gfortran --help=common atp. (Možnosti pro volbu --help= naznačuje gfortran --help.) Čtenář těchto helpů se tak např. dozví, že chce-li vytvořit spustitelný soubor jiného jména než a.exe, musí si o to říct explicitně; pro překlad souboru becko.f90 do becko.exe tak je třeba psát gfortran -o becko becko.f90. Free Pascal na rozdíl od GCC překládá do exe-souborů stejnojmenných jako je zdroj. Své volby vypisuje po samotném fpc nebo fpc -? a bez stránkování po fpc -h. Stránkování výpisu nápovědy pro gfortran lze dosáhnout univerzálním způsobem (funkčním ve Windows i v Linuxu), tzv. přesměrováním výstupu na vstup tzv. filtru more, konkrétně: gfortran --help | more. Vypíše se stránka, zobrazí se slůvko More a dál můžeme mačkat Enter, mezerník nebo pro ukončení klávesu q nebo kombinaci Ctrl+C (Cancel). Kombinace Ctrl+C je značně univerzální řešení při potřebě přerušit výpis v příkazovém okně.

Budiž nám nadějí, že problémy s nastavením proměnné PATH člověk obvykle řeší jednorázově a pak už nemusí. (Dokud si s tím nesedne k jinému počítači.)