Cvičení 1: Rychlý úvod do SageMath

Tomáš Kalvoda, KAM FIT ČVUT, 2019

SageMath, zkráceně Sage, je open-source počítačový algebraický systém. Poněkud zjednodušeně lze říci, že jde o distribuci matematického softwaru spojenou jednotným Pythonovským rozhraním. První verze byla zveřejněna v roce 2005 Williamem Steinem (tou dobou profesor matematiky na Univerzitě ve Washingtonu).

Jak Sage získat?

K dispozici je několik v podstatě ekvivalentních možností jak získat svou vlastní kopii Sage:

instalátor pro Windows je k dispozici na Githubu,
binárky pro Linux,
binárky pro Apple Mac OSX,
vedle toho je možné využít Docker kontejnery nebo Virtualbox obrazy.

V neposlední řadě lze využít online službu COCALC, kterou lze v základní verzi používat zdarma. Možnosti COCALC jdou daleko za hranice Sage samotného a nebudeme se zde jimi podrobněji zabývat.

Notebooky

Se Sage si lze hrát pomocí několika různých rozhraní. Pro jistotu zde nejprve všechny zmíníme, aby se případný čtenář měl možnost zorientovat:

příkazová řádka (vhodné pro malé experimenty),
skripty (vhodné pro neinteraktivní výpočty),
Sage Notebook (aktuálně zastaralé) rozhraní používající webový prohlížeč,
Sage Worksheet (analog Sage Notebooku používaný v COCALC),
Jupyter Notebook.

Během semestru se budu soustředit prakticky pouze na poslední možnost.

Jupyter Notebook

Ze všech těchto možností jsem se rozhodl v tomto kurzu používat Jupyter notebook (zkráceně Jupyter; dříve též IPython notebook). Zde je stručný seznam jeho vlastností:

Jupyter opět běží v prohlížeči, má čisté a jednoduché uživatelské rozhraní.
Jupyter notebook v podstatě není nic jiného než interaktivní záznam vyhodnocování kusů kódu, zaznamenávání a zobrazování výstupů a anotovaný pomocí Markdown textu. Na pozadí jde o prostý JSON soubor.
Jupyter není vázaný na Sage, nebo Python. V tomto prostředí můžete pomocí "jader" přistupovat v podstatě k libovolnému jinému prostředí. Jupyter je v tomto směru agnostický. Velmi intenzivně je používaný v oblastech jako datové vědy (kde se pracuje s R, Python, Julia,...). Klidně byste ho ale mohli použít třeba na experimentování s C nebo C++. Oficiální seznam jader najdete na Githubu.
Fakultní Gitlab (ale samozřejmě i Github) umí Jupyter notebooky pěkně a automaticky zobrazovat, notebooky lze exportovat do spousty dalších formátů, např. pdf a html.

Vedle toho na Jupyter notebooky můžete narazit i v nečekaných oblastech (viz např. zde a zde).

Pokud máte na svém počítači vlastní instalaci Sage, můžete Jupyter notebook server spustit z příkazové řádky příkazem

$ sage --notebook=jupyter

Číselné množiny v Sage

Ústředními objekty, které studuje lineární algebra, jsou vektorové prostory a lineární zobrazení mezi těmito prostory. Vektorový prostor si lze představit jako množinu objektů (vektorů), které lze sčítat a násobit čísly, přičemž tyto operace mají jisté garantované vlastnosti.

V tomto cvičení se podrobněji podíváme na výše zmíněná čísla, resp. jejich zobecnění. V lineární algebře zavádíme pojem tělesa, jako množiny objektů, které můžeme vzájemně sčítat a násobit, přičemž pro obě operace platí asociativní a komutativní zákon, distributivita násobení vůči sčítání, máme speciální objekty " $0$ " a " $1$ " se známými vlastnostmi ( $a + 0 = a$ a $1 \cdot a = a$ ), dále ke každému prvku $a$ existuje prvek k němu opačný, ozn. $-a$ , splňující $a + (-a) = 0$ a ke každému nenulovému prvku $a$ existuje prvek k němu inverzní, ozn. $a^{-1}$ , splňující $a \cdot a^{-1} = 1$ .

Nezapomněli jsme na operace odčítání a dělení? Ty bychom s čísly přece také rádi prováděli. Tyto operace jsou ale definovány pomocí opačných a inverzních prvků. Přesněji, pro libovolné dva prvky $a$ a $b$ tělesa klademe $a - b := a + (-b)$ a je-li $b$ nenulový pak bychom mohli položit $a / b := a \cdot b^{-1}$ .

Jinak řečeno, v tělese můžeme počítat přesně tak jak jsme doposud byli zvyklí s reálnými (nebo komplexními) čísly. Existují ovšem i další tělesa, která jsou velmi důležitá v aplikacích. Jak uvidíme později, pojem vektorového prostoru lze vybudovat nad libovolným tělesem. Většina úvah (nebo postupů, či algoritmů) je na této vyšší úrovni shodná, jen aritmetika na nižší úrovni probíhá v jistém zadaném tělese.

V tomto souboru si projdeme a připomeneme různé číselné množiny, které jsou nám dostupné v Sage. Důležitým aspektem tohoto výkladu je objektový pohled na věc (viz BI-PA2), v tomto směru velmi blízký matematické abstrakci popsané v předchozím odstavci.

Dost ale povídání, pusťme se do toho.

Celá čísla

Množina celých čísel $\mathbb{Z}$ spolu s jejich standardním sčítáním a násobením netvoří těleso, drtivá většina nenulových prvků (vyjma $1$ a $-1$ ) totiž mezi celými čísly nemá inverzní prvek (tj. např. pro $x = 2 \in \mathbb{Z}$ neexistuje celé číslo $y$ splňující $x \cdot y = 1$ ). Přesto ale touto množinou začneme, jde o velmi jednoduchou a známou číselnou množinu.

V Sage celá čísla najdete pod symbolem ZZ (notace inspirovaná matematickým symbolem dvojitého "Z", tj. $\mathbb{Z}$ )

ZZ

Integer Ring

Poznámka: Kód v buňce spustíte stisknutím SHIFT+ENTER (stejně jako v Mathematica). Funkce print vypíše výstup, který se zobrazuje pod buňkou (jinak kód v buňce vypíše výstup posledního příkazu, tj. poslední print je redundantní).

Poznámka: Zmínit jak vyvolat nápovědu pomocí ? a zdrojový kód pomocí ??.

print(ZZ.is_field())  # není to těleso
print(ZZ.is_finite()) # není konečné, tedy nemá konečný počet prvků

False
False

Mimo to libovolný "celočíselný literál" je v Sage automaticky interpretován jako celé číslo, resp. prvek ZZ. Nejde tedy o Pythonovský int, ale o Sageovské objekty, které mají (na rozdíl od Pythonovského integeru) spoustu zajímavých metod (a na rozdíl od obyčejného integeru podporují opravdu velká celá čísla, co paměť dovolí).

# Sage je inteligentní a rozpozná celé číslo
x = 60
print(x in ZZ)

# Občas může být nutné Sage donutit interpretovat objekt tak jak chceme
x = ZZ(60.0)
print(x in ZZ)

True
True

# Další metody můžete prozkoumat stisknutím klávesy TAB za tečkou.
print(x.factor())         # symbolická faktorizace
print(x.divisors())       # dělitelé
print(x.is_prime())       # je to prvočíslo?
print(x.bits())           # bitová binární reprezentace
print(x.digits())         # cifry
print(x.previous_prime()) # předcházející prvočíslo

2^2 * 3 * 5
[1, 2, 3, 4, 5, 6, 10, 12, 15, 20, 30, 60]
False
[0, 0, 1, 1, 1, 1]
[0, 6]
59

Mezi celými čísly můžeme přirozeně provádět známé algebraické operace. Kdybychom požadovanou operací celá čísla opustili, tak Sage nevrátí chybu, ale vrátí výsledek v nejmenším možném rámci, kde daná operace má smysl (coertion model).

print(2 + 3)
print(5 - 3)
print(7 * 2)
print(8 / 2)
print(9 / 2)
print(9 / 2 in ZZ)
print(sqrt(4) in ZZ)
print(sqrt(5))
print(sqrt(5) in ZZ)

5
2
14
4
9/2
False
True
sqrt(5)
False

Tím se přirozeně dostáváme k následující kapitole.

Racionální čísla

Množina racionálních čísel spolu se známými operacemi jejich sčítání a násobení už tvoří těleso. k libovolnému nenulovému racionálnímu číslo $x = p/q$ existuje jeho inverzní prvek vůči násobení, $y = q/p$ , splňující $x \cdot y = 1$ .

Množina celých čísel v Sage se značí QQ, analogicky k matematickému zápisu $\mathbb{Q}$ .

QQ

Rational Field

Sage o této (mimojiné) například ví:

print(QQ.is_field())   # je to těleso
print(QQ.is_finite())  # není konečné, nemá nekonečný počet prvků

True
False

Podobně jako celá čísla mají racionální čísla tu pěknou vlastnost, že je lze v principu reprezentovat exaktně pomocí počítače (tj. bez ztráty přesnosti). Jsme jen omezeni množstvím paměti, kterou máme k uložení čitatele a jmenovatele k dispozici. To je činní atraktivní -- výpočet v tomto tělese je absolutně přesný (dokud máme paměť). Na typické školní příklady, kde se vyskytují pěkná celá čísla a zlomky si vystačíme právě s $\mathbb{Q}$ .

Podobně jako u celých čísel Sage "inteligentně" rozezná racionínální čísla.

print(4/3 in QQ)
print(4/3)

True
4/3

Je důležité si uvědomit, že tyto čtyři třetiny máme nyní reprezentované absolutně přesně. Neuchylujeme se zde k desetinným číslům s plovoucí čárkou, v kterých toto číslo nelze přesně reprezentovat.

S racionálními čísly můžeme opět přirozeně počítat tak jak jsme zvyklí.

print(4/3 + 1/2)
print(4/3 * 2/3)

11/6
8/9

print(sqrt(4/9) in QQ) # pokud to lze, vrátí Sage výsledek jako racionální číslo
print(sqrt(2/3) in QQ)

True
False

Reálná a komplexní čísla

Reálná čísla $\mathbb{R}$ a komplexní čísla $\mathbb{C}$ s příslušnými operacemi sčítání a násobení tvoří tělesa. Na rozdíl od celých a racionálních čísel už ale s nimi nejsme pomocí současných počítačů schopni skutečně pracovat. Známá čísla s pohyblivou desetinnou čárkou (floating point numbers) jsou pouze konečnou podmnožinou množiny racionálních čísel. Tento nedostatek lze částečně obejít pomocí symbolických výrazů (viz níže). Řada problémů numerických výpočtů (na které později narazíme) právě pramení z této principiální nemožnosti reálná čísla přesně reprezentovat.

Je dobré si uvědomit, že při počítání s těmito tzv. strojovými čísly strácíme řadu důležitých vlastností. Nejen, že většinu čísel pomocí nich nemůže přesně reprezentovat, ale i při každé operaci s nimi (sčítání, násobení, atd.) děláme většinou chybu způsobenou zaokrouhlováním, tedy převodem skutečného výsledku na číslo reprezentovatelné v dané přesnosti. Strojová čísla navíc nesplňují ani asociativní zákony...

V Sage lze reálná čísla s pohyblivou desetinnou čárkou (a i komplexní -- dvojice reálných) reprezentovatel v libovolné přesnosti (tj. nemáme jenom single a double, můžeme přesně říct kolikabitové reprezentace čísel budeme uvažovat a potom v tomto rámci počítat, se všemi s tím spojenými chybami). Ve výchozím nastavení se chovají prakticky stejně jako 64 bitový double, který má 53 bitovou mantisu.

print(RR)
print(RealField(666))

Real Field with 53 bits of precision
Real Field with 666 bits of precision

Veškerá aritmetika a zaokrouhlování se pak odehrává v dané přesnosti. Například, jednoduchá ukázka:

T = RealField(4)
x = T(pi)
y = T(1/3)
print(x)
print(y)

3.2
0.34

Co přesně máme v proměnné x?

x.sign_mantissa_exponent()

(1, 13, -2)

n(13 * 2^(-2))

3.25000000000000

Ukázkové násobení:

x * y

1.1

Srovnejte tento výsledek s lepší aproximací

n(pi / 3, digits=100)

1.047197551196597746154214461093167628065723133125035273658314864102605468762069666209344941780705689

Při našich školních výpočtech v zásadě platí, že reálným a komplexním číslům, resp. jejich nepřesné reprezentaci v počítači pomocí strojových čísel, se chceme vyhnout. Problematikou výpočtů pomocí strojových čísel se zabývá numerická matematika. V našich příkladech takřka vždy můžeme počítat v racionálních číslech, nebo využít symbolických výrazů.

Symbolické výrazy

Sage, jakožto správný počítačový algebraický systém, umožňuje pracovat se symbolickými výrazy, tj. v absolutní přesnosti. Čtenáři je asi jasné co reprezentují následující výrazy:

Poznámka: Příkaz show je sofistikovanější verze příkazu print. show dokáže pěkně vysázet symbolické výrazy, matice, atd.

show(sqrt(2))
show(sin(pi/7))
show(e^7 - 10)
show(sqrt(pi)/3)

\newcommand{\Bold}[1]{\mathbf{#1}}\sqrt{2}

\newcommand{\Bold}[1]{\mathbf{#1}}\sin\left(\frac{1}{7} \, \pi\right)

\newcommand{\Bold}[1]{\mathbf{#1}}e^{7} - 10

\newcommand{\Bold}[1]{\mathbf{#1}}\frac{1}{3} \, \sqrt{\pi}

Také bychom mohli vytvářet symbolické proměnné, tím se budeme zabývat později.

Je důležité uvědomit si, rozdíl ve výše uvedených výrazech a počítáním v konečné přesnosti:

print(sqrt(2.0))
print(sin(pi/7.0))
print(e^7.0 - 10)
print(sqrt(pi)/3.0)

1.41421356237310
sin(0.142857142857143*pi)
1086.63315842846
0.333333333333333*sqrt(pi)

Tj. jakmile do našeho výrazu zamotáme strojové číslo, tak ztrácíme přesnost. Tomu bychom se při našem počítání měli vyhnout, tj. místo případného 0.5 je vždy lepší napsat 1/2, či dělit dvojkou.

Konečná tělesa prvočíselného řádu (modulární aritmetika)

Novým typem tělesa, s kterým se milý čtenář ještě asi blížeji nesetkal, jsou tzv. konečná tělesa. Tedy tělesa s konečným počtem prvků. Tato tělesa mají tu podstatnou výhodu, že je můžeme v celé obecnosti reprezentovat v počítači (tj. všechny prvky udržet v paměti v absolutní přesnosti). Zároveň se jedná o tělesa v pravém slova smyslu, jsou splněny všechny požadavky na operace v tělese kladené.

Existují dva typy konečných těles: konečná tělesa, která mají prvočíselný počet prvků, a konečná tělesa mající počet prvků roven mocnině prvočísla (např. $2^8$ ). V BI-LIN si vystačíme s prvním uvedeným typem (konstrukce druhého typu vyžaduje techničtější manipulace s polynomy).

Nyní ke konstrukci tělesa s $p$ prvky ( $p$ je prvočíslo): je to množina $\{0,1,\ldots,p-1\}$ s operacemi sčítání a násobení modulo $p$ .

Například vezměme $p = 5$ :

T = GF(5) # Galois Field, na počest francouzského matematika Evariste Galois.
print(T)

Finite Field of size 5

Výpis prvků a vlastností.

print(T.is_field())  # je to těleso
print(T.is_finite()) # je konečné
print(T.order())     # řád, tj. počet prvků
print(list(T))       # výpis všech prvků

True
True
5
[0, 1, 2, 3, 4]

Pokud chceme z nějakého důvodu iterovat přes všechny prvky tělesa, můžeme k tomu použít Pythonovský for cyklus. Vypočtěme si třeba všechny kvadráty prvků z tohoto tělesa.

for x in T:
    print(x^2)

Chceme-li s "čísly" $0$ až $4$ pracovat skutečně v rámci tohoto tělesa, musíme k tomu použít konstruktor T (jinak je Sage přirozeně bude chápat jako celočíselné literály).

x = T(2)
y = T(3)
print(x + y)
print(x * y)

0
1

Všimněte si, že používáme obyčejné operátory + a *, to jak se tyto operátory chovají je dáno typem objektů mezi kterými působí. Nikde není potřeba explicitně psát jaké modulo používáme, to je dáno v definici T.

T je těleso a tedy každý nenulový prvek by v něm měl mít svůj inverzní prvek. Zkusme ho ručně najít pro $x = 2$ : hledáme takový prvek, který po vynásobení s $x$ dá $1$ .

for t in T:
    if t * x == 1:
        print(t)

Skutečně:

T(3) * x

Můžeme samozřejmě použít přímo i metody v Sage:

print(x^(-1))
print(1 / x)

3
3

U takto malých těles (řádu malého prvočísla) typicky není problematické příslušný inverzní prvek metodou pokus-omyl uhodnout (a to nám v BI-LIN bude stačit, podrobněji se touto úlohou budete zabývat v BI-ZDM a BI-BEZ). Existuje ale efektivní způsob jak tyto inverze hledat: rozšířený Euklidův algoritmus, EEA.

Přesněji, inverzi $x$ modulo $p$ nalezneme tak, že pomocí EEA spočteme Bézoutovy koeficienty $k,\ell\in\mathbb{Z}$ splňující $x \cdot k + p \cdot \ell = 1$ , což znamená $x \cdot k \equiv 1 \pmod p$ a proto hledanou inverzí je $k \bmod p$ . Zpět k předchozímu příkladu:

_, k, l = xgcd(2, 5)
print(2*k + 5*l)
print(mod(k, 5))

1
3

Podrobněji se těmito tělesy zabýváme ve vedlejším noteboooku.

Příklad

K pochopení jakým způsobem se tělesa v Lineární algebře používají použiejem jednoduchý programátorský příklad. Pointou je, že řada lineárně algebraických konstrukcí vlastně nezávísí na tom, v jakém tělese počítáme. Stačí nám jistota, že dané algebraické operace mají požadované vlastnosti.

Uvažme následující funkci.

def f(x, y):
    if y != 0:
        return 2 * x / y + y
    else:
        return 3*x

Python je dynamicky typovaný jazyk. Zavoláme-li funkci f na argumentech, které podporují uvedené operace, tak vše dobře projde (duck typing). Například:

# racionální čísla
x = 2/3
y = 4/3
f(x, y)

7/3

x = 1/2
y = 0
f(x, y)

3/2

# strojová čísla (netvoří těleso!)
x = RR(pi)
y = RR(e)
f(x, y)

5.02973652804089

x = RR(2.3)
y = 0
f(x, y)

6.90000000000000

# symbolické výrazy
x = sqrt(2)
y = pi
f(x, y)

pi + 2*sqrt(2)/pi

x = sin(3)
y = sin(0)
f(x, y)

3*sin(3)

# konečná tělesa
T = GF(5)
x = T(2)
y = T(3)
f(x, y)

x = T(3)
y = T(0)
f(x, y)