Hlavní navigace

Aritmetické operace s hodnotami ve formátu plovoucí řádové čárky

7. 6. 2006
Doba čtení: 13 minut

Sdílet

Ve třetí části seriálu bude popsán způsob provádění základních aritmetických operací s formáty plovoucí čárky podle normy IEEE 754. Také si ukážeme implementaci výpočtů složitějších funkcí, například druhé odmocniny, goniometrických funkcí apod. Nezapomeneme ani na velmi významný algoritmus CORDIC, který byl použit i v některých matematických koprocesorech.

Obsah

1. Operace prováděné s hodnotami v FP reprezentaci podle normy IEEE 754
2. Součet a rozdíl dvou FP hodnot
3. Součin dvou FP hodnot
4. Podíl dvou FP hodnot
5. Iterační výpočet druhé odmocniny
6. Demonstrační příklad na výpočet druhé odmocniny
7. Goniometrické funkce počítané pomocí číselných řad
8. Úvodní informace o iteračním algoritmu CORDIC
9. Obsah dalšího pokračování tohoto seriálu

1. Operace prováděné s hodnotami v FP reprezentaci podle normy IEEE 754

V předchozích dvou pokračováních tohoto seriálu jsme si podrobně popsali formáty uložení číselných hodnot definovaných normou IEEE 754. Jednalo se o formát jednoduché přesnosti (single), dvojité přesnosti (double) a konečně o rozšířený formát (temporary/ex­tended), který je implementován zejména v matematických koprocesorech firmy Intel řady i87 a samozřejmě také v kompatibilních koprocesorech (AMD apod.). Dnešní část začneme popisem provádění základních aritmetických operací, tj. operací sčítání a odčítání (kapitola druhá), násobení (třetí kapitola) a dělení (čtvrtá kapitola). Po popisu těchto operací si v páté a šesté kapitole řekneme, jakým způsobem je možné implementovat výpočet druhé odmocniny kladných reálných čísel. V sedmé kapitole bude popsán princip výpočtů goniometrických funkcí pomocí součtu členů číselných řad a nakonec si v osmé kapitole přiblížíme principy fungování algoritmu CORDIC, pomocí nějž je možné zjednodušit, urychlit a/nebo zpřesnit výpočty složitějších matematických funkcí, zejména funkcí goniometrických.

2. Součet a rozdíl dvou FP hodnot

Při sečítání dvou hodnot uložených ve formátu pohyblivé řádové binární čárky podle normy IEEE 754 (může se přitom jednat o libovolný formát popsaný v předchozí části tohoto seriálu, zejména však o formát jednoduché a dvojité přesnosti) by se mělo postupovat podle následujícího algoritmu:

  1. Nejprve je zjištěno, zda se sčítání či odečítání neprovádí se speciálními hodnotami, tj. s hodnotou NaN (not a number) či s nekonečny. Výsledek součtu či rozdílu provedeného se speciálními hodnotami je uveden v první tabulce (pod popisem algoritmu).
  2. Poté je provedeno nalezení většího z obou sčítaných čísel. Porovnávají se však pouze hodnoty exponentů, nikoli mantis (zde se ukazuje výhodnost použití exponentů posunutých o bias, neboť při jejich porovnávání nemusíme brát v úvahu znaménko). Pokud mají obě sčítaná či odečítaná čísla stejné exponenty, přeskočí se následující bod.
  3. Dále se vypočítá číselný rozdíl exponentů, tj. diff=e1-e2. Mantisu menšího čísla je nutné posunout doprava o tolik bitů, kolik činí rozdíl exponentů (samozřejmě po vydělení rozdílu hodnotou 2n, ve skutečnosti se pouze porovnávají exponenty bez zbytečného převodu na jejich reálnou hodnotu). Bitový posun mantisy odpovídá násobení či dělení dané hodnoty mocninou čísla 2, stejně tak zvýšení či snížení hodnoty exponentu. Pro jiné báze (radixy) to však neplatí.
  4. Podle znaménkových bitů obou operandů je vyčíslen buď součet nebo rozdíl jejich mantis. Zde se jedná o „klasický“ součet pomocí vícebitové sčítačky s postupným či zrychleným přenosem. Při odečítání je nutné vzít v úvahu, že se nejedná o čísla uložená ve dvojkovém doplňku, tj. odečítání je poněkud složitější než v případě příště popsaných numerických hodnot uložených v FX formátu.
  5. Pokud hodnota výsledku přeteče (tj. výsledek se nevejde do počtu bitů mantisy), posunou se bity v mantise doprava o jeden bit a exponent se zvýší o jedničku. Z toho vyplývá, že součet je nutné provádět na sčítačce minimálně o jeden bit širší, než je bitová šířka mantisy (postačí přenos CARRY).
  6. V případě, že mantisa nemá nastaven nejvyšší bit na jedničku, posouvá se (i několikrát) doleva a exponent se přitom snižuje. To je takzvaná normalizace, která zajistí, že nejvyšší bit mantisy je vždy jedničkový a tudíž nemusí být v mantise uložen – to platí zejména pro formáty single a double, u formátu temporary je nejvyšší bit mantisy uložen v bitu číslo 63.
  7. Dále se musí vyjádřit znaménko výsledku: v případě FPU sčítání se jednoduše zkopíruje znaménko většího čísla, u odečítání se musí vzít v úvahu i pozice čísla v operaci, tj. zda se jedná o menšenec či menšitel.
  8. Po vyčíslení znaménka se provede zaokrouhlení hodnoty podle právě nastaveného zaokrouhlovacího režimu. V normě jsou specifikovány čtyři způsoby zaokrouhlení (tzv. rounding modes): zaokrouhlení směrem k nule, zaokrouhlení směrem ke kladnému nekonečnu, zaokrouhlení směrem k zápornému nekonečnu a konečně zaokrouhlení k nejbližšímu sudému reprezentovatelnému číslu (pojmenování posledního zaokrouhlovacího režimu není přesné, protože rozhodnutí, zda se jedná o sudé číslo, je zapotřebí provést pouze pro hraniční hodnoty).
  9. Pokud během výpočtu dojde k tomu, že exponent dosáhne své maximální hodnoty, dojde k přetečení výsledku. To již není možné nijak obejít, pouze zvýšením přesnosti a/nebo rozsahu, tj. volbou jiného FP formátu. V každém případě je při přetečení jako výsledek operace sčítání či odečítání vráceno kladné nebo záporné nekonečno, podle samostatně vypočítaného znaménka.
  10. Pokud naopak hodnota exponentu dosáhne své minimální hodnoty, dochází k podtečení. Při něm se buď generuje výjimka, nebo se hodnota výsledku nastaví na kladnou či zápornou nulu, opět podle samostatně vypočítaného znaménka výsledku.

Kvůli přetékání, podtékání a zaokrouhlování hodnot během sečítání/odečítání přestávají platit některá matematická pravidla, například A+B-A==B či (A+B)+C=A+(B+C) a naopak začínají platit pravidla nová, třeba při reálném programování velmi nebezpečné a záludné A+B=A (platí samozřejmě pouze pro některé hodnoty A a B). V prvním bodu algoritmu pro součet či rozdíl dvou FP hodnot je zjišťováno, zda se výpočet neprovádí se speciálními hodnotami. Výsledek operace v případě použití speciálních hodnot je ukázán v následující tabulce:

Hodnota X Hodnota Y Výsledek operace X+Y
konečná +∞ +∞
konečná -∞ -∞
+∞ konečná +∞
-∞ konečná -∞
+∞ +∞ +∞
-∞ -∞ -∞
+∞ -∞ NaN
-∞ +∞ NaN
NaN libovolná NaN
libovolná NaN NaN

3. Součin dvou FP hodnot

Součin dvou numerických hodnot uložených ve formátu plovoucí řádové čárky je opět – vzhledem k rozdělení hodnot na mantisu a exponent – komplikovanější než prosté poslání obsahů dvou registrů do násobičky. Nejprve si teoreticky ukažme, jakým způsobem se vynásobí dvě hodnoty reprezentované svou mantisou a exponentem. Vstupními hodnotami jsou čísla x1 a x2:

x1=(-1)s1 × 2e1 × m1
x2=(-1)s2 × 2e2 × m2

Součin těchto hodnot můžeme vyjádřit vzorcem:

x1x2=(-1)s1×(-1)s2 × 2e1+e2 × m­1×m2

Prakticky to znamená, že (pokud nebudeme uvažovat nutnou normalizaci a zaokrouhlení výsledků) se mantisy obou hodnot navzájem vynásobí a exponenty se sečtou. Nejjednodušší je operace se znaménkovými bity – na ně se aplikuje bitová operace nonekvivalence (XOR:⊕). Vzhledem k tomu, že se musíme v co největší míře vyhnout přetečení, podtečení a denormalizaci hodnot, je násobení dvou FP hodnot složitější a provádí se následujícím postupem:

  1. Podobně jako u operace sečítání či odečítání, i při násobení je nejprve zjištěno, zda do operace nevstupují speciální hodnoty. S těmi je zacházeno tak, jak ukazuje druhá tabulka (umístěná pod popisem algoritmu).
  2. Exponent výsledku se rovná součtu obou exponentů, tj. provede se operace e=e1+e2. Od tohoto výsledku je nutné odečíst posun exponentů (bias), a to z toho důvodu, že v mantise jsou uloženy hodnoty s binární čárkou umístěnou za prvním bitem a prostým vynásobením mantis by vlastně vznikl výsledek posunutý o několik (binárních) řádů směrem doleva.
  3. Výsledná mantisa vznikne součinem obou mantis, tj. provede se operace m=m1×m2.
  4. Pokud výsledná mantisa přeteče, je bitově posunuta doprava o n bitů a exponent je o tuto hodnotu n zvýšen.
  5. Když nemá mantisa nejvyšší bit roven 1, posouvá se doleva a exponent se naopak snižuje (normalizace). Tato operace se může provádět i několikrát.
  6. Pokud se po součtu exponentů překročí maximální hodnota (dojde k přetečení), je výsledkem nekonečno (kladné či záporné, to záleží na znaménkových bitech operandů).
  7. Pokud se po součtu exponentů překročí minimální hodnota (dojde k podtečení), je výsledkem nula (kladná či záporná, opět záleží na znaménkových bitech obou operandů)
  8. Výsledný znaménkový bit je vyjádřen nonekvivalencí obou znaménkových bitů, tj. s=s1 ⊕ s2.
  9. Po vyčíslení znaménka se provede zaokrouhlení hodnoty podle právě nastaveného zaokrouhlovacího režimu.
Hodnota X Hodnota Y Výsledek operace X×Y
nenulová kladná +∞ +∞
nenulová kladná -∞ -∞
nenulová záporná +∞ -∞
nenulová záporná -∞ +∞
nula (0,0) +∞ NaN
nula (0,0) -∞ NaN
+∞ nenulová kladná +∞
-∞ nenulová kladná -∞
+∞ nenulová záporná -∞
-∞ nenulová záporná +∞
+∞ nula (0,0) NaN
-∞ nula (0,0) NaN
+∞ +∞ +∞
+∞ -∞ -∞
-∞ +∞ -∞
-∞ -∞ +∞
NaN libovolná NaN
libovolná NaN NaN

4. Podíl dvou FP hodnot

Při výpočtu podílu dvou FP hodnot se využívá následující vztah:

x1/x2=(-1)s1×(-1)s2 × 2e1-e2 × m1/m2

Vlastní vydělení se může provádět buď jako samostatná operace (tak to činí dnešní FPU), která vyžaduje děličku, nebo je – v případě nutnosti šetření logickými obvody – možné použít postupné odečítání dělitele od dělence (algoritmus známý ze základních škol). Třetí možností, kterou implementují některé FPU a GPU, je výpočet převrácené hodnoty dělitele pomocí specializovaných postupů a následné vynásobení dělence převrácenou hodnotou dělitele. Všechny tři postupy však musí zaručit, že se při dělení nulou bude zachovávat znaménko výsledku (tj. kladné či záporné nekonečno) a také se zajistí, aby se detekovala nepovolená operace 0/0:

x/0+=+∞
x/0-=-∞
0/0=NaN

5. Iterační výpočet druhé odmocniny

Výpočet druhé odmocniny je v praxi velmi často používaný, proto je implementován ve většině FPU či GPU. My si zde ukážeme jednu z možností implementace, která spočívá v iterativním zpřesňování odhadu výsledku. Nyní si tedy povězme, jakým způsobem se může iterativní výpočet druhé odmocniny provádět. Vstupní hodnotou algoritmu je původní číslo y, výstupem je (po i-té iteraci) odhad výsledku xi. V každé iteraci se provádí poměrně jednoduchá operace:

xi+1=1/2 (xi+y/xi)

Zbývá pouze nastavit počáteční odhad x1. Většinou se volí hodnota x1=y/2, protože je ji možné velmi jednoduše spočítat prostým bitovým posunem mantisy nebo snížením hodnoty exponentu. Přesnější (a tím pádem i rychlejší, jelikož se provede méně iterací) je však odhad vypočtený z mantisy a exponentu vstupující hodnoty y:

x1=m×2e/2

tj. hodnota exponentu se sníží na polovinu. Důkaz pro korektnost takto vytvořeného prvního odhadu je poměrně přímočarý. Vstupní hodnota y může být rozepsána na mantisu a exponent:

y=m×2e

Znaménko můžeme ignorovat, protože odmocniny se počítají pouze pro kladné hodnoty. Prvotní odhad vychází z aplikace odmocniny na výše uvedený rozepsaný vztah:

y1/2=m1/2×2e/2

Hodnotu mantisy je možné (resp. nutné) ignorovat, protože normalizovaná mantisa představuje hodnoty z intervalu <1, 2-ε> a tento interval se aplikací odmocniny nemění.

6. Demonstrační příklad na výpočet druhé odmocniny

Pro ilustraci iterativního výpočtu druhé odmocniny a zejména faktu, že řešení rychle spěje ke korektnímu výsledku, jsem vytvořil demonstrační příklad v programovacím jazyku C. Po překladu a spuštění tohoto příkladu se zobrazí tabulka s postupně zpřesňovaným výpočtem druhé odmocniny čísla 10 000. Kromě průběžných výsledků je zobrazována i absolutní a relativní chyba výpočtu. Všimněte si, že již po deseti iteracích klesá relativní chyba pod setinu procenta (alespoň v rámci dané přesnosti FP formátu), takže iterační metoda je zvolena vhodně, protože řešení konverguje velmi rychle.

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

// tato funkce provede jednu iteraci
// se zlepšeným odhadem výsledku x_i
float sqrt_step(float xi, float y)
{
    return 1.0/2.0*(xi+y/xi);
}

int main(void)
{
    // vstupní hodnota, ze které se počítá odmocnina
    double y=10000;
    // postupně zpřesňovaný odhad výsledku
    double xi=y;
    // hodnota pro porovnání výsledků
    double sqr=sqrt(y);
    // počitadlo iterací
    int i;
    for (i=0; i<20; i++) {
        xi=sqrt_step(xi, y);
        printf("%d\t%11.5f\t%11.5f\t%10.2f%%\n",
                i+1,
                xi,                       // vypočtená hodnota
                fabs(xi-sqr),             // absolutní chyba
                100.0f*fabs(xi-sqr)/sqr); // relativní chyba
    }
    return 0;
}

// finito 
iterace  odhad x_i       absolutní chyba   relativní chyba
1        5000.50000      4900.50000        4900.50%
2        2501.25000      2401.25000        2401.25%
3        1252.62402      1152.62402        1152.62%
4         630.30365       530.30365         530.30%
5         323.08450       223.08450         223.08%
6         177.01808        77.01808          77.02%
7         116.75475        16.75475          16.75%
8         101.20219         1.20219           1.20%
9         100.00714         0.00714           0.01%
10        100.00000         0.00000           0.00%
11        100.00000         0.00000           0.00%
12        100.00000         0.00000           0.00%
13        100.00000         0.00000           0.00%
14        100.00000         0.00000           0.00%
15        100.00000         0.00000           0.00%
16        100.00000         0.00000           0.00%
17        100.00000         0.00000           0.00%
18        100.00000         0.00000           0.00%
19        100.00000         0.00000           0.00%
20        100.00000         0.00000           0.00% 

7. Goniometrické funkce počítané pomocí číselných řad

Goniometrické funkce je možné počítat více způsoby. V další části tohoto seriálu si ukážeme výpočet goniometrických funkcí pomocí algoritmu CORDIC, dnes si řekneme, jakým způsobem je možné použít číselné řady. Uvažujme funkci sin(), která je, jak víme už ze střední školy, funkcí periodickou s periodou 2π (pro další goniometrické funkce platí podobné vztahy a závěry, tj. v dalším textu se jimi nebudu zabývat).

Na první pohled to vypadá, že se můžeme omezit pouze na výpočet této funkce pro vstupní hodnoty z intervalu <0, 2π>. Ve skutečnosti však můžeme využít velké symetrie této funkce, a to dokonce dvakrát. V intervalu <π, 2π> je funkce zrcadlově otočena vůči vzoru z intervalu <0, π>. Kromě toho je funkce na intervalu <0, π/2> vertikálně symetrická s intervalem <π/2, π>. Z obou těchto symetrií vyplývá, že je zapotřebí přímým výpočtem zjistit hodnoty funkce sin() pouze na intervalu <0, π/2>, hodnoty v ostatních intervalech se dopočítají buď změnou znaménka výsledku, nebo změnou vstupní hodnoty. Daný interval „kvadrantu“ má i další přednost: funkce sin() je v něm neustále rostoucí bez nežádoucích zákmitů a lokálních minim či maxim.

Všechny goniometrické funkce je možné na určitém intervalu (typicky kvadrantu) rozepsat pomocí Taylorova rozvoje. V případě funkce sin() dostaneme nekonečný rozvoj, který vypadá následovně:

sin(x)=x1/1! – x3/3! + x5/5! – x7/7! + …

Tento rozvoj sice na první pohled vypadá výpočetně velmi náročně, ve skutečnosti je však zapotřebí poměrně malého množství členů, o čemž nás přesvědčí následující demonstrační příklad, který pro výpočet funkce sin() na intervalu <0, π/2> používá pouze první tři nenulové členy posloupnosti:

#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif

#define EPSILON 10.0e-20

// tato funkce provede výpočet funkce
// sin() pomocí tří nenulových členů
// Taylorova rozvoje
double computeSin(double x)
{
    double t1=x/1.0;
    double t3=x*x*x/6.0;
    double t5=x*x*x*x*x/120.0;
    return t1-t3+t5;
}

int main(void)
{
    // počitadlo iterací
    double alfa;
    for (alfa=0.0; alfa<=M_PI/2.0; alfa+=M_PI/40.0) {
        double sin1=sin(alfa);
        double sin2=computeSin(alfa);
        printf("%5.3f\t%8.6f\t%8.6f\t%8.6f\t%6.2f%%\n",
                alfa,
                sin1,
                sin2,
                fabs(sin2-sin1),
                fabs(sin1)<EPSILON ? 0.0 : 100.0*fabs(sin2-sin1)/sin1);
    }
    return 0;
}

// finito 

Po překladu a spuštění tohoto příkladu se zobrazí následující tabulka (bez prvního řádku) s vyčíslenými hodnotami funkce sin() pro 21 úhlů z prvního kvadrantu. Ve druhém sloupci je výsledek standardní matematické operace sin(), třetí sloupec obsahuje hodnoty vypočtené pomocí tří členů Taylorova rozvoje a v posledních dvou sloupcích je spočtena a zobrazena absolutní a relativní chyba. Všimněte si, že chyba postupně stoupá od prakticky nulové hodnoty (přesný výsledek) až po cca 0,5%. To je jedna z vlastností Taylorova rozvoje; při vzdalování se od počáteční hodnoty chyba roste. Pokud by bylo zapotřebí použít vyšší přesnosti, není nic jednoduššího, než přidat další jeden či dva členy rozvoje.

úhel    sin()           Taylorův rozvoj absolutní chyba   relativní chyba
0.000   0.000000    0.000000    0.000000      0.00%
0.079   0.078459    0.078459    0.000000      0.00%
0.157   0.156434    0.156434    0.000000      0.00%
0.236   0.233445    0.233445    0.000000      0.00%
0.314   0.309017    0.309017    0.000000      0.00%
0.393   0.382683    0.382684    0.000000      0.00%
0.471   0.453990    0.453992    0.000001      0.00%
0.550   0.522499    0.522502    0.000003      0.00%
0.628   0.587785    0.587793    0.000008      0.00%
0.707   0.649448    0.649465    0.000017      0.00%
0.785   0.707107    0.707143    0.000036      0.01%
0.864   0.760406    0.760477    0.000071      0.01%
0.942   0.809017    0.809146    0.000129      0.02%
1.021   0.852640    0.852866    0.000226      0.03%
1.100   0.891007    0.891386    0.000379      0.04%
1.178   0.923880    0.924493    0.000613      0.07%
1.257   0.951057    0.952017    0.000961      0.10%
1.335   0.972370    0.973834    0.001464      0.15%
1.414   0.987688    0.989867    0.002178      0.22%
1.492   0.996917    1.000088    0.003170      0.32%
1.571   1.000000    1.004525    0.004525      0.45% 

8. Uvodní informace o iteračním algoritmu CORDIC

CORDIC neboli COordinate ROtation DIgital Computer je výpočetní metoda používající specializovaný iterativní algoritmus, který slouží v první řadě k výpočtu trigonometrických funkcí s předem známou přesností, tj. sin(), cos(), tan() atd. Po malých úpravách je možné tento algoritmus využít i pro další výpočty, například fáze (úhlu) a velikosti komplexních čísel, vyčíslení logaritmických funkcí, hyperbolických funkcí (sinh(), cosh(), tanh()) atd. Jednou ze základních vlastností metody CORDIC, která se snad nejvíce v minulosti postarala o jeho velké rozšíření, je jednoduchost operací, které se v každé iteraci provádí – používá se pouze sčítání, odečítání a bitové posuny, žádné další operace nejsou zapotřebí; dokonce není ani nutné implementovat násobičku (přesněji řečeno, u některých „odvozených“ funkcí, například tan() je zapotřebí dělička k vydělení vypočtených hodnot sin() a cos()). Díky této vlastnosti bylo možné CORDIC využít i ve velmi omezených zařízeních, například kalkulačkách či osmibitových mikrořadičích (omezením zde mám na mysli relativně malý počet logických členů). Podrobnější informace o této velmi zajímavé metodě budou uvedeny v následujícím pokračování seriálu.

UX DAy - tip 2

9. Obsah dalšího pokračování tohoto seriálu

Ve čtvrté části tohoto seriálu dokončíme tři díly věnované numerickým formátům, které používají plovoucí řádovou čárku. Podrobněji si popíšeme princip algoritmu CORDIC, který je implementován v mnoha digitálních zařízeních, včetně kalkulaček, digitálních signálových procesorů a specializovaných obvodů vytvořených například v FPGA.

Byl pro vás článek přínosný?

Autor článku

Vystudoval VUT FIT a v současné době pracuje na projektech vytvářených v jazycích Python a Go.