Knihovna mpmath: práce s numerickými hodnotami s plovoucí řádovou čárkou
Co se dozvíte v článku
- Knihovna mpmath: práce s numerickými hodnotami s plovoucí řádovou čárkou
- Reprezentace numerických hodnot v systému plovoucí řádové čárky
- Role báze: použít dvojku nebo desítku?
- Reprezentace numerických hodnot v knihovně mpmath
- Ukázka přednastavených numerických hodnot
- Vytvoření projektu, který bude využívat knihovnu mpmath
- Základní balíček mpmath
- Práce s reálnými čísly
- Kontext, nastavení přesnosti a/nebo počtu platných cifer
- Praktické příklady změny přesnosti či počtu platných cifer
- Kontext a změna zaokrouhlovacího režimu
- Funkce pro zpracování reálných hodnot
- Detekce situace, kdy je výsledkem komplexní číslo
- Iterativní výpočet hodnoty π
- Nekonečna a hodnoty typu NaN
- Dělení nulou
- Rychlost výpočtů prováděných knihovnou mpmath
- Jednoduchý benchmark
- Repositář s demonstračními příklady
- Odkazy na Internetu
Připomeňme si, že v programovacím jazyku Python existují tři základní numerické datové typy, které jsou „vhodně“ pojmenovány takovým způsobem, že mohou mnohé programátory poměrně úspěšně zmást. Jedná se o datové typy int, float a complex. Numerický datový typ int umožňuje práci s celočíselnými typy s prakticky neomezeným rozsahem, na rozdíl od podobně pojmenovaného typu v jiných jazycích (C, C++, Java, Go, …), kde se jedná o celočíselný typ s předem specifikovaným omezeným rozsahem (například 32bitů či 64bitů).
Druhý standardní numerický datový typ se v Pythonu jmenuje float, což může být opět matoucí, protože interně se jedná o céčkovský typ double, nikoli float. A konečně numerický typ complex sestává z dvojice hodnot typu float, které reprezentují reálnou a imaginární složku komplexního čísla. Komplexní jednotka se v Pythonu zapisuje znakem j a nikoli i:
>>> type(10)
<class 'int'>
>>> type(10**10000000)
<class 'int'>
>>> type(3.14)
<class 'float'>
>>> type(1+2j)
<class 'complex'>
>>> type(1+2i)
File "<stdin>", line 1
type(1+2i)
SyntaxError: invalid decimal literal
Jak již ovšem bylo napsáno v perexu tohoto článku, existují oblasti, v nichž se využívají takové numerické hodnoty, které není vhodné či dokonce vůbec možné reprezentovat s využitím základních (primitivních) datových typů, mezi něž typicky patří výše zmíněný celočíselný typ či typy (různé varianty int) a standardní typy numerických hodnot s plovoucí řádovou čárkou (různé varianty float, v Pythonu se konkrétně jedná o typ odpovídající dvojité přesnosti podle normy IEEE 754).
Z tohoto důvodu se v některých vyšších programovacích jazycích setkáme s podporou takových numerických datových typů, které umožňují reprezentovat celočíselné hodnoty s libovolným rozsahem (to je právě případ Pythonu a jeho typu int) popř. hodnoty s desetinnou tečkou, které budou mít jak neomezený rozsah, tak i neomezenou (či volitelnou) přesnost (podporováno v Clojure i mnoha implementacích Scheme a LISPu). To ovšem není vše, protože se mnohdy setkáme i s datovým typem zlomek, kde jak čitatel, tak i jmenovatel mají neomezený rozsah (opět podporováno v Clojure, Scheme, LISPu a dalších vysokoúrovňových programovacích jazycích).
Kromě základní trojice typů popsaných v předchozím textu je však možné ve standardní knihovně Pythonu nalézt realizaci dvou dalších numerických datových typů. Jedná se o typ nazvaný Decimal a o typ Fraction (pozor, je umístěný v modulu fractions, tedy s „s“ na konci). My si dnes představíme knihovnu mpmath a její datový typ mpf, jenž umožňuje pracovat s hodnotami, které mají plovoucí řádovou čárku, ovšem obecně neomezený rozsah a přesnost.
Reprezentace numerických hodnot v systému plovoucí řádové čárky
V této kapitole si v rychlosti připomeneme, jak vlastně interně vypadá reprezentace numerických hodnot v systému plovoucí řádové čárky (anglicky floating point, protože se v angličtině při zápisu čísel namísto čárky používá tečka). Systém plovoucí řádové čárky je založen na tom, že vybraná podmnožina reálných čísel (konkrétně se jedná o vybraná racionální čísla) může být vyjádřena vztahem se třemi proměnnými a dvojicí konstant:
XFP = (-1)s × bexp-bias × m
přičemž význam jednotlivých symbolů použitých v tomto vztahu je následující:
- XFP značí reprezentovanou numerickou hodnotu z podmnožiny racionálních čísel (ta je zase podmnožinou čísel reálných). Díky vyhrazeným (speciálním) hodnotám je většinou možné rozlišit kladnou a zápornou nulu i kladné a záporné nekonečno, což je jeden z důležitých rozdílů oproti běžnému způsobu reprezentace celých čísel v dvojkovém doplňku. Také se většinou může uložit nečíselná hodnota: NaN – (Not a Number), která je výsledkem některých matematicky nedefinovaných operací, například 0/0 nebo 00.
- b je báze, někdy také poněkud nesprávně nazývaná radix, u normy IEEE 754 je to vždy dvojka, protože výpočty s bází dvě jsou pro číslicové obvody nejjednodušší. V minulosti se však používaly i jiné báze, například 8, 16 nebo i 10, s nimi se však již dnes prakticky nesetkáme (o to relevantnější jsou však v dnešním článku, v němž se kromě jiného zmíníme i o balíčku decimal, jenž používá bázi 10).
- exp je vždy kladná hodnota exponentu posunutého o hodnotu bias.
- bias je hodnota, díky které je uložený exponent vždy kladný. Tato hodnota se většinou volí dle vztahu:
bias=2eb-1-1, kde eb je počet bitů vyhrazených pro exponent. Pro specifické účely je však možné zvolit i jinou hodnotu. - m je mantisa, která je u formátů dle normy IEEE 754 vždy kladná (ovšem ne vždy tomu tak musí být).
- s je znaménkový bit nabývající hodnoty 0 nebo 1. Pokud je tento bit nulový, je reprezentovaná hodnota XFP kladná, v opačném případě se jedná o zápornou hodnotu. Vzhledem k tomu, že je jeden bit vyhrazen na uložení znaménka, je možné rozlišit kladnou a zápornou nulu (některé systémy reprezentace hodnot ovšem bit znaménka používají v negované podobě, jen výjimečně se setkáme s hodnotami –1 a 1).
Do tabulky zobrazené pod tímto odstavcem jsem se pokusil napsat všechny důležité formáty reprezentace exponentu a mantisy současně s uvedením báze. Ta je většinou dvojková, ovšem výjimkou je IBM 360 se šestnáctkovou bází a FP rutiny u osmibitových domácích mikropočítačů Atari s desítkovou bází:
| Počítač/norma/systém | Šířka (b) | Báze | Exponent (b) | Mantisa (b) |
|---|---|---|---|---|
| IEEE 754 half | 16 | 2 | 5 | 10+1 |
| IEEE 754 single | 32 | 2 | 8 | 23+1 |
| IEEE 754 double | 64 | 2 | 11 | 52+1 |
| IEEE 754 double extended | 80 | 2 | 15 | 64 |
| IEEE 754 quadruple | 128 | 2 | 15 | 112+1 |
| IEEE 754 octuple | 256 | 2 | 19 | 236+1 |
| IBM řady 7×x | 36 | 2 | 8 | 27 |
| IBM 360 single | 32 | 16 | 7 | 24 |
| IBM 360 double | 64 | 16 | 7 | 56 |
| HP 3000 single | 32 | 2 | 9 | 22 |
| HP 3000 double | 64 | 2 | 9 | 54 |
| CDC 6000, 6600 | 60 | 2 | 11 | 48+1 |
| Cray-1 | 64 | 2 | 15 | 48 |
| Strela | 43 | 2 | 7 | 35 |
| Apple II | 40 | 2 | 8 | 31+1 |
| ZX Spectrum | 40 | 2 | 8 | 31+1 |
| Atari (FP rutiny) | 48 | 10 | 7 | 40 |
| Turbo Pascal real | 48 | 2 | 8 | 39 |
Role báze: použít dvojku nebo desítku?
Většina formátů s plovoucí řádovou čárkou používá bázi s hodnotou 2 nebo 16. To ovšem v praxi vede k některým problémům, protože není možné přesně reprezentovat některá desetinná čísla. Ostatně si to můžete sami otestovat v interpretru jazyka Python na jednoduchém výpočtu:
1.2-1.0 0.19999999999999996
Mj. i kvůli tomuto problému se v některých oblastech používá báze nastavená na hodnotu 10. Tento formát velmi často nalezneme i u kalkulaček, navíc kombinovaný s mantisou uloženou v BCD (Binary Coded Decimal). Například slavná HP-35 používá vlastně velmi podobný formát, jako již zmíněné osmibitové mikropočítače Atari, ovšem s mnohem větší přesností:
- Exponent je uložen jako dvojice BCD číslic, může tedy mít hodnotu 0 .. 99 (to je obvyklé u většiny kalkulaček).
- Mantisa je uložena jako deset BCD číslic (jednodušší kalkulačky mají jen osm BCD číslic, popř. počítají s deseti číslicemi, ale zobrazují jich jen osm).
- Další dva čtyřbitové „nibbly“ jsou použity pro uložení znaménka mantisy a exponentu (zde se ovšem poněkud plýtvá místem).
Celkově je jedna numerická hodnota uložena v 56 bitech (14×4), což sice může vypadat jako plýtvání místem, ovšem musíme si uvědomit, že kalkulačky pracují pouze se dvěma či třemi registry, popř. v případě HP se zásobníkem se čtyřmi prvky (ale už ne s rozsáhlými vektory nebo maticemi). Použití báze nastavené na 10 a BCD má své výhody, především není nutné složitě kódovat vstup od uživatele na interní reprezentaci a naopak výsledky dekódovat zpět na displej (nehledě na to, že u báze rovné deseti je zcela přesně známá přesnost vyjádřená v počtu desetinných míst).
Na DPD (Densely Packed Decimal) jsou založeny i standardizované formáty decimal32, decimal64 a decimal128, které jsou taktéž popsány v normě IEEE 754–2008. Základní vlastnosti těchto formátů jsou shrnuty v následující tabulce:
| Formát | Šířka | Mantisa | Exponent |
|---|---|---|---|
| decimal32 | 32 bitů | 7 cifer | −95 až +96 |
| decimal64 | 64 bitů | 16 cifer | −383 až +384 |
| decimal128 | 128 bitů | 34 cifer | −6143 až +6144 |
Datový typ Decimal v Pythonu, který je dostupný ve standardním balíčku decimal, je postaven na podobném principu. Jedná se tedy o vhodným způsobem uložené hodnoty s plovoucí řádovou čárkou, přičemž základem exponentu (báze) je hodnota 10 a nikoli dnes obvyklejší hodnota 2. Navíc není typ Decimal omezen pouze na podporu HW (tedy na FPU jednotku), takže umožňuje specifikovat jak způsoby zaokrouhlení, tak i přesnost (počet cifer), chování při dělení nulou atd.
Reprezentace numerických hodnot v knihovně mpmath
Ústředním tématem dnešního článku je popis základních vlastností knihovny nazvané mpmath. V této knihovně se pracuje s několika reprezentacemi numerických hodnot. Základem, který nás bude dnes zajímat především, je typ pojmenovaný mpf, jenž je založen na reprezentaci reálné hodnoty ve formátu plovoucí řádové čárky. Ovšem k dispozici jsou i další numerické typy, především komplexní čísla mpc, což je dvojice hodnot typu mpf a dále mpi určený pro intervalovou aritmetiku (základní vlastnosti i způsoby použití si ukážeme příště).
Abychom pochopili vlastnosti knihovny mpmath, je důležité vědět, jak jsou vlastně interně uloženy numerické hodnoty typu mpf. Skutečně je použit formát plovoucí řádové čárky, tj. numerická hodnota je reprezentována několika složkami. Konkrétně se jedná o čtveřici (sign, man, exp, bc). Kulaté závorky jsou použity naschvál, protože se skutečně jedná o běžnou Pythonovskou n-tici (tuple):
| Název | Typ | Stručný popis |
|---|---|---|
| sign | int | znaménko reprezentované hodnotou 0 nebo 1 (pozor: nikoli –1 a 1) |
| man | int | mantisa, celočíselná kladná hodnota s neomezeným rozsahem |
| exp | int | exponent, celočíselná hodnota s neomezeným rozsahem |
| bc | int | počet bitů nutných pro uložení mantisy (bit count) |
Na tomto místě je vhodné si znovu uvědomit, že datový typ int má v Pythonu neomezený rozsah, takže (alespoň teoreticky) máme k dispozici i neomezený rozsah a přesnost hodnot typu mpf.
Ukázka přednastavených numerických hodnot
Podívejme se na příklady způsobu uložení některých numerických hodnot. Následující tabulka je přímo převzata z dokumentace ke knihovně mpmath, doplnil jsem jen způsob výpočtu:
| Hodnota | Mantisa | Exponent | Výpočet |
|---|---|---|---|
| 3 | 3 | 0 | 3×20=3×1 |
| 10 | 5 | 1 | 5×21=5×2 |
| 1,25 | 5 | –2 | 5×2-2=5×0,25 |
| –16 | –1 | 4 | –1×24=-1×16 |
Přímo ve zdrojových kódech knihovny mpmath je definováno několik základních konstant a taktéž speciálních hodnot. U speciálních hodnot je typické, že mají nulovou mantisu, ovšem exponent je nenulový (nulový exponent má pouze nula):
| Jméno | Hodnota | Sign | Mantisa | Exponent | Bit count |
|---|---|---|---|---|---|
| fzero | 0 | 0 | 0 | 0 | 0 |
| fone | 1 | 0 | 1 | 0 | 1 |
| fnone | –1 | 1 | 1 | 0 | 1 |
| ftwo | 2 | 0 | 1 | 1 | 1 |
| ften | 10 | 0 | 5 | 1 | 3 |
| fhalf | 1/2 | 0 | 1 | –1 | 1 |
| nan | NaN | 0 | 0 | –123 | –1 |
| finf | ∞ | 0 | 0 | –456 | –2 |
| fninf | -∞ | 1 | 0 | –789 | –3 |
Vytvoření projektu, který bude využívat knihovnu mpmath
V praktické části dnešního článku si ukážeme základní způsoby použití balíčku mpmath. Zejména se zaměříme na způsob práce s numerickými hodnotami typu mpf. V případě balíčku mpmath se, na rozdíl od výše zmíněného balíčku decimal, nejedná o balíček ze základní knihovny Pythonu, což znamená, že je nutné ho doinstalovat, například s využitím pip, pdm, uv atd. V článku použiji nástroj uv, ovšem naprosto stejným způsobem (jen záměnou příkazu) lze použít i pdm.
Nejprve si necháme vygenerovat nový projekt:
$ uv init mpmath-lib Initialized project `mpmath-lib` at `/tmp/ramdisk/mpmath-lib`
V adresáři s novým projektem by měl vzniknout mj. i soubor nazvaný pyproject.toml obsahující jak konfiguraci projektu, tak i všechny přímé závislosti (dependencies). Seznam závislostí je prozatím prázdný:
[project] name = "mpmath-lib" version = "0.1.0" description = "Add your description here" readme = "README.md" requires-python = ">=3.11" dependencies = []
V dalším kroku do projektu přidáme jedinou závislost, kterou je pochopitelně právě balíček mpmath:
$ uv add mpmath Using CPython 3.13.9 interpreter at: /usr/bin/python3.13 Creating virtual environment at: .venv Resolved 2 packages in 231ms Prepared 1 package in 318ms Installed 1 package in 7ms + mpmath==1.4.1
Projektový soubor pyproject.toml by měl nyní vypadat následovně:
[project]
name = "mpmath-lib"
version = "0.1.0"
description = "Add your description here"
readme = "README.md"
requires-python = ">=3.11"
dependencies = [
"mpmath>=1.4.1",
]
Ještě se pro úplnost přesvědčíme o tom, že tento balíček si s sebou nepřinesl žádné další závislosti:
$ uv tree Resolved 2 packages in 1ms mpmath-lib v0.1.0 └── mpmath v1.4.1
$ uv run mpf_add.py
I samotný interpret Pythonu se bude spouštět přes příkaz uv:
$ uv run python Python 3.13.9 (main, Oct 14 2025, 00:00:00) [GCC 15.2.1 20250808 (Red Hat 15.2.1-1)] on linux Type "help", "copyright", "credits" or "license" for more information. >>>
Samozřejmě je možné provést i „globální“ instalaci balíčku mpmath, popř. instalaci pouze pro aktuálně přihlášeného uživatele:
$ sudo pip install mpmath
popř.:
$ pip install --user mpmath
Základní balíček mpmath
Knihovna mpmath obsahuje jak základní balíček nazvaný mpmath, tak i několik podbalíčků, například calculus, matrices atd. Ovšem prozatím se věnujme právě balíčku mpmath:
import mpmath help(mpmath)
Předchozí dva příkazy ověří, je že možné balíček mpmath naimportovat a zobrazit si jeho nápovědu:
Help on package mpmath:
NAME
mpmath
PACKAGE CONTENTS
__main__
_interactive
_version
calculus (package)
ctx_base
ctx_fp
ctx_iv
ctx_mp
ctx_mp_python
function_docs
functions (package)
identification
libfp
libmp (package)
math2
matrices (package)
rational
usertools
visualization
Užitečnější bude zobrazení dostupných funkcí (mezi ně se započítají i konstruktory, což je vlastně správně):
import mpmath
for item in dir(mpmath):
if item[0] != "_" and callable(getattr(mpmath, item)):
print(item)
Výsledek:
FPContext MPContext MPIntervalContext absmax absmin acos acosh acot acoth acsc acsch agm airyai airyaizero airybi airybizero almosteq altzeta ... ... ... unitvector upper_gamma webere whitm whitw workdps workprec zeros zeta zetazero
Práce s reálnými čísly
V úvodních kapitolách jsme si řekli, že v knihovně mpmath jsou numerické hodnoty (konkrétně podmnožina reálných čísel) reprezentovány čtveřicí (sign, man, exp, bc). Ovšem v praktických příkladech se s touto čtveřicí nepracuje přímo. Namísto toho je definována třída mpf, jejíž instance se používají namísto „holé“ čtveřice:
import mpmath help(mpmath.mpf)
Zobrazená nápověda v tomto případě neukazuje, jak vlastně probíhá konstrukce objektu s jeho inicializací:
Help on class mpf in module mpmath.ctx_mp_python: class mpf(_mpf) | mpf(val=(0, 0, 0, 0), **kwargs) | | Method resolution order: | mpf | _mpf | mpnumeric | builtins.object | | Data and other attributes defined here: | | context = <mpmath.ctx_mp.MPContext object> | | ---------------------------------------------------------------------- | Methods inherited from _mpf: | | __abs__(self) | | __add__(self, other) | | __bool__(self) | | __divmod__(self, other) | | __eq__(self, other) | Return self==value. | | __float__(self)
Způsob konstrukce a inicializace hodnot typu mpf ukazuje následující příklad, ve kterém jsou použity všechny dostupné metody předávání inicializační hodnoty:
from mpmath import mpf
# výchozí hodnota
a = mpf()
print(type(a))
print(a)
# získání hodnoty z řetězce
b = mpf("1.23456")
print(type(b))
print(b)
# získání hodnoty z typu float
c = mpf(1.23456)
print(type(c))
print(c)
# získání hodnoty z typu int
d = mpf(42)
print(type(d))
print(d)
Výsledky:
<class 'mpmath.ctx_mp_python.mpf'> 0.0 <class 'mpmath.ctx_mp_python.mpf'> 1.23456 <class 'mpmath.ctx_mp_python.mpf'> 1.23456 <class 'mpmath.ctx_mp_python.mpf'> 42.0
S hodnotami typu mpf je možné provádět všechny základní aritmetické operace, což vlastně znamená, že jsou přetíženy příslušné operátory:
from mpmath import mpf
x = mpf("1.00")
y = mpf(2.00)
z = x + y
print(type(z))
print(z)
Výsledky:
<class 'mpmath.ctx_mp_python.mpf'> 3.0
Dokonce je možné ve výpočtech zkombinovat hodnoty typu mpf a běžné numerické hodnoty Pythonu:
from mpmath import mpf
x = mpf("1.00")
y = 2.0
z = x + y
print(type(z))
print(z)
w = y + x
print(type(w))
print(w)
Výsledky:
<class 'mpmath.ctx_mp_python.mpf'> 3.0 <class 'mpmath.ctx_mp_python.mpf'> 3.0
A právě zde lze s výhodou využít teoreticky neomezeného rozsahu hodnot typu mpf. Ověřme si komutativitu:
from mpmath import mpf
x = mpf("1.23e10000000")
y = 42
z = x + y
print(type(z))
print(z)
w = y + x
print(type(w))
print(w)
Výsledky budou v tomto případě korektní:
<class 'mpmath.ctx_mp_python.mpf'> 1.23e+10000000 <class 'mpmath.ctx_mp_python.mpf'> 1.23e+10000000
Kontext, nastavení přesnosti a/nebo počtu platných cifer
Většina matematických operací, s nimiž se v knihovně mpmath setkáme, je implementována formou metod objektu typu context. Dnes se seznámíme s kontextem nazvaným mp (instance třídy MPContext), jenž pracuje s hodnotami typu mpf, s nimiž jsme se setkali v předchozích kapitolách. Tento kontext je dostupný jako globální hodnota (pro přístup k ní je jen zapotřebí ho naimportovat) a umožňuje například nastavit počet cifer použitých při výpočtech atd.:
from mpmath import mp help(mp)
Nápověda ke třídě MPContext je poměrně dlouhá, ovšem zajímavá je informace o konstruktoru MPContext, která obsahuje výchozí hodnoty atributů kontextu (přesnost a taktéž zaokrouhlovací režim):
Help on MPContext in module mpmath.ctx_mp object: class MPContext(mpmath.ctx_mp_python.PythonMPContext, mpmath.ctx_base.StandardBaseContext) | MPContext(prec=53, rounding='n', trap_complex=False) | | Context for multiple precision floatng-point arithmetic. | | Method resolution order: | MPContext | mpmath.ctx_mp_python.PythonMPContext | mpmath.ctx_base.StandardBaseContext | mpmath.ctx_base.Context | mpmath.functions.functions.SpecialFunctions | mpmath.functions.rszeta.RSCache | mpmath.calculus.quadrature.QuadratureMethods | mpmath.calculus.inverselaplace.LaplaceTransformInversionMethods | mpmath.calculus.calculus.CalculusMethods | mpmath.matrices.matrices.MatrixMethods | mpmath.matrices.calculus.MatrixCalculusMethods | mpmath.matrices.linalg.LinearAlgebraMethods | mpmath.matrices.eigen.Eigen | mpmath.identification.IdentificationMethods | mpmath.calculus.optimization.OptimizationMethods | mpmath.calculus.odes.ODEMethods | mpmath.visualization.VisualizationMethods | builtins.object
Nastavení kontextu si pochopitelně můžeme zobrazit:
from mpmath import mp print(mp)
Výchozí hodnoty by měly vypadat následovně:
Mpmath settings: mp.prec = 53 [default: 53] mp.dps = 15 [default: 15] mp.rounding = 'n' [default: 'n'] mp.trap_complex = False [default: False]
Význam jednotlivých atributů je shrnut v tabulce:
| Atribut | Stručný popis |
|---|---|
| prec | přesnost vyjádřená počtem bitů mantisy (binární přesnost) |
| dps | přesnost vyjádřená počtem desetinných míst (decimální přesnost) |
| rounding | zaokrouhlovací režim (bude ukázán dále) |
| trap_complex | chování v případě, že výsledkem operace je komplexní hodnota |
Nejzajímavější je změna přesnosti. Můžeme modifikovat buď atribut prec nebo dps. Druhý z atributů se automaticky dopočítá (je asi logické, že oba atributy na sobě závisí). Otestujme si změnu přesnosti, napřed z 53 bitů na 100 desítkových číslic a poté na 100 bitů mantisy:
from mpmath import mp print(mp) print() mp.dps = 100 print(mp) print() mp.prec = 100 print(mp) print()
Z výsledků je patrné, jakým způsobem se dopočítá buď dps nebo naopak prec:
Mpmath settings: mp.prec = 53 [default: 53] mp.dps = 15 [default: 15] mp.rounding = 'n' [default: 'n'] mp.trap_complex = False [default: False] Mpmath settings: mp.prec = 336 [default: 53] mp.dps = 100 [default: 15] mp.rounding = 'n' [default: 'n'] mp.trap_complex = False [default: False] Mpmath settings: mp.prec = 100 [default: 53] mp.dps = 29 [default: 15] mp.rounding = 'n' [default: 'n'] mp.trap_complex = False [default: False]
Praktické příklady změny přesnosti či počtu platných cifer
Otestujme si nyní, jak se nastavení kontextu projeví na prováděných výpočtech. Nejdříve si necháme zobrazit výsledek podílu 1/7, a to při výchozím nastavení kontextu:
from mpmath import mpf
x = mpf("1.00")
y = mpf("7.00")
z = x / y
print(z)
Výsledek vypadá stejně, jako bychom pro výpočet použiti standardní typ float. To je ovšem očekávané chování, protože výchozí kontext má nastavenu přesnost na 53 bitů resp. 15 decimálních cifer:
0.142857142857143
V dalším příkladu počet dekadických cifer postupně zvyšujeme:
from mpmath import mp, mpf
x = mpf("1.00")
y = mpf("7.00")
for digits in range(20):
mp.dps = digits
z = x / y
print(z)
Výsledky ukazují jak zvyšující se počet cifer výsledku, tak i snahu knihovny mpmath o zaokrouhlování:
0.2 0.1 0.14 0.143 0.1429 0.14286 0.142857 0.1428571 0.14285714 0.142857143 0.1428571429 0.14285714286 0.142857142857 0.1428571428571 0.14285714285714 0.142857142857143 0.1428571428571429 0.14285714285714286 0.142857142857142857 0.1428571428571428571
Kontext je možné měnit i lokálně s využitím standardního manažeru kontextu (context manager). Kontext bude v tomto případě změněn je v rámci bloku with a poté bude obnovena jeho původní konfigurace:
from mpmath import workdps, mpf
x = mpf("1.00")
y = mpf("7.00")
for digits in range(20):
with workdps(digits):
z = x / y
print(z)
Výsledky by měly být stejné, jako v příkladu předchozím demonstračním příkladu:
0.2 0.1 0.14 0.143 0.1429 0.14286 0.142857 0.1428571 0.14285714 0.142857143 0.1428571429 0.14285714286 0.142857142857 0.1428571428571 0.14285714285714 0.142857142857143 0.1428571428571429 0.14285714285714286 0.142857142857142857 0.1428571428571428571
Větší počet cifer se projeví i v případě, že se počítá s velkými hodnotami. Příkladem je výpočet faktoriálu s výchozí konfigurací kontextu:
from mpmath import mpf
def factorial(n):
"""Výpočet faktoriálu ve smyčce."""
f = mpf("1")
for x in range(1, n + 1):
f *= x
return f
for n in range(31):
print(n, factorial(n))
Zobrazené výsledky do značné míry odpovídají výsledkům získaným při použití standardních hodnot typu float:
0 1.0 1 1.0 2 2.0 3 6.0 4 24.0 5 120.0 6 720.0 7 5040.0 8 40320.0 9 362880.0 10 3628800.0 11 39916800.0 12 479001600.0 13 6227020800.0 14 87178291200.0 15 1307674368000.0 16 20922789888000.0 17 355687428096000.0 18 6.402373705728e+15 19 1.21645100408832e+17 20 2.43290200817664e+18 21 5.10909421717094e+19 22 1.12400072777761e+21 23 2.5852016738885e+22 24 6.20448401733239e+23 25 1.5511210043331e+25 26 4.03291461126606e+26 27 1.08888694504184e+28 28 3.04888344611714e+29 29 8.8417619937397e+30 30 2.65252859812191e+32
Nic nám ovšem nebrání ve zvětšení počtu dekadických cifer, například na hodnotu 100:
from mpmath import mp, mpf
def factorial(n):
"""Výpočet faktoriálu ve smyčce."""
f = mpf("1")
for x in range(1, n + 1):
f *= x
return f
mp.dps = 100
for n in range(31):
print(n, factorial(n))
V tomto případě budou výsledky vypadat následovně:
0 1.0 1 1.0 2 2.0 3 6.0 4 24.0 5 120.0 6 720.0 7 5040.0 8 40320.0 9 362880.0 10 3628800.0 11 39916800.0 12 479001600.0 13 6227020800.0 14 87178291200.0 15 1307674368000.0 16 20922789888000.0 17 355687428096000.0 18 6402373705728000.0 19 121645100408832000.0 20 2432902008176640000.0 21 51090942171709440000.0 22 1124000727777607680000.0 23 25852016738884976640000.0 24 620448401733239439360000.0 25 15511210043330985984000000.0 26 403291461126605635584000000.0 27 10888869450418352160768000000.0 28 304888344611713860501504000000.0 29 8841761993739701954543616000000.0 30 265252859812191058636308480000000.0
Kontext a změna zaokrouhlovacího režimu
V rámci kontextu je nastaven i zaokrouhlovací režim. Knihovna mpmath rozlišuje pět zaokrouhlovacích režimů:
| Označení | Slovně | Význam |
|---|---|---|
| „n“ | nearest | zaokrouhlení k nejbližší hodnotě |
| „f“ | floor | zaokrouhlení směrem k −∞ |
| „c“ | ceiling | zaokrouhlení směrem k +∞ |
| „d“ | down | zaokrouhlení k nule |
| „u“ | up | zaokrouhlení od nuly |
Výběr zaokrouhlovacího režimu je snadný, což si ukážeme na příkladu, ve kterém se počítá hodnota 1/3 a taktéž –1/3:
from mpmath import mp, mpf
mp.dps = 2
x = mpf("1.00")
y = mpf("3.00")
print("mode\t 1/3\t -1/3")
print("-" * 22)
for rounding in "n", "f", "c", "d", "u":
mp.rounding = rounding
z1 = x / y
z2 = -x / y
print(rounding, "\t", z1, "\t", z2)
Výsledná tabulka ukazuje některé rozdíly mezi zaokrouhlovacími režimy:
mode 1/3 -1/3 ---------------------- n 0.33 -0.33 f 0.33 -0.34 c 0.34 -0.33 d 0.33 -0.33 u 0.34 -0.34
Funkce pro zpracování reálných hodnot
V rámci kontextu mp je možné provádět výpočty mnoha funkcí, které jako svůj vstup akceptují reálné hodnoty (s libovolným rozsahem a nebo přesností). Příkladem může být výpočet druhé odmocniny provedený na zvolený počet desetinných míst (povšimněte si, že funkci sqrt voláme tak, jakoby se jednalo o globálně definovanou funkci):
from mpmath import workdps, mpf, sqrt
x = mpf("2.00")
for digits in range(30):
with workdps(digits):
y = sqrt(x)
print(y)
Výsledky:
2.0 1.0 1.4 1.41 1.414 1.4142 1.41421 1.414214 1.4142136 1.41421356 1.414213562 1.4142135624 1.41421356237 1.414213562373 1.4142135623731 1.4142135623731 1.414213562373095 1.414213562373095 1.41421356237309505 1.414213562373095049 1.4142135623730950488 1.4142135623730950488 1.414213562373095048802 1.4142135623730950488017 1.41421356237309504880169 1.414213562373095048801689 1.4142135623730950488016887 1.41421356237309504880168872 1.414213562373095048801688724 1.4142135623730950488016887242
I když funkce sqrt jako svůj parametr akceptuje reálné číslo, může být výsledkem hodnota komplexní:
from mpmath import workdps, mpf, sqrt
x = mpf("-2.00")
for digits in range(30):
with workdps(digits):
y = sqrt(x)
print(y)
Výsledky:
(0.0 + 2.0j) (0.0 + 1.0j) (0.0 + 1.4j) (0.0 + 1.41j) (0.0 + 1.414j) (0.0 + 1.4142j) (0.0 + 1.41421j) (0.0 + 1.414214j) (0.0 + 1.4142136j) (0.0 + 1.41421356j) (0.0 + 1.414213562j) (0.0 + 1.4142135624j) (0.0 + 1.41421356237j) (0.0 + 1.414213562373j) (0.0 + 1.4142135623731j) (0.0 + 1.4142135623731j) (0.0 + 1.414213562373095j) (0.0 + 1.414213562373095j) (0.0 + 1.41421356237309505j) (0.0 + 1.414213562373095049j) (0.0 + 1.4142135623730950488j) (0.0 + 1.4142135623730950488j) (0.0 + 1.414213562373095048802j) (0.0 + 1.4142135623730950488017j) (0.0 + 1.41421356237309504880169j) (0.0 + 1.414213562373095048801689j) (0.0 + 1.4142135623730950488016887j) (0.0 + 1.41421356237309504880168872j) (0.0 + 1.414213562373095048801688724j) (0.0 + 1.4142135623730950488016887242j)
V knihovně mpmath najdeme i další užitečné funkce, například funkce goniometrické:
from mpmath import mpf, sin, cos, pi
x = pi / 4
print("π\t", pi)
print("π/4\t", x)
print("sin π/4\t", sin(x))
print("cos π/4\t", cos(x))
Výsledky:
π 3.14159265358979 π/4 0.785398163397448 sin π/4 0.707106781186547 cos π/4 0.707106781186548
Výpočty logaritmů atd.:
from mpmath import mpf, ln, log10
x = 100
print("x\t", x)
print("ln x\t", ln(x))
print("log x\t", log10(x))
print("ln 0\t", ln(0))
print("log 0\t", log10(0))
Výsledky vypadají následovně:
x 100 ln x 4.60517018598809 log x 2.0 ln 0 -inf log 0 -inf
Umocnění se ovšem nemusí provádět přes funkci, protože je přetížen operátor **:
from mpmath import mpf
x = mpf("2")
y = mpf("3")
print("x^y\t", x**y)
print("y^x\t", y**x)
print("x^-y\t", x**-y)
print("y^-x\t", y**-x)
Opět se podívejme na výsledky:
x^y 8.0 y^x 9.0 x^-y 0.125 y^-x 0.111111111111111
Nejsme ovšem omezeni pouze na celočíselné hodnoty:
from mpmath import mpf, workdps
x = mpf("1.2")
y = mpf("2.3")
with workdps(50):
print("x^y\t", x**y)
Výsledky:
x^y 1.5209567545525315862936214273099533693700312594907
Detekce situace, kdy je výsledkem komplexní číslo
V případě, že je v kontextu nastaveno zachytávání komplexních výsledků, dopadne výpočet druhé odmocniny odlišně:
from mpmath import workdps, mp, mpf, sqrt
x = mpf("-2.00")
mp.trap_complex = True
y = sqrt(x)
print(y)
Odmocnina ze záporné hodnoty je komplexní číslo, ovšem my jsme povolili zachytávání komplexních hodnot:
Traceback (most recent call last):
File "/tmp/ramdisk/mpmath-lib/func_sqrt_3.py", line 7, in <module>
y = sqrt(x)
File "/tmp/ramdisk/mpmath-lib/.venv/lib64/python3.13/site-packages/mpmath/ctx_mp_python.py", line 1166, in f
return ctx.make_mpf(mpf_f(x._mpf_, prec, rounding))
~~~~~^^^^^^^^^^^^^^^^^^^^^^^^^
File "/tmp/ramdisk/mpmath-lib/.venv/lib64/python3.13/site-packages/mpmath/libmp/libmpf.py", line 1794, in mpf_sqrt
raise ComplexResult("square root of a negative number")
mpmath.libmp.libmpf.ComplexResult: square root of a negative number
Iterativní výpočet hodnoty π
Zkusme si nyní provést nějaký výpočet, v němž by se do jisté míry využila konfigurovatelná přesnost výpočtů s hodnotami typu mpf. Pro jednoduchost jsem vybral jeden z dnes již nepoužívaných výpočtů konstanty π z nekonečné řady (nepoužívaný je tento postup z toho důvodu, že konverguje velmi pomalu, což ovšem pro nás bude výhoda). Jedná se o takzvaný Wallis product, což je forma řady, která vypadá následovně:
První realizace tohoto výpočtu bude založena na standardním numerickém datovém typu float, tedy na hodnotách s plovoucí čárkou a exponentem, jehož báze je rovna 2, ovšem počet bitů odpovídá přesně normě IEEE 754 (dvojitá přesnost, ovšem v Pythonu označená termínem float):
from math import pi
result = 2
for n in range(1, 1000):
m = 4 * n * n
u = m / (m-1)
result *= u
abs_error = pi - result
rel_error = 100.0 * abs_error / pi
print(result, "\t", abs_error, "\t", rel_error)
Po spuštění tohoto skriptu se bude vypisovat jak postupně se zpřesňující výsledek počítané nekonečně řady, tak i absolutní a relativní chyba vypočtená vůči konstantě π (ovšem typu float). Prvních 100 výsledků vypadá (ve zkrácené podobě) takto:
2.6666666666666665 0.4749259869231266 15.11736368432249 2.844444444444444 0.297148209145349 9.458521263277328 2.9257142857142853 0.21587836787550785 6.87162187079954 2.972154195011337 0.16943845857845607 5.393393646526529 3.0021759545569062 0.1394166990328869 4.437771360127803 3.0231701920013605 0.11842246158843261 3.769504026981831 3.038673628883419 0.10291902470637426 3.276014304043273 3.0505899960555105 0.0910026575342826 2.896704556215998 3.06003454712689 0.08155810646290318 2.5960751585572193 3.0677038066434985 0.0738888469462946 2.3519550461726566 3.074055160280442 0.06753749330935133 2.1497851808438146 3.0794013431678864 0.06219131042190673 1.9796108942017927 3.083963419231839 0.057629234357954306 1.8343955029339434 3.087902069831113 0.05369058375868008 1.7090243605366735 3.09133688859622 0.050255764993573315 1.5996906835183655 3.0943587232869296 0.04723393030286349 1.503502697871758 3.09703782174865 0.04455483184114328 1.4182243452292251 3.0994293567461395 0.04216329684365361 1.3420994219436762 3.101577263438271 0.04001539015152211 1.2737294284730982 3.103516961539233 0.03807569205056005 1.2119869202982834 3.105277322833356 0.036315330756437 1.1559528799808176 3.1068821173154406 0.0347105362743525 1.1048706850867482 3.108351092311807 0.033241561277986165 1.058111758696728 3.1097007888347385 0.031891864755054566 1.0151495840370264 3.110945166901499 0.0306474866882942 0.9755397999570167 3.1120960900117107 0.02949656357808239 0.9389047795352351 3.1131637044508227 0.028428949138970427 0.9049215564750451 3.1141567391252885 0.027435914464504663 0.8733122810544696 3.115082744697434 0.026509908892359046 0.8438366082269468 3.1159482858879586 0.025644367701834536 0.8162855764426229 3.1167590973076535 0.02483355628213957 0.7904766473706606 3.1175202106403295 0.024072442949463646 0.7662496575409568 3.1182360591387543 0.02335659445103877 0.7434634921351108 3.1189105640185164 0.02268208957127671 0.7219933349843635 3.119547206305518 0.022045447284275266 0.7017283815928417 3.12014908691642 0.021443566673373216 0.6825699267175955 3.120718977160606 0.02087367642918725 0.6644297568411868 3.1212593613990753 0.020333292190717778 0.6472287922969135 3.121772473245434 0.019820180344359173 0.6308959349555173 3.1222603264214372 0.019332327168355867 0.6153670860627161 3.1227247411658103 0.01886791242398278 0.6005843056203689 3.1231673669264293 0.01842528666336385 0.5864950900719064 3.1235897019320986 0.0180029516576945 0.5730517493133022 ... ... ... 3.132311790781417 0.009280862808376256 0.2954190384221622 3.1324201790229056 0.009172474566887523 0.2919689335409682 3.132526064841755 0.009066588748038118 0.2885984832463251 3.132629533910784 0.00896311967900898 0.28530496048770426 3.132730668036172 0.008861985553620944 0.28208576129354795 3.1328295453731676 0.008763108216625515 0.27893839790503083 3.132926240627508 0.008666412962285097 0.27586049236466975 3.133020825243655 0.008571828346138233 0.27284977052462517 3.133113367580835 0.008479286008958198 0.26990405644312926 3.133203933077802 0.008388720511991021 0.2670212671399492 3.1332925844071484 0.008300069182644698 0.26419940768452227 3.133379381619936 0.008213271969857328 0.2614365665921804 3.133464382281348 0.008128271308445179 0.25873091150621563 3.1335476415980024 0.008045011991790751 0.25608068514541454 3.1336292125375205 0.007963441052272646 0.25348420149802325 3.1337091459408963 0.007883507648896781 0.25093984224493776 3.133787490628162 0.00780516296163114 0.24844605339627468
Přepis do podoby využívající hodnoty typu mpf s nastavitelnou přesností může vypadat následovně (nejedná se ovšem o nejideálnější způsob implementace, protože neustále vytváříme nové objekty):
from mpmath import pi, mpf, mp
mp.dps = 100
result = mpf("2")
for n in range(1, 1000):
m = mpf(4 * n * n)
u = m / (m-1)
result *= u
abs_error = pi - result
rel_error = mpf("100.0") * abs_error / pi
print(result, "\t", abs_error, "\t", rel_error)
Výsledky, opět ve zkrácené podobě:
2.666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666666667 0.4749259869231265717959767166128362175305027327084391543082779256411497396195423319613681586754504013 15.11736368432248758992865953465900691495485560508989400124408316858837459507918128527263852466502082 2.844444444444444444444444444444444444444444444444444444444444444444444444444444444444444444444444444 0.2971482091453487940181989388350584397527249549306613765305001478633719618417645541835903808976726235 9.45852126327732009592390350363627404261851264542922026799368871316093290141779337095748109297602221 2.925714285714285714285714285714285714285714285714285714285714285714285714285714285714285714285714286 0.2158783678755075241769290975652171699114551136608201066892303065935306920004947129137491110564027822 6.871621870799529241521729318025881872407613006727197989936365533536959555744016038699123409918194273 2.972154195011337868480725623582766439909297052154195011337868480725623582766439909297052154195011338 0.16943845857845536998191775969673644428787234722091080963707611158219282351976908933098267114710573 5.393393646526505896149058354819943489429956070326042402475038002323260501073286134551490448170864023 3.00217595455690693785931881169976408071646166884262122357360452598547836643074738312833550928789024 0.1394166990328863006033245715797388034807077305324845974013400663223380398554616154996993160542268276 4.437771360127783733483897328100953019626218252854588285328321214467939900074026398536859038556428306 3.023170192001360832529663698494867326036157204988373819542650711481880313049144217975386806555637725 0.1184224615884324059329796847846355581610121943867320014322938808259360932370647806526480187864793434 3.769504026981824179172595910814945698085142856021403588022924859324359060214404205519634276588291441 3.038673628883419093209302999512789722579932370142160351950561740771530981321191111400901610691820482 0.1029190247063741452533403837667131616172370292329454690243828515362854249650178872271332146502965859 3.276014304043269431373481018049894137562502562975359503858939858602945516933452432214606760057974986 3.050589996055510932790515952452055564629265438260364902350367865245144828306764409876983577792572719 0.09100265753428230567212743082744731956790396111474091862457672706267157797944458875105124754954434865 2.896704556215988135025926041650089800847061396555655031325053348244525695431230677046820904215065083 ... ... ... 3.140806174053955615400178245660413819682477073377639887081848146490285046624665679563038899041276056 0.0007864795358376230624651376190890645146923259974659338930964458175313596615433190649959263008410116468 0.02503442115383543135665572790250708437519219017084062575656839466993861458089907038143020605426414408 3.140806960828458052887228080086451490599512877828957528285997312171810458876371558708340033093254234 0.0007856927613351855754153031930513935976565215461482926889472801360059474098374399196947922488628339596 0.0250093773436031128527446227895419644860248661964323326868051095642898137003923625706070131433885898
Nekonečna a hodnoty typu NaN
Při práci s reálnými hodnotami typu mpf je nutné počítat s tím, že výsledkem některých aritmetických a jiných operací bude kladné nekonečno, záporné nekonečno nebo i hodnota NaN (Not a Number). A tyto speciální hodnoty lze předávat do funkcí namísto „normálních“ reálných hodnot. Některé možnosti, které mohou v praxi nastat, jsou ukázány v dalším demonstračním příkladu:
from mpmath import workdps, mpf, log, sin
x = mpf("0.00")
y = mpf("0.00")
z = log(x)
print("log 0=", z)
inf = mpf("inf")
print("inf+ -inf=", inf+z)
print("inf/inf=", inf/inf)
print("log inf=", log(inf))
print("sin inf=", sin(inf))
Podívejme se na výsledky:
log 0= -inf inf+ -inf= nan inf/inf= nan log inf= inf sin inf= nan
Dělení nulou
V normě IEEE 754 jsou mj. definovány i výsledky v případě, že dělitelem je nula. Měl by se nastavit příznak dělení nulou, ovšem výsledkem by mělo být kladné nekonečno, záporné nekonečno nebo hodnota NaN (podle prováděné operace). Python při použití datového typu float ovšem v případě, že se dělí nulou, vyhodí výjimku typu ZeroDivisionError:
try:
1/0
except Exception as e:
print(e)
try:
-1/0
except Exception as e:
print(e)
try:
0/0
except Exception as e:
print(e)
Chování skriptu po spuštění:
division by zero division by zero division by zero
Podobně se chová i knihovna mpmath, ovšem s tím rozdílem, že výjimka po převodu na řetězec vrátí prázdný text. Nicméně výjimka je vždy vyhozena:
from mpmath import mpf
one = mpf("1")
zero = mpf("0")
try:
x = one/zero
except Exception as e:
print(type(e))
try:
x = -one/zero
except Exception as e:
print(type(e))
try:
x = zero/zero
except Exception as e:
print(type(e))
Výsledky:
<class 'ZeroDivisionError'> <class 'ZeroDivisionError'> <class 'ZeroDivisionError'>
Rychlost výpočtů prováděných knihovnou mpmath
Knihovna mpmath dokáže při výpočtech využít funkce poskytované knihovnou gmpy resp. gmpy2. Ta realizuje rozhraní mezi Pythonem a nativním céčkovým kódem, ve kterém se výpočty budou provádět. Pokud ovšem tato knihovna není nainstalována, jsou všechny výpočty prováděny přímo v Pythonu, což ovšem není ani jednoduché ani rychlé. Ostatně se o tom můžeme přesvědčit například pohledem na realizaci „pouhého“ součtu:
def mpf_add(s, t, prec=0, rnd=round_fast, _sub=0):
"""
Add the two raw mpf values s and t.
With prec=0, no rounding is performed. Note that this can
produce a very large mantissa (potentially too large to fit
in memory) if exponents are far apart.
"""
ssign, sman, sexp, sbc = s
tsign, tman, texp, tbc = t
tsign ^= _sub
# Standard case: two nonzero, regular numbers
if sman and tman:
offset = sexp - texp
if offset:
if offset > 0:
# Outside precision range; only need to perturb
if offset > 100 and prec:
delta = sbc + sexp - tbc - texp
if delta > prec + 4:
offset = prec + 4
sman <<= offset
if tsign == ssign: sman += 1
else: sman -= 1
return normalize(ssign, sman, sexp-offset,
sman.bit_length(), prec, rnd)
# Add
if ssign == tsign:
man = tman + (sman << offset)
# Subtract
else:
if ssign: man = tman - (sman << offset)
else: man = (sman << offset) - tman
if man >= 0:
ssign = 0
else:
man = -man
ssign = 1
bc = man.bit_length()
return normalize(ssign, man, texp, bc, prec or bc, rnd)
elif offset < 0:
# Outside precision range; only need to perturb
if offset < -100 and prec:
delta = tbc + texp - sbc - sexp
if delta > prec + 4:
offset = prec + 4
tman <<= offset
if ssign == tsign: tman += 1
else: tman -= 1
return normalize(tsign, tman, texp-offset,
tman.bit_length(), prec, rnd)
# Add
if ssign == tsign:
man = sman + (tman << -offset)
# Subtract
else:
if tsign: man = sman - (tman << -offset)
else: man = (tman << -offset) - sman
if man >= 0:
ssign = 0
else:
man = -man
ssign = 1
bc = man.bit_length()
return normalize(ssign, man, sexp, bc, prec or bc, rnd)
# Equal exponents; no shifting necessary
if ssign == tsign:
man = tman + sman
else:
if ssign: man = tman - sman
else: man = sman - tman
if man >= 0:
ssign = 0
else:
man = -man
ssign = 1
bc = man.bit_length()
return normalize(ssign, man, texp, bc, prec or bc, rnd)
# Handle zeros and special numbers
if _sub:
t = mpf_neg(t)
if not sman:
if sexp:
if s == t or tman or not texp:
return s
return fnan
if tman:
return normalize(tsign, tman, texp, tbc, prec or tbc, rnd)
return t
if texp:
return t
if sman:
return normalize(ssign, sman, sexp, sbc, prec or sbc, rnd)
return s
Jednoduchý benchmark
Rychlost nebo možná spíše pomalost výpočtů prováděných knihovnou mpmath si ověříme na velmi jednoduchém benchmarku, který je založen na výše zmíněném iterativním výpočtu hodnoty π. Porovnávat budeme rychlost výpočtů s typy float (standard), mpf (výchozí kontext) a mpf s přesností na 100 a 1000 desítkových cifer:
from time import time
from mpmath import mpf, mp, workdps
def wallis_pi_float(maxiter):
result = 2
for n in range(1, maxiter):
m = 4 * n * n
u = m / (m - 1)
result *= u
return result
def wallis_pi_mpf(maxiter):
result = mpf("2")
for n in range(1, maxiter):
m = mpf(4 * n * n)
u = m / (m - 1)
result *= u
return result
def benchmark(function, maxiter):
t1 = time()
function(maxiter)
t2 = time()
return t2 - t1
for maxiter in (100, 1000, 10000, 100000, 1000000):
dt1 = benchmark(wallis_pi_float, maxiter)
dt2 = benchmark(wallis_pi_mpf, maxiter)
with workdps(100):
dt3 = benchmark(wallis_pi_mpf, maxiter)
with workdps(1000):
dt4 = benchmark(wallis_pi_mpf, maxiter)
print(f"{maxiter:7} {dt1:9.7f} {dt2:9.7f} {dt3:9.7f} {dt4:9.7f}")
Výsledky v číselné podobě:
100 0.0000169 0.0007770 0.0007627 0.0019510 1000 0.0001135 0.0068927 0.0075638 0.0189612 10000 0.0011377 0.0634780 0.0710256 0.1901760 100000 0.0139859 1.1295485 1.1711011 2.7714369 1000000 0.1247005 8.7620215 10.3339429 23.6119809
V grafické podobě:
Obrázek 1: Rychlost výpočtu hodnoty π s různými typy i proměnným počtem iterací.
Z výsledků je patrné, že rychlost výpočtů s hodnotami typu float je podle očekávání mnohem (řádově) rychlejší. A u typu mpf (opět podle očekávání) rychlost závisí na požadované přesnosti a rozsahu. Čím větší máme požadavky na přesnost a rozsah, tím pomalejší budou výpočty, a oproti float se mohou lišit i o několik řádů(!).
Repositář s demonstračními příklady
Demonstrační příklady, s nimiž jsme se dnes seznámili a které jsou určeny pro Python 3.11 (a libovolnou vyšší verzi Pythonu) a knihovnu mpmath, jsou dostupné, jak je zvykem, na GitHubu. V tabulce níže jsou uvedeny odkazy na jednotlivé zdrojové kódy:
Odkazy na Internetu
- Číselné hodnoty s neomezeným rozsahem a přesností v programovacím jazyku Go (1)
https://www.root.cz/clanky/ciselne-hodnoty-s-neomezenym-rozsahem-a-presnosti-v-programovacim-jazyku-go-1/ - Číselné hodnoty s neomezeným rozsahem a přesností v programovacím jazyku Go (2)
https://www.root.cz/clanky/ciselne-hodnoty-s-neomezenym-rozsahem-a-presnosti-v-programovacim-jazyku-go-2/ - Číselné hodnoty s neomezeným rozsahem a přesností v programovacím jazyku Go: typ Decimal
https://www.root.cz/clanky/ciselne-hodnoty-s-neomezenym-rozsahem-a-presnosti-v-programovacim-jazyku-go-typ-decimal/ - Balíček big pro jazyk Go
https://pkg.go.dev/math/big - Zdrojové kódu pro balíček big
https://cs.opensource.google/go/go/+/master:src/math/big/ - Arbitrary-precision arithmetic
https://en.wikipedia.org/wiki/Arbitrary-precision_arithmetic - Floating-point error mitigation
https://en.wikipedia.org/wiki/Floating-point_error_mitigation - Beating Floating Point at its Own Game: Posit Arithmetic
http://www.johngustafson.net/pdfs/BeatingFloatingPoint.pdf - Unum (number format)
https://en.wikipedia.org/wiki/Unum_(number_format) - The GNU MPFR Library
https://www.mpfr.org/ - GMP: Arithmetic without limitations
https://gmplib.org/ - GNU MP 6.2.1 manual
https://gmplib.org/manual/index - Anatomy of a posit number
https://www.johndcook.com/blog/2018/04/11/anatomy-of-a-posit-number/ - Better floating point: posits in plain language
http://loyc.net/2019/unum-posits.html - Posits, a New Kind of Number, Improves the Math of AI: The first posit-based processor core gave a ten-thousandfold accuracy boost
https://spectrum.ieee.org/floating-point-numbers-posits-processor - Posit Standard Document (2022)
https://posithub.org/khub_widget - Standard for Posit™ Arithmetic (2022)
https://posithub.org/docs/posit_standard-2.pdf - Posit Calculator
https://posithub.org/widget/calculator/ - SoftPosit
https://gitlab.com/cerlane/SoftPosit - PySigmoid
https://github.com/mightymercado/PySigmoid - sgpositpy
https://github.com/xman/sgpositpy - SoftPosit.jl
https://github.com/milankl/SoftPosit.jl - SigmoidNumbers.jl
https://github.com/MohHizzani/SigmoidNumbers.jl - How many digits can float8, float16, float32, float64, and float128 contain?
https://stackoverflow.com/questions/56514892/how-many-digits-can-float8-float16-float32-float64-and-float128-contain - 15. Floating Point Arithmetic: Issues and Limitations (Python documentation)
https://docs.python.org/3/tutorial/floatingpoint.html - Number limits, overflow, and roundoff
https://www.khanacademy.org/computing/computers-and-internet/xcae6f4a7ff015e7d:digital-information/xcae6f4a7ff015e7d:limitations-of-storing-numbers/a/number-limits-overflow-and-roundoff - The upper and lower limits of IEEE-754 standard
https://math.stackexchange.com/questions/2607697/the-upper-and-lower-limits-of-ieee-754-standard - Norma IEEE 754 a příbuzní: formáty plovoucí řádové tečky
https://www.root.cz/clanky/norma-ieee-754-a-pribuzni-formaty-plovouci-radove-tecky/ - IEEE-754 Floating-Point Conversion
http://babbage.cs.qc.cuny.edu/IEEE-754.old/32bit.html - Small Float Formats
https://www.khronos.org/opengl/wiki/Small_Float_Formats - Art of Assembly language programming: The 80×87 Floating Point Coprocessors
https://courses.engr.illinois.edu/ece390/books/artofasm/CH14/CH14–3.html - Art of Assembly language programming: The FPU Instruction Set
https://courses.engr.illinois.edu/ece390/books/artofasm/CH14/CH14–4.html - INTEL 80387 PROGRAMMER'S REFERENCE MANUAL
http://www.ragestorm.net/downloads/387intel.txt - Floating-Point Formats
http://www.quadibloc.com/comp/cp0201.htm - Representation of numbers
https://mpmath.org/doc/current/technical.html - Does Python not follow IEEE-754 in case of division by zero?
https://stackoverflow.com/questions/78064239/does-python-not-follow-ieee-754-in-case-of-division-by-zero - gmpy2 2.3.0
https://pypi.org/project/gmpy2/
ČLÁNKY DO MAILU