Režim fast math v překladačích: přednosti, zápory a možné pasti

Obsah

1. Režim fast math v překladačích: přednosti, zápory a možné pasti

2. Podmínky, které musí překladač ošetřit v režimu IEEE 754 („slow math“)

3. Operace, které může překladač provádět v režimu fast math

4. Nastavení proměnné errno v průběhu výpočtů

5. Výsledky získané při použití režimu „slow math“ i fast math

6. Práce s hodnotami NaN (Not a Number)

7. Výsledky získané při použití režimu „slow math“ i fast math

8. Test na speciální hodnoty (NaN atd.)

9. Výsledek překladu demonstračního příkladu do mezikódu

10. A jak vypadá výsledný strojový kód?

11. Výpočet sumy prvků pole

12. Úprava předchozího příkladu: omezení pseudonáhodných hodnot

13. Vliv přeskládání prvků na výsledek může být enormní

14. Jak se vlastně změnil mezikód výpočtu sumy?

15. A co výsledný strojový kód?

16. Přesné nebo rychlé výpočty trigonometrických funkcí

17. Jaká funkce sinf se vlastně volá?

18. Příloha: Makefile pro překlad všech demonstračních příkladů

19. Repositář s demonstračními příklady

20. Odkazy na Internetu

1. Režim fast math v překladačích: přednosti, zápory a možné pasti

Většina moderních překladačů, ať již se jedná o překladače jazyka C, C++, nebo například i Javy, při zpracování hodnot s plovoucí řádovou čárkou dodržují pravidla a podmínky specifikované v normě IEEE 754, a to včetně dalších rozšíření této normy (IEEE 754–2008). V této normě jsou specifikovány nejenom vlastní formáty uložení numerických hodnot v systému pohyblivé řádové čárky (FP formátu), ale i pravidla implementace základních aritmetických operací s těmito hodnotami, aplikace zaokrouhlovacích režimů, způsoby některých konverzí apod. Konkrétně je v této normě popsáno:

1. Základní (basic) a rozšířený (extended) formát uložení numerických hodnot.
2. Způsob provádění základních matematických operací:
   • součet
   • rozdíl
   • součin
   • podíl
   • zbytek po dělení
   • druhá odmocnina
   • porovnání
3. Režimy zaokrouhlování.
4. Způsob práce s denormalizovanými hodnotami.
5. Pravidla konverze mezi celočíselnými formáty (integer bez a se znaménkem) a formáty s plovoucí řádovou čárkou.
6. Způsob konverze mezi různými formáty s plovoucí řádovou čárkou (single → double atd.).
7. Způsob konverze základního formátu s plovoucí řádovou čárkou na řetězec číslic (včetně nekonečen a nečíselných hodnot).
8. Práce s hodnotami NaN (not a number) a výjimkami, které mohou při výpočtech za určitých předpokladů vzniknout.

IEEE 754 je zaměřena na provádění stabilních a opakovatelných výpočtů nezávislých na použité architektuře. Ovšem některé její požadavky mohou být v praxi zbytečně „pedantické“, protože existuje poměrně velké množství aplikací, ve kterých je pouze nutné provádět aritmetické a jiné výpočty a neřešit například nekonečna, kladnou a zápornou nulu či hodnotu NaN. Příkladem může být zpracování rastrových obrázků pomocí konvolučních fitrů. Zde nikdy nebudou vstupem nekonečna ani NaN, nedojde k dělení nulou a výsledkem bude vždy běžné reálné číslo (tj. nikoli kladné či záporné nekonečno). A podobných případů existuje větší množství. Z tohoto důvodu překladače nabízejí optimalizace, které bývají shrnuty pod pojmem fast math. Těmito optimalizacemi se budeme zabývat v dnešním článku.

2. Podmínky, které musí překladač ošetřit v režimu IEEE 754 („slow math“)

Překladače programovacího jazyka C většinou při výpočtech s hodnotami s plovoucí řádovou čárkou dodržují normu IEEE 754 a popř. i její rozšíření. Co to ovšem znamená v praxi? U výpočtů je nutné detekovat různé stavy či podmínky a reagovat na ně přesně takovým způsobem, který je v IEEE 754 předepsán. Uveďme si některé případy těchto stavů nebo podmínek:

1. Operace s hodnotami s plovoucí řádovou čárkou obecněnejsou asociativní, tj. (a + b) + c ≠ a + (b + c). Překladač by tedy (pokud má odpovídat IEEE 754) neměl měnit pořadí operací atd. To, co se může zdát být maličkost, ovšem znamená praktickou nemožnost vektorizovaných výpočtů (SIMD).
2. Poměrně přesně jsou definovány výsledky operací v případě, že jedním operandem je kladné nebo záporné nekonečno. Překladač programovacího jazyka C by měl v těchto případech splňovat požadavky IEEE 754.
3. Výsledky některých výpočtů buď nejsou definovány (například podíl 0/0 nebo rozdíl nekonečen) nebo je není možné reprezentovat reálným číslem (odmocnina ze záporné hodnoty). V takových případech norma IEEE 754 předepisuje, že se má vrátit hodnota NaN (Not a Number).
4. Taktéž je nutné korektně reagovat tehdy, pokud výsledek překročí maximální rozsah hodnot. Například pro formát s dvojitou přesností (double) by se při překročení hodnoty 10³⁰⁸ mělo vrátit nekonečno a pro hodnotu –10³⁰⁸ naopak záporné nekonečno.
5. Měla by se taktéž rozlišovat kladná a záporná nulová hodnota. To se sice může zdát jako maličkost, ovšem při dělení zápornou nulou bude výsledkem záporné nekonečno, tj. hodnota zcela odlišná od kladného nekonečna.
6. Norma IEEE 754 taktéž předepisuje několik zaokrouhlovacích režimů. Ty je možné v C konfigurovat; ovšem touto problematikou se dnes nebudeme zabývat.

Všechny stavy popsané v těchto bodech jsou splněny za předpokladu, že se NEpoužije režim fast math nebo nějaká podmnožina tohoto režimu. Co je důležité: IEEE 754 předepisuje všechny vlastnosti z toho důvodu, aby byly výpočty snadno reprodukovatelné, a to i při použití odlišného překladače nebo dokonce při spuštění výpočtů na zcela odlišné architektuře. Za tyto vlastnosti (reprodukovatelnost, stabilita) ovšem mnohdy platíme delšími časy výpočtu; s využitím fast math lze mnohdy dosáhnout i řádového urychlení.

3. Operace, které může překladač provádět v režimu fast math

V případě, že je režim fast math povolen, nebo je povolena jen některá jeho část, může překladač generovat takový kód, který porušuje všechny výše uvedené vlastnosti. Všechny dostupné optimalizační metody jsou pro přehled vypsány v následující tabulce:

4. Nastavení proměnné errno v průběhu výpočtů

Některé matematické operace s hodnotami s plovoucí řádovou čárkou v případě chyby (neplatná vstupní hodnota atd.) nastavují proměnnou errno. Tato proměnná je vytvořena zvlášť pro každé vlákno, není tedy globální:

/* The error code set by various library functions.  */
extern int *__errno_location (void) __THROW __attribute_const__;
# define errno (*__errno_location ())

Hodnotu této proměnné lze později (tj. po provedení výpočtů) přečíst a zjistit, zda předchozí výpočet probíhal v pořádku. To je sice obecně velmi dobrá vlastnost, ovšem ne často se v praxi využívá, takže je možné její nastavování zakázat, a to buď již známým přepínačem -ffast-math (ten provádí více optimalizací) nebo přepínačem -fno-math-errno.

Příkladem funkcí, která proměnnou errno dokážou nastavovat, jsou funkce pro výpočet druhé odmocniny. Viz (zde zkrácený) pohled na jejich manuálovou stránku:

NAME
       sqrt, sqrtf, sqrtl - square root function
 
LIBRARY
       Math library (libm, -lm)
 
SYNOPSIS
       #include <;math.h>
 
       double sqrt(double x);
       float sqrtf(float x);
       long double sqrtl(long double x);

Zde také nalezneme informaci o errno:

ERRORS
       See math_error(7) for information on how to determine whether an error has occurred when calling these functions.
 
       The following errors can occur:
 
       Domain error: x less than -0
              errno is set to EDOM.  An invalid floating-point exception (FE_INVALID) is raised.

V dalším demonstračním příkladu si ověříme chování funkce sqrtf při výpočtu druhé odmocniny z kladného čísla a poté i čísla záporného. Program nejdříve vypíše hodnotu konstanty EDOM (viz man stránku) a poté u každého výsledku vypisuje hodnotu proměnné errno:

#include <math.h>
#include <errno.h>
#include <stdio.h>
 
int main() {
    printf("errno code for 'EDOM: Mathematics argument out of domain of function: %d\n", EDOM);
 
    float x = 2;
    printf("x=%f\terrno: %d\n", x, errno);
 
    float y = sqrtf(x);
    printf("y=%f\terrno: %d\n", y, errno);
 
    float z = sqrtf(-x);
    printf("z=%f\terrno: %d\n", z, errno);
    return 0;
}

Očekávané výsledky (za předpokladu, že nepoužijeme fast math):

5. Výsledky získané při použití režimu „slow math“ i fast math

V dalším kroku si vyzkoušíme, zda a jak vůbec bude proměnná errno nastavena v případě, že se snažíme o výpočet druhé odmocniny ze záporné vstupní hodnoty. Výpočet bude proveden celkem šestkrát, protože provedeme překlad jak s využitím Clangu (LLVM), tak i GCC. A překládat budeme postupně bez optimalizací (-O0), se zapnutými optimalizacemi na rychlost (-O3) a taktéž s překladem s optimalizacemi na rychlost a současně i s rychlými matematickými výpočty:

$ gcc -O0 -lm no_math_errno.c -o no_math_errno_gcc_O0
$ gcc -O3 -lm no_math_errno.c -o no_math_errno_gcc_O3
$ gcc -O3 -ffast-math -lm no_math_errno.c -o no_math_errno_gcc_fast_math
$ clang -O0 -lm no_math_errno.c -o no_math_errno_clang_O0
$ clang -O3 -lm no_math_errno.c -o no_math_errno_clang_O3
$ clang -O3 -ffast-math -lm no_math_errno.c -o no_math_errno_clang_fast_math

Výsledky výpočtů:

$ ./no_math_errno_clang_O0
errno code for 'EDOM: Mathematics argument out of domain of function: 33
x=2.000000      errno: 0
y=1.414214      errno: 0
z=-nan  errno: 33
 
$ ./no_math_errno_clang_O3
errno code for 'EDOM: Mathematics argument out of domain of function: 33
x=2.000000      errno: 0
y=1.414214      errno: 0
z=-nan  errno: 33
 
$ ./no_math_errno_clang_fast_math
errno code for 'EDOM: Mathematics argument out of domain of function: 33
x=2.000000      errno: 0
y=1.414214      errno: 0
z=-nan  errno: 0
 
$ ./no_math_errno_gcc_O0
errno code for 'EDOM: Mathematics argument out of domain of function: 33
x=2.000000      errno: 0
y=1.414214      errno: 0
z=-nan  errno: 33
 
$ ./no_math_errno_gcc_O3
errno code for 'EDOM: Mathematics argument out of domain of function: 33
x=2.000000      errno: 0
y=1.414214      errno: 0
z=-nan  errno: 33
 
$ ./no_math_errno_gcc_fast_math
errno code for 'EDOM: Mathematics argument out of domain of function: 33
x=2.000000      errno: 0
y=1.414214      errno: 0
z=-nan  errno: 0

Z těchto výsledků vyplývá, že oba překladače Clang i GCC proměnnou errno nastavují korektně, a to i při zapnutí optimalizací, ovšem ve chvíli, kdy se použije přepínač -ffast-math, je proměnná ignorována (což v důsledku může znamenat mnohem kratší vykonávání výpočtů).

6. Práce s hodnotami NaN (Not a Number)

Norma IEEE 754 (a popřípadě i její dodatky) poměrně přesně specifikuje, ve kterých případech je výsledkem nějaké operace hodnota NaN neboli Not a Number. Jedná se například o tyto operace:

1. Podíl 0/0 nebo libovolná kombinace –0/0, 0/-0 atd.
2. Součin nuly a nekonečna (a pochopitelně i další kombinace s různými znaménky)
3. Výpočet zbytku po dělení, pokud je dělencem nekonečno
4. Výpočet zbytku po dělení, pokud je dělitelem nula
5. Součet kladného a záporného nekonečna
6. Rozdíl dvou nekonečen
7. Sinus nekonečna (a další podobné operace)
8. Odmocnina ze záporného čísla
9. Logaritmus záporného čísla

Otestujme si, zda například rozdíl dvou (kladných) nekonečen skutečně povede k tomu, že výsledkem bude hodnota NaN. V příkladu nejdříve vypočítáme nekonečno jako podíl 1/0 a posléze se pokusíme od sebe obě nekonečna odečíst:

#include <errno.h>
#include <math.h>
#include <stdio.h>
 
int main() {
    float x = 1.0;
    float y = x / 0;
    float z = y - y;
 
    printf("x=%f\n", x);
    printf("y=%f\n", y);
    printf("z=%f\n", z);
    return 0;
}

7. Výsledky získané při použití režimu „slow math“ i fast math

Pokusme se nyní zdrojový kód demonstračního příkladu uvedeného v šesté kapitole přeložit bez povolení optimalizací, s povolením běžných optimalizací i s povolením „rychlých“ výpočtů. Zajímavé je, že překladač GCC při překladu vypisuje varování o tom, že se bude provádět dělení nulou, nicméně překlad se provede:

$ gcc -O0 -lm no_nan.c -o no_nan_gcc_O0
no_nan.c: In function ‘main’:
no_nan.c:7:17: warning: division by zero [-Wdiv-by-zero]
    7 |     float y = x / 0;
      |                 ^
 
$ gcc -O3 -lm no_nan.c -o no_nan_gcc_O3
no_nan.c: In function ‘main’:
no_nan.c:7:17: warning: division by zero [-Wdiv-by-zero]
    7 |     float y = x / 0;
      |                 ^
 
$ gcc -O3 -ffast-math -lm no_nan.c -o no_nan_gcc_fast_math
no_nan.c: In function ‘main’:
no_nan.c:7:17: warning: division by zero [-Wdiv-by-zero]
    7 |     float y = x / 0;
      |                 ^

Překladač Clang v případě, že nepovolíme přísnější režim překladu, tato varování nevypíše a překlad provede „potichu“:

$ clang -O0 -lm no_nan.c -o no_nan_clang_O0
$ clang -O3 -lm no_nan.c -o no_nan_clang_O3
$ clang -O3 -ffast-math -lm no_nan.c -o no_nan_clang_fast_math

Ovšem v tomto případě jsou mnohem důležitější vypočtené výsledky, které jednotlivé přeložené varianty vypíšou:

$ ./no_nan_gcc_O0
x=1.000000
y=inf
z=-nan
 
$ ./no_nan_gcc_O3
x=1.000000
y=inf
z=-nan
 
$ ./no_nan_gcc_fast_math
x=1.000000
y=inf
z=0.000000
 
$ ./no_nan_clang_O0
x=1.000000
y=inf
z=-nan
 
$ ./no_nan_clang_O3
x=1.000000
y=inf
z=nan
 
$ ./no_nan_clang_fast_math
x=1.000000
y=inf
z=nan

Nyní se výsledky získané GCC a Clangem odlišují, což je zajímavé. Opět si je můžeme shrnout do tabulky:

Z výsledků je zřejmé, že i výpočty s nekonečny se mohou odlišovat podle toho, zda je povolený režim fast math či nikoli.

8. Test na speciální hodnoty (NaN atd.)

U hodnoty NaN (Not a Number) ještě na chvíli zůstaneme. V následujícím (zdánlivě triviálním) demonstračním příkladu se ve funkci nazvané nan_test testuje, zda je hodnota parametru f rovna NaN či nikoli. Tento test zajišťuje standardní knihovní funkce nazvaná isnan:

DESCRIPTION
       Floating point numbers can have special values, such as infinite or NaN.
       With the macro fpclassify(x) you can find out what type x is.  The macro
       takes any floating-point expression as argument.  The result is  one  of
       the following values:
 
       FP_NAN        x is "Not a Number".
       ...
       ...
       ...
 
isnan(x)      returns a nonzero value if (fpclassify(x) == FP_NAN)

Zdrojový kód tohoto příkladu je krátký a zdánlivě bezproblémový a zcela deterministický:

#include <math.h>
 
int nan_test(float f) {
    return isnan(f);
}

9. Výsledek překladu demonstračního příkladu do mezikódu

V dalším kroku necháme demonstrační příklad z předchozí kapitoly přeložit do mezikódu LLVM IR, aby bylo patrné, jaké instrukce (na úrovni LLVM IR) budou použity. A pochopitelně bude překlad opět proveden se zákazem optimalizací, podruhé s povolením běžných optimalizací a potřetí s povolením rychlých matematických výpočtů:

$ clang -O0 -S -emit-llvm nan_test.c -o nan_test_O0.ll
 
$ clang -O3 -S -emit-llvm nan_test.c -o nan_test_O3.ll
 
$ clang -O3 -ffast-math -S -emit-llvm nan_test.c -o nan_test_fast_math.ll
nan_test.c:4:12: warning: use of NaN is undefined behavior due to the currently enabled floating-point options [-Wnan-infinity-disabled]
    4 |     return isnan(f);
      |            ^~~~~~~~
/usr/include/math.h:980:20: note: expanded from macro 'isnan'
  980 | #  define isnan(x) __builtin_isnan (x)
      |                    ^~~~~~~~~~~~~~~~~~~
1 warning generated.

Překlad se zákazem optimalizací volá funkci llvm.is.fpclass.f32:

; Function Attrs: noinline nounwind optnone uwtable
define dso_local i32 @nan_test(float noundef %0) #0 {
  %2 = alloca float, align 4
  store float %0, ptr %2, align 4
  %3 = load float, ptr %2, align 4
  %4 = call i1 @llvm.is.fpclass.f32(float %3, i32 3)
  %5 = zext i1 %4 to i32
  ret i32 %5
}

Při povolení bezpečných optimalizací se přímo používá instrukce nazvaná fcmp určená pro porovnání dvou hodnot s plovoucí řádovou čárkou. Používá se porovnání specifikované slovem uno neboli unordered (either nans), tedy test, zda je jeden z operandů NaN (a druhý operand je nulový, ten NaN zcela jistě není):

; Function Attrs: mustprogress nofree norecurse nosync nounwind willreturn memory(none) uwtable
define dso_local range(i32 0, 2) i32 @nan_test(float noundef %0) local_unnamed_addr #0 {
  %2 = fcmp uno float %0, 0.000000e+00
  %3 = zext i1 %2 to i32
  ret i32 %3
}

Ovšem v režimu fast-math je situace zcela jiná, protože se vrací konstantní hodnota 0 nebo false:

; Function Attrs: mustprogress nofree norecurse nosync nounwind willreturn memory(none) uwtable
define dso_local noundef range(i32 0, 2) i32 @nan_test(float noundef nofpclass(nan inf) %0) local_unnamed_addr #0 {
  ret i32 0
}

10. A jak vypadá výsledný strojový kód?

Nejpřesnější a možná i nejvíce překvapivé bude prozkoumání výsledného strojového kódu funkce nan_test. Na rozdíl od předchozí kapitoly nyní pro překlad použijeme GCC, aby bylo patrné, že se oba překladače chovají podobně. Céčkovskou funkci nan_test si necháme přeložit do assembleru poprvé bez optimalizací, podruhé s optimalizacemi na rychlost a potřetí s optimalizacemi i s „rychlými“ výpočty:

$ gcc -O0 -S -masm=intel nan_test.c -o nan_test_O0.asm
$ gcc -O3 -S -masm=intel nan_test.c -o nan_test_O3.asm
$ gcc -O3 -ffast-math -S -masm=intel nan_test.c -o nan_test_fast_math.asm

Výsledek překladu bez optimalizací ukazuje, že se volá instrukce ucomiss porovnávající dvě FP hodnoty (viz též popis této instrukce):

nan_test:
        push    rbp
        mov     rbp, rsp
        movss   DWORD PTR [rbp-4], xmm0
        movss   xmm0, DWORD PTR [rbp-4]
        ucomiss xmm0, DWORD PTR [rbp-4]
        setp    al
        movzx   eax, al
        pop     rbp
        ret

Po zapnutí optimalizací se provede porovnání pomocí instrukce ucomiss, která mj. nastavuje i příznak parity flag. Pokud jsou hodnoty neporovnatelné (unordered), což je případ porovnání dvou NaN, je tento příznak nastaven na jedničku, takže následuje instrukce setp, která na základě tohoto příznaku změní návratovou hodnotu funkce:

nan_test:
        xor     eax, eax
        ucomiss xmm0, xmm0
        setp    al
        ret

Ovšem při zapnutí -ffast-math je funkce přeložena takovým způsobem, že vrací konstantní nulu, nezávisle na předaném parametru! To je zcela zřejmé při pohledu na vygenerovaný kód této funkce:

nan_test:
        xor     eax, eax
        ret

11. Výpočet sumy prvků pole

Zajímavé situace nastávají u programového kódu, ve kterém se FP operace provádí ve smyčce nad prvky pole či polí. Takové operace je teoreticky možné „vektorizovat“, tj. realizovat je SIMD operacemi. Ovšem má to jeden háček – vektorizace typicky předpokládá asociativitu operací, což je ovšem vlastnost, kterou norma IEEE 754 (zcela správně) omezuje. Pokud tedy budeme chtít, aby se například výpočet sumy prvků pole prováděl zrychleným způsobem (součet například osmi prvků jedinou instrukcí), je nutné povolit režim fast-math či alespoň -fassociative-math. Ověřme si, jak se tato optimalizační metoda projeví ve výsledcích výpočtu.

První příklad byl převzat z článku GCC -ffast-math Explained: What It Actually Does, Key Details & Practical Examples. Provádí se v něm výpočet sumy všech prvků pole naplněných pseudonáhodnými hodnotami (vždy stejnými). Vypisuje se jak výsledek výpočtů, tak i jejich doba:

#include <stdlib.h>
#include <stdio.h>
#include <time.h>
 
#define N 1000000
 
float array[N];
 
float sum_array() {
    float sum = 0.0f;
    for (int i = 0; i < N; i++) {
        sum += array[i]; // Simple addition
    }
    return sum;
}
 
int main() {
    // Initialize array with random values (0.0f to 1.0f)
    for (int i = 0; i < N; i++) {
        array[i] = (float)rand() / RAND_MAX;
    }
 
    // Time the sum
    clock_t start = clock();
    float   result = sum_array();
    clock_t end = clock();
 
    printf("Sum: %.4f\n", result);
    printf("Time: %.4f ms\n", (double)(end - start) / CLOCKS_PER_SEC * 1000);
    return 0;
}

Překlad tohoto příkladu provedeme oběma překladači a s již známými přepínači:

$ gcc -O0 -lm array_sum_1.c -o array_sum_1_gcc_O0
$ gcc -O3 -lm array_sum_1.c -o array_sum_1_gcc_O3
$ gcc -O3 -ffast-math -lm array_sum_1.c -o array_sum_1_gcc_fast_math
$ clang -O0 -lm array_sum_1.c -o array_sum_1_clang_O0
$ clang -O3 -lm array_sum_1.c -o array_sum_1_clang_O3
$ clang -O3 -ffast-math -lm array_sum_1.c -o array_sum_1_clang_fast_math

Výsledky výpočtů i s jejich časy:

$ ./array_sum_1_clang_O0
Sum: 500003.7500
Time: 4.5270 ms
 
$ ./array_sum_1_clang_O3
Sum: 500003.7500
Time: 1.0720 ms
 
$ ./array_sum_1_clang_fast_math
Sum: 500004.9688
Time: 0.2700 ms
 
$ ./array_sum_1_gcc_O0
Sum: 500003.7500
Time: 4.3110 ms
 
$ ./array_sum_1_gcc_O3
Sum: 500003.7500
Time: 1.2460 ms
 
$ ./array_sum_1_gcc_fast_math
Sum: 500007.9688
Time: 0.3490 ms

Z výsledků plyne, že režim fast math vede skutečně k mnohem rychlejším výpočtům (cca čtyřikrát rychlejším oproti „běžné“ optimalizované variantě), ovšem výsledky se budou odlišovat. Rozdíl je sice v tomto případě relativně malý, ale zřetelný.

DESCRIPTION
       The  rand()  function  returns  a  pseudo-random  integer  in  the  range  0 to RAND_MAX inclusive (i.e., the mathematical range
       [0, RAND_MAX]).
 
       The srand() function sets its argument as the seed for a new sequence of pseudo-random integers to be returned by rand().  These
       sequences are repeatable by calling srand() with the same seed value.
 
       If no seed value is provided, the rand() function is automatically seeded with a value of 1.
       ...
       ...
       ...

12. Úprava předchozího příkladu: omezení pseudonáhodných hodnot

Použití pseudonáhodných, kterými je pole inicializováno, ve skutečnosti demonstrační příklad zbytečně komplikuje. Pokusme se tedy jeho zdrojový kód upravit tak, aby se například sčítala nějaká aritmetická řada, a to ideálně s takovými hodnotami, aby výsledkem nebylo nekonečno:

#include <stdio.h>
#include <time.h>
 
#define N 10000000
 
float array[N];
 
float sum_array() {
    float sum = 0.0f;
    for (int i = 0; i < N; i++) {
        sum += array[i];
    }
    return sum;
}
 
int main() {
    for (int i = 0; i < N; i++) {
        array[i] = i/10.0;
    }
 
    clock_t start = clock();
    float   result = sum_array();
    clock_t end = clock();
 
    printf("Sum: %.0f\n", result);
    printf("Time: %.4f ms\n", (double)(end - start) / CLOCKS_PER_SEC * 1000);
    return 0;
}

Překlad:

$ gcc -O0 -lm array_sum_2.c -o array_sum_2_gcc_O0
$ gcc -O3 -lm array_sum_2.c -o array_sum_2_gcc_O3
$ gcc -O3 -ffast-math -lm array_sum_2.c -o array_sum_2_gcc_fast_math
$ clang -O0 -lm array_sum_2.c -o array_sum_2_clang_O0
$ clang -O3 -lm array_sum_2.c -o array_sum_2_clang_O3
$ clang -O3 -ffast-math -lm array_sum_2.c -o array_sum_2_clang_fast_math

Výsledky výpočtů i s jejich časy:

$ ./array_sum_2_clang_O0
Sum: 5081496813568
Time: 36.8590 ms
 
$ ./array_sum_2_clang_O3
Sum: 5081496813568
Time: 11.0620 ms
 
$ ./array_sum_2_clang_fast_math
Sum: 5000135704576
Time: 3.1710 ms
 
$ ./array_sum_2_gcc_O0
Sum: 5081496813568
Time: 39.7770 ms
 
$ ./array_sum_2_gcc_O3
Sum: 5081496813568
Time: 11.1300 ms
 
$ ./array_sum_2_gcc_fast_math
Sum: 4988595077120
Time: 4.0580 ms

13. Vliv přeskládání prvků na výsledek může být enormní

Mohlo by se zdát, že využití asociativity operací (zde konkrétně operace součtu) má jen relativně malý vliv na vypočtené výsledky. Ukažme si ovšem, co se stane ve chvíli, kdy je pole, jehož prvky sčítáme, naplněno sekvencí hodnot, přičemž sousední prvky mají opačná znaménka:

    for (int i = 0; i < N; i++) {
        array[i] = i % 2 ? i : -i;
    }

Úplný zdrojový kód takto upraveného demonstračního příkladu vypadá následovně:

#include <stdio.h>
#include <time.h>
 
#define N 100000
 
float array[N];
 
float sum_array() {
    float sum = 0.0f;
    for (int i = 0; i < N; i++) {
        sum += array[i];
    }
    return sum;
}
 
int main() {
    for (int i = 0; i < N; i++) {
        array[i] = i % 2 ? i : -i;
    }
 
    clock_t start = clock();
    float   result = sum_array();
    clock_t end = clock();
 
    printf("Sum: %.4f\n", result);
    printf("Time: %.4f ms\n", (double)(end - start) / CLOCKS_PER_SEC * 1000);
    return 0;
}

Výsledky jsou nyní v případě fast math značně odlišné od očekávané hodnoty 50000 (tedy N/2):

$ ./array_sum_3_clang_O0
Sum: 50000.0000
Time: 0.9900 ms
 
$ ./array_sum_3_clang_O3
Sum: 50000.0000
Time: 0.2540 ms
 
$ ./array_sum_3_clang_fast_math
Sum: 76800.0000
Time: 0.0570 ms
 
$ ./array_sum_3_gcc_O0
Sum: 50000.0000
Time: 1.0070 ms
 
$ ./array_sum_3_gcc_O3
Sum: 50000.0000
Time: 0.2030 ms
 
$ ./array_sum_3_gcc_fast_math
Sum: 88064.0000
Time: 0.0690 ms

14. Jak se vlastně změnil mezikód?

A jak vypadá mezikód LLVM IR získaný překladem demonstračního příkladu z předchozí kapitoly? Neoptimalizovaná verze nepočítá s asociativitou operací:

; Function Attrs: noinline nounwind optnone uwtable
define dso_local float @sum_array() #0 {
  %1 = alloca float, align 4
  %2 = alloca i32, align 4
  store float 0.000000e+00, ptr %1, align 4
  store i32 0, ptr %2, align 4
  br label %3
 
3:                                                ; preds = %13, %0
  %4 = load i32, ptr %2, align 4
  %5 = icmp slt i32 %4, 100000
  br i1 %5, label %6, label %16
 
6:                                                ; preds = %3
  %7 = load i32, ptr %2, align 4
  %8 = sext i32 %7 to i64
  %9 = getelementptr inbounds [100000 x float], ptr @array, i64 0, i64 %8
  %10 = load float, ptr %9, align 4
  %11 = load float, ptr %1, align 4
  %12 = fadd float %11, %10
  store float %12, ptr %1, align 4
  br label %13
 
13:                                               ; preds = %6
  %14 = load i32, ptr %2, align 4
  %15 = add nsw i32 %14, 1
  store i32 %15, ptr %2, align 4
  br label %3, !llvm.loop !4
 
16:                                               ; preds = %3
  %17 = load float, ptr %1, align 4
  ret float %17
}

Výpočet uvnitř těla smyčky vypadá takto:

  %10 = load float, ptr %9, align 4
  %11 = load float, ptr %1, align 4
  %12 = fadd float %11, %10
  store float %12, ptr %1, align 4

Překlad se zapnutím optimalizací, ale nikoli v režimu fast math:

; Function Attrs: nofree norecurse nosync nounwind memory(read, argmem: none, inaccessiblemem: none) uwtable
define dso_local float @sum_array() local_unnamed_addr #0 {
  br label %2
 
1:                                                ; preds = %2
  ret float %35
 
2:                                                ; preds = %2, %0
  %3 = phi i64 [ 0, %0 ], [ %36, %2 ]
  %4 = phi float [ 0.000000e+00, %0 ], [ %35, %2 ]
  %5 = getelementptr inbounds nuw [100000 x float], ptr @array, i64 0, i64 %3
  %6 = load float, ptr %5, align 16, !tbaa !3
  %7 = fadd float %4, %6
  %8 = or disjoint i64 %3, 1
  %9 = getelementptr inbounds nuw [100000 x float], ptr @array, i64 0, i64 %8
  %10 = load float, ptr %9, align 4, !tbaa !3
  %11 = fadd float %7, %10
  %12 = or disjoint i64 %3, 2
  %13 = getelementptr inbounds nuw [100000 x float], ptr @array, i64 0, i64 %12
  %14 = load float, ptr %13, align 8, !tbaa !3
  %15 = fadd float %11, %14
  %16 = or disjoint i64 %3, 3
  %17 = getelementptr inbounds nuw [100000 x float], ptr @array, i64 0, i64 %16
  %18 = load float, ptr %17, align 4, !tbaa !3
  %19 = fadd float %15, %18
  %20 = or disjoint i64 %3, 4
  %21 = getelementptr inbounds nuw [100000 x float], ptr @array, i64 0, i64 %20
  %22 = load float, ptr %21, align 16, !tbaa !3
  %23 = fadd float %19, %22
  %24 = or disjoint i64 %3, 5
  %25 = getelementptr inbounds nuw [100000 x float], ptr @array, i64 0, i64 %24
  %26 = load float, ptr %25, align 4, !tbaa !3
  %27 = fadd float %23, %26
  %28 = or disjoint i64 %3, 6
  %29 = getelementptr inbounds nuw [100000 x float], ptr @array, i64 0, i64 %28
  %30 = load float, ptr %29, align 8, !tbaa !3
  %31 = fadd float %27, %30
  %32 = or disjoint i64 %3, 7
  %33 = getelementptr inbounds nuw [100000 x float], ptr @array, i64 0, i64 %32
  %34 = load float, ptr %33, align 4, !tbaa !3
  %35 = fadd float %31, %34
  %36 = add nuw nsw i64 %3, 8
  %37 = icmp eq i64 %36, 100000
  br i1 %37, label %1, label %2, !llvm.loop !7
}

Zde došlo k rozbalení smyčky, ale stále se provádí sekvence:

  %26 = load float, ptr %25, align 4, !tbaa !3
  %27 = fadd float %23, %26
  %29 = getelementptr inbounds nuw [100000 x float], ptr @array, i64 0, i64 %28
  %30 = load float, ptr %29, align 8, !tbaa !3
  %31 = fadd float %27, %30

Ovšem při zapnutí režimu fast-math je vše jinak:

; Function Attrs: nofree norecurse nosync nounwind memory(read, argmem: none, inaccessiblemem: none) uwtable
define dso_local nofpclass(nan inf) float @sum_array() local_unnamed_addr #0 {
  br label %1
 
1:                                                ; preds = %1, %0
  %2 = phi i64 [ 0, %0 ], [ %39, %1 ]
  %3 = phi <4 x float> [ zeroinitializer, %0 ], [ %37, %1 ]
  %4 = phi <4 x float> [ zeroinitializer, %0 ], [ %38, %1 ]
  %5 = getelementptr inbounds nuw [100000 x float], ptr @array, i64 0, i64 %2
  %6 = getelementptr inbounds nuw i8, ptr %5, i64 16
  %7 = load <4 x float>, ptr %5, align 16, !tbaa !3
  %8 = load <4 x float>, ptr %6, align 16, !tbaa !3
  %9 = fadd fast <4 x float> %7, %3
  %10 = fadd fast <4 x float> %8, %4
  %11 = add nuw nsw i64 %2, 8
  %12 = getelementptr inbounds nuw [100000 x float], ptr @array, i64 0, i64 %11
  %13 = getelementptr inbounds nuw i8, ptr %12, i64 16
  %14 = load <4 x float>, ptr %12, align 16, !tbaa !3
  %15 = load <4 x float>, ptr %13, align 16, !tbaa !3
  %16 = fadd fast <4 x float> %14, %9
  %17 = fadd fast <4 x float> %15, %10
  %18 = add nuw nsw i64 %2, 16
  %19 = getelementptr inbounds nuw [100000 x float], ptr @array, i64 0, i64 %18
  %20 = getelementptr inbounds nuw i8, ptr %19, i64 16
  %21 = load <4 x float>, ptr %19, align 16, !tbaa !3
  %22 = load <4 x float>, ptr %20, align 16, !tbaa !3
  %23 = fadd fast <4 x float> %21, %16
  %24 = fadd fast <4 x float> %22, %17
  %25 = add nuw nsw i64 %2, 24
  %26 = getelementptr inbounds nuw [100000 x float], ptr @array, i64 0, i64 %25
  %27 = getelementptr inbounds nuw i8, ptr %26, i64 16
  %28 = load <4 x float>, ptr %26, align 16, !tbaa !3
  %29 = load <4 x float>, ptr %27, align 16, !tbaa !3
  %30 = fadd fast <4 x float> %28, %23
  %31 = fadd fast <4 x float> %29, %24
  %32 = add nuw nsw i64 %2, 32
  %33 = getelementptr inbounds nuw [100000 x float], ptr @array, i64 0, i64 %32
  %34 = getelementptr inbounds nuw i8, ptr %33, i64 16
  %35 = load <4 x float>, ptr %33, align 16, !tbaa !3
  %36 = load <4 x float>, ptr %34, align 16, !tbaa !3
  %37 = fadd fast <4 x float> %35, %30
  %38 = fadd fast <4 x float> %36, %31
  %39 = add nuw nsw i64 %2, 40
  %40 = icmp eq i64 %39, 100000
  br i1 %40, label %41, label %1, !llvm.loop !7
 
41:                                               ; preds = %1
  %42 = fadd fast <4 x float> %38, %37
  %43 = tail call fast float @llvm.vector.reduce.fadd.v4f32(float 0.000000e+00, <4 x float> %42)
  ret float %43
}

Můžeme zde vidět jak rozbalení smyčky, tak i vektorizované součty:

  %30 = fadd fast <4 x float> %28, %23

Ovšem kvůli vektorizaci se prvky sčítají v odlišném pořadí, než je předepsáno ve zdrojovém kódu.

15. A co výsledný strojový kód?

Pro úplnost si ještě ukážeme výsledný strojový kód, který překladače dokážou vygenerovat. V tomto případě použiji GCC (Clang byl využit pro LLVM IR).

Varianta přeložená bez optimalizací:

sum_array:
        push    rbp
        mov     rbp, rsp
        pxor    xmm0, xmm0
        movss   DWORD PTR [rbp-4], xmm0
        mov     DWORD PTR [rbp-8], 0
        jmp     .L2
.L3:
        mov     eax, DWORD PTR [rbp-8]
        cdqe
        movss   xmm0, DWORD PTR array[0+rax*4]
        movss   xmm1, DWORD PTR [rbp-4]
        addss   xmm0, xmm1
        movss   DWORD PTR [rbp-4], xmm0
        add     DWORD PTR [rbp-8], 1
.L2:
        cmp     DWORD PTR [rbp-8], 99999
        jle     .L3
        movss   xmm0, DWORD PTR [rbp-4]
        pop     rbp
        ret

Zde se opakovaně provádí součet mezivýsledku s prvkem přečteným z pole:

        movss   xmm0, DWORD PTR array[0+rax*4]
        movss   xmm1, DWORD PTR [rbp-4]
        addss   xmm0, xmm1
        movss   DWORD PTR [rbp-4], xmm0

Kód získaný po zapnutí optimalizací, ovšem bez fast math:

sum_array:
        mov     eax, OFFSET FLAT:array
        pxor    xmm0, xmm0
.L2:
        addss   xmm0, DWORD PTR [rax]
        add     rax, 16
        addss   xmm0, DWORD PTR [rax-12]
        addss   xmm0, DWORD PTR [rax-8]
        addss   xmm0, DWORD PTR [rax-4]
        cmp     rax, OFFSET FLAT:array+400000
        jne     .L2
        ret

Programová smyčka je sice optimalizována a rozbalena, ovšem stále se provádí součty ve správném pořadí. Když si odmyslíme zvýšení offsetu RAX, provádí se čtyři součty v jedné iteraci:

        addss   xmm0, DWORD PTR [rax]
        addss   xmm0, DWORD PTR [rax+4]
        addss   xmm0, DWORD PTR [rax+8]
        addss   xmm0, DWORD PTR [rax+12]

A konečně výsledek překladu v režimu fast math:

sum_array:
        mov     eax, OFFSET FLAT:array
        mov     edx, OFFSET FLAT:array+400000
        pxor    xmm0, xmm0
.L2:
        addps   xmm0, XMMWORD PTR [rax]
        add     rax, 32
        addps   xmm0, XMMWORD PTR [rax-16]
        cmp     rdx, rax
        jne     .L2
        movaps  xmm1, xmm0
        movhlps xmm1, xmm0
        addps   xmm1, xmm0
        movaps  xmm0, xmm1
        shufps  xmm0, xmm1, 85
        addps   xmm0, xmm1
        ret

Povšimněte si vektorizované (SIMD) instrukce addps, která navíc ukazuje, jak se prvky sčítají (prvek N+prvek N+8 atd.). Předpokládá se zde asociativita operace součtu!

16. Přesné nebo rychlé výpočty trigonometrických funkcí

V závěrečné části dnešního článku se ještě zmiňme o další optimalizaci (na rychlost), kterou překladače provádí v případě povolení optimalizačního režimu fast math. Dochází resp. může docházet k náhradě přesných, ale pomalých trigonometrických funkcí za jejich rychlejší varianty. Tyto rychlé varianty mohou být založeny například na výpočtu jen několika členů Taylorova rozvoje, na pomocných tabulkách s mezivýsledky atd.

V dalším demonstračním příkladu se sčítá milion hodnot trigonometrické funkce sinus, přičemž vstupní hodnoty jsou v intervalu 0 až 2π. Teoreticky by součet hodnot měl vyjít nulový, ovšem v tomto případě je nutné počítat s mnoha chybami (zaukrouhlení, omezená přesnost atd.):

#include <stdio.h>
#include <math.h>
#include <time.h>
 
#define N 1000000
 
int main() {
    float result = 0.0f;
    clock_t start = clock();
 
    for (int i = 0; i < N; i++) {
        result += sinf(M_PI*2.0*i/N);
    }
 
    clock_t end = clock();
 
    printf("Sum:  %.4f\n", result);
    printf("Time: %.4f ms\n", (double)(end - start) / CLOCKS_PER_SEC * 1000);
    return 0;
}

Jak je již zvykem, provedeme překlad oběma překladači jazyka C a to pro všechny testované nastavení optimalizací:

$ gcc -O0 -lm sin.c -o sin_gcc_O0
$ gcc -O3 -lm sin.c -o sin_gcc_O3
$ gcc -O3 -ffast-math -lm sin.c -o sin_gcc_fast_math
$ clang -O0 -lm sin.c -o sin_clang_O0
$ clang -O3 -lm sin.c -o sin_clang_O3
$ clang -O3 -ffast-math -lm sin.c -o sin_clang_fast_math

Výsledky, včetně časů výpočtu, vypadají takto:

$ ./sin_clang_O0
Sum:  -0.0212
Time: 45.0080 ms
 
$ ./sin_clang_O3
Sum:  -0.0212
Time: 13.4320 ms
 
$ ./sin_clang_fast_math
Sum:  0.0006
Time: 11.0540 ms
 
$ ./sin_gcc_O0
Sum:  -0.0212
Time: 17.2820 ms
 
$ ./sin_gcc_O3
Sum:  -0.0212
Time: 14.1390 ms
 
$ ./sin_gcc_fast_math
Sum:  0.0111
Time: 4.8700 ms

Podobný příklad, ovšem postavený na funkci sin a datovém typu double, je taktéž součástí příkladů, takže si můžete otestovat i jeho chování.

17. Jaká funkce sinf se vlastně volá?

Po překladu zdrojového kódu, ve kterém se volá funkce sinf, do assembleru můžeme snadno zjistit, jaká varianta funkce sinf se vlastně volá. Pokud při překladu neuvedeme přepínač -ffast-math, bude se volat funkce ze standardní knihovny. To si můžeme snadno ověřit pohledem do assembleru vygenerovaného překladačem:

main:
        push    rbp
        push    rbx
        mov     ebx, 1
        sub     rsp, 24
        call    clock
        mov     DWORD PTR [rsp+12], 0x00000000
        mov     rbp, rax
.L2:
        pxor    xmm0, xmm0
        cvtsi2sd        xmm0, ebx
        mulsd   xmm0, QWORD PTR .LC1[rip]
        add     ebx, 1
        divsd   xmm0, QWORD PTR .LC2[rip]
        cvtsd2ss        xmm0, xmm0
        call    sinf
        addss   xmm0, DWORD PTR [rsp+12]
        movss   DWORD PTR [rsp+12], xmm0
        cmp     ebx, 1000000
        jne     .L2
        call    clock
        mov     edi, OFFSET FLAT:.LC3
        pxor    xmm0, xmm0
        mov     rbx, rax
        mov     eax, 1
        cvtss2sd        xmm0, DWORD PTR [rsp+12]
        call    printf
        sub     rbx, rbp
        pxor    xmm0, xmm0
        mov     edi, OFFSET FLAT:.LC5
        cvtsi2sd        xmm0, rbx
        divsd   xmm0, QWORD PTR .LC2[rip]
        mov     eax, 1
        mulsd   xmm0, QWORD PTR .LC4[rip]
        call    printf
        add     rsp, 24
        xor     eax, eax
        pop     rbx
        pop     rbp
        ret

Ovšem v případě, že je při překladu použit přepínač -ffast-math bude se v případě překladače GCC volat funkce, která bude mít vygenerované unikátní jméno, zde konkrétně _ZGVbN4v_sinf:

main:
        push    rbp
        push    rbx
        xor     ebx, ebx
        sub     rsp, 56
        call    clock
        pxor    xmm6, xmm6
        movdqa  xmm2, XMMWORD PTR .LC0[rip]
        mov     rbp, rax
        mov     eax, 4
        movaps  XMMWORD PTR [rsp], xmm6
        movd    xmm6, eax
        pshufd  xmm7, xmm6, 0
        movaps  XMMWORD PTR [rsp+32], xmm7
.L2:
        pshufd  xmm1, xmm2, 238
        cvtdq2pd        xmm0, xmm2
        movaps  XMMWORD PTR [rsp+16], xmm2
        add     ebx, 1
        mulpd   xmm0, XMMWORD PTR .LC5[rip]
        cvtdq2pd        xmm1, xmm1
        mulpd   xmm1, XMMWORD PTR .LC5[rip]
        cvtpd2ps        xmm0, xmm0
        cvtpd2ps        xmm1, xmm1
        movlhps xmm0, xmm1
        call    _ZGVbN4v_sinf
        addps   xmm0, XMMWORD PTR [rsp]
        movdqa  xmm2, XMMWORD PTR [rsp+16]
        paddd   xmm2, XMMWORD PTR [rsp+32]
        movaps  XMMWORD PTR [rsp], xmm0
        cmp     ebx, 249999
        jne     .L2
        call    clock
        movaps  xmm5, XMMWORD PTR [rsp]
        mov     edi, OFFSET FLAT:.LC2
        mov     rbx, rax
        mov     eax, 1
        movaps  xmm1, xmm5
        sub     rbx, rbp
        movhlps xmm1, xmm5
        addps   xmm1, xmm5
        movaps  xmm0, xmm1
        shufps  xmm0, xmm1, 85
        addps   xmm0, xmm1
        subss   xmm0, DWORD PTR .LC1[rip]
        cvtss2sd        xmm0, xmm0
        call    printf
        pxor    xmm0, xmm0
        mov     edi, OFFSET FLAT:.LC4
        cvtsi2sd        xmm0, rbx
        mulsd   xmm0, QWORD PTR .LC3[rip]
        mov     eax, 1
        call    printf
        add     rsp, 56
        xor     eax, eax
        pop     rbx
        pop     rbp
        ret

18. Příloha: Makefile pro překlad všech demonstračních příkladů

Všechny demonstrační příklady využívající překladač Clang a GCC, které byly použity v dnešním článku, je možné přeložit s využitím souboru Makefile, jehož obsah je vypsán pod tímto odstavcem:

outputs := \
           no_math_errno_gcc_O0 no_math_errno_gcc_O3 no_math_errno_gcc_fast_math \
           no_math_errno_clang_O0 no_math_errno_clang_O3 no_math_errno_clang_fast_math \
           no_nan_gcc_O0 no_nan_gcc_O3 no_nan_gcc_fast_math \
           no_nan_clang_O0 no_nan_clang_O3 no_nan_clang_fast_math \
           nan_test_O0.ll nan_test_O3.ll nan_test_fast_math.ll \
           nan_test_O0.asm nan_test_O3.asm nan_test_fast_math.asm \
           sin_gcc_O0 sin_gcc_O3 sin_gcc_fast_math \
           sin_clang_O0 sin_clang_O3 sin_clang_fast_math \
           sin_fast.asm sin_slow.asm \
           array_sum_3_O0.ll array_sum_3_O3.ll array_sum_3_fast_math.ll \
           array_sum_3_O0.asm array_sum_3_O3.asm array_sum_3_fast_math.asm \
           array_sum_1_gcc_O0 array_sum_1_gcc_O3 array_sum_1_gcc_fast_math \
           array_sum_1_clang_O0 array_sum_1_clang_O3 array_sum_1_clang_fast_math \
           array_sum_2_gcc_O0 array_sum_2_gcc_O3 array_sum_2_gcc_fast_math \
           array_sum_2_clang_O0 array_sum_2_clang_O3 array_sum_2_clang_fast_math \
           array_sum_3_gcc_O0 array_sum_3_gcc_O3 array_sum_3_gcc_fast_math \
           array_sum_3_clang_O0 array_sum_3_clang_O3 array_sum_3_clang_fast_math
 
all:    $(outputs)
 
clean:
        rm -f $(outputs)
 
no_math_errno_gcc_O0:   no_math_errno.c
        gcc -O0 -lm $< -o $@
 
no_math_errno_gcc_O3:   no_math_errno.c
        gcc -O3 -lm $< -o $@
 
no_math_errno_gcc_fast_math:    no_math_errno.c
        gcc -O3 -ffast-math -lm $< -o $@
 
no_math_errno_clang_O0: no_math_errno.c
        clang -O0 -lm $< -o $@
 
no_math_errno_clang_O3: no_math_errno.c
        clang -O3 -lm $< -o $@
 
no_math_errno_clang_fast_math:  no_math_errno.c
        clang -O3 -ffast-math -lm $< -o $@
 
no_nan_gcc_O0:  no_nan.c
        gcc -O0 -lm $< -o $@
 
no_nan_gcc_O3:  no_nan.c
        gcc -O3 -lm $< -o $@
 
no_nan_gcc_fast_math:   no_nan.c
        gcc -O3 -ffast-math -lm $< -o $@
 
no_nan_clang_O0:        no_nan.c
        clang -O0 -lm $< -o $@
 
no_nan_clang_O3:        no_nan.c
        clang -O3 -lm $< -o $@
 
no_nan_clang_fast_math: no_nan.c
        clang -O3 -ffast-math -lm $< -o $@
 
nan_test_O0.ll: nan_test.c
        clang -O0 -S -emit-llvm $< -o $@
 
nan_test_O3.ll: nan_test.c
        clang -O3 -S -emit-llvm $< -o $@
 
nan_test_fast_math.ll:  nan_test.c
        clang -O3 -ffast-math -S -emit-llvm $< -o $@
 
nan_test_O0.asm:        nan_test.c
        gcc -O0 -S -masm=intel $< -o $@
 
nan_test_O3.asm:        nan_test.c
        gcc -O3 -S -masm=intel $< -o $@
 
nan_test_fast_math.asm: nan_test.c
        gcc -O3 -ffast-math -S -masm=intel $< -o $@
 
sin_gcc_O0:     sin.c
        gcc -O0 -lm $< -o $@
 
sin_gcc_O3:     sin.c
        gcc -O3 -lm $< -o $@
 
sin_gcc_fast_math:      sin.c
        gcc -O3 -ffast-math $< -o $@ -lm
 
sin_clang_O0:   sin.c
        clang -O0 -lm $< -o $@
 
sin_clang_O3:   sin.c
        clang -O3 -lm $< -o $@
 
sin_clang_fast_math:    sin.c
        clang -O3 -ffast-math -lm $< -o $@
 
sin_slow.asm:   sin.c
        gcc -S -masm=intel -O3 $< -o $@ -lm
 
sin_fast.asm:   sin.c
        gcc -S -masm=intel -O3 -ffast-math $< -o $@ -lm
 
array_sum_3_O0.ll:      array_sum_3.c
        clang -O0 -S -emit-llvm $< -o $@
 
array_sum_3_O3.ll:      array_sum_3.c
        clang -O3 -S -emit-llvm $< -o $@
 
array_sum_3_fast_math.ll:       array_sum_3.c
        clang -O3 -ffast-math -S -emit-llvm $< -o $@
 
array_sum_3_O0.asm:     array_sum_3.c
        gcc -O0 -S -masm=intel $< -o $@
 
array_sum_3_O3.asm:     array_sum_3.c
        gcc -O3 -S -masm=intel $< -o $@
 
array_sum_3_fast_math.asm:      array_sum_3.c
        gcc -O3 -ffast-math -S -masm=intel $< -o $@
 
array_sum_1_gcc_O0:     array_sum_1.c
        gcc -O0 -lm $< -o $@
 
array_sum_1_gcc_O3:     array_sum_1.c
        gcc -O3 -lm $< -o $@
 
array_sum_1_gcc_fast_math:      array_sum_1.c
        gcc -O3 -ffast-math -lm $< -o $@
 
array_sum_1_clang_O0:   array_sum_1.c
        clang -O0 -lm $< -o $@
 
array_sum_1_clang_O3:   array_sum_1.c
        clang -O3 -lm $< -o $@
 
array_sum_1_clang_fast_math:    array_sum_1.c
        clang -O3 -ffast-math -lm $< -o $@
 
array_sum_2_gcc_O0:     array_sum_2.c
        gcc -O0 -lm $< -o $@
 
array_sum_2_gcc_O3:     array_sum_2.c
        gcc -O3 -lm $< -o $@
 
array_sum_2_gcc_fast_math:      array_sum_2.c
        gcc -O3 -ffast-math -lm $< -o $@
 
array_sum_2_clang_O0:   array_sum_2.c
        clang -O0 -lm $< -o $@
 
array_sum_2_clang_O3:   array_sum_2.c
        clang -O3 -lm $< -o $@
 
array_sum_2_clang_fast_math:    array_sum_2.c
        clang -O3 -ffast-math -lm $< -o $@
 
array_sum_3_gcc_O0:     array_sum_3.c
        gcc -O0 -lm $< -o $@
 
array_sum_3_gcc_O3:     array_sum_3.c
        gcc -O3 -lm $< -o $@
 
array_sum_3_gcc_fast_math:      array_sum_3.c
        gcc -O3 -ffast-math -lm $< -o $@
 
array_sum_3_clang_O0:   array_sum_3.c
        clang -O0 -lm $< -o $@
 
array_sum_3_clang_O3:   array_sum_3.c
        clang -O3 -lm $< -o $@
 
array_sum_3_clang_fast_math:    array_sum_3.c
        clang -O3 -ffast-math -lm $< -o $@

19. Repositář s demonstračními příklady

Demonstrační příklady, s nimiž jsme se v dnešním článku seznámili a které jsou určeny pro překlad s využitím Clangu nebo GCC, jsou dostupné, jak je v tomto seriálu zvykem, na GitHubu:

20. Odkazy na Internetu

1. Clang: -ffast-math
   https://clang.llvm.org/doc­s/UsersManual.html#cmdopti­on-ffast-math
2. Gcc: Options That Control Optimization
   https://gcc.gnu.org/online­docs/gcc/Optimize-Options.html
3. GCC -ffast-math Explained: What It Actually Does, Key Details & Practical Examples
   https://www.codestudy.net/blog/what-does-gcc-s-ffast-math-actually-do/
4. What does gcc's ffast-math actually do?
   https://stackoverflow.com/qu­estions/7420665/what-does-gccs-ffast-math-actually-do
5. Beware of fast-math
   https://simonbyrne.github­.io/notes/fastmath/
6. Understanding fast-math
   https://www.nutrient.io/blog/un­derstanding-fast-math/
7. „Fast math“ optimization (Wikipedia)
   https://en.wikipedia.org/wiki/Floating-point_arithmetic#%22Fast_mat­h%22_optimization
8. 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/
9. IEEE-754 Floating-Point Conversion
   http://babbage.cs.qc.cuny.edu/IEEE-754.old/32bit.html
10. Small Float Formats
    https://www.khronos.org/o­pengl/wiki/Small_Float_For­mats
11. Binary-coded decimal
    https://en.wikipedia.org/wiki/Binary-coded_decimal
12. Why Intel is betting on BFLOAT16 to be a game changer for deep learning training? Hint: Range trumps Precision
    https://hub.packtpub.com/why-intel-is-betting-on-bfloat16-to-be-a-game-changer-for-deep-learning-training-hint-range-trumps-precision/
13. half-rs (pro Rust)
    https://github.com/starkat99/half-rs
14. float16 (pro Go)
    https://github.com/x448/float16
15. bfloat16 – Hardware Numerics Definition
    https://software.intel.com/en-us/download/bfloat16-hardware-numerics-definition
16. Intel Prepares To Graft Google’s Bfloat16 Onto Processors
    https://www.nextplatform.com/2019/07/15/in­tel-prepares-to-graft-googles-bfloat16-onto-processors/
17. A Study of BFLOAT16 for Deep Learning Training
    https://arxiv.org/pdf/1905.12322.pdf
18. BFloat16s.jl
    https://github.com/JuliaCom­puting/BFloat16s.jl
19. Half Precision Arithmetic: fp16 Versus bfloat16
    https://nhigham.com/2018/12/03/half-precision-arithmetic-fp16-versus-bfloat16/
20. bfloat16 floating-point format (Wikipedia)
    https://en.wikipedia.org/wi­ki/Bfloat16_floating-point_format
21. Unum (number format)
    https://en.wikipedia.org/wi­ki/Unum_(number_format)#Po­sit
22. Performance Benefits of Half Precision Floats
    https://software.intel.com/en-us/articles/performance-benefits-of-half-precision-floats
23. 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/
24. IEEE-754 Floating-Point Conversion
    http://babbage.cs.qc.cuny.edu/IEEE-754.old/32bit.html
25. Small Float Formats
    https://www.khronos.org/o­pengl/wiki/Small_Float_For­mats
26. Binary-coded decimal
    https://en.wikipedia.org/wiki/Binary-coded_decimal
27. Chen–Ho encoding
    https://en.wikipedia.org/wi­ki/Chen%E2%80%93Ho_encoding
28. Densely packed decimal
    https://en.wikipedia.org/wi­ki/Densely_packed_decimal
29. A Summary of Chen-Ho Decimal Data encoding
    http://speleotrove.com/decimal/chen-ho.html
30. Art of Assembly language programming: The 80×87 Floating Point Coprocessors
    https://courses.engr.illi­nois.edu/ece390/books/arto­fasm/CH14/CH14–3.html
31. Art of Assembly language programming: The FPU Instruction Set
    https://courses.engr.illi­nois.edu/ece390/books/arto­fasm/CH14/CH14–4.html
32. INTEL 80387 PROGRAMMER'S REFERENCE MANUAL
    http://www.ragestorm.net/dow­nloads/387intel.txt
33. Floating-Point Formats
    http://www.quadibloc.com/com­p/cp0201.htm
34. Data types (SciPy)
    https://docs.scipy.org/doc/numpy-1.13.0/user/basics.types.html
35. Is there a way to use errno safely in a multi-threaded application?
    https://stackoverflow.com/qu­estions/449778/is-there-a-way-to-use-errno-safely-in-a-multi-threaded-application
36. FCOMI/FCOMIP/FUCOMI/FUCOMIP — Compare Floating-Point Values and Set EFLAGS
    https://www.felixcloutier­.com/x86/fcomi:fcomip:fuco­mi:fucomip
37. Operations generating NaN
    https://en.wikipedia.org/wi­ki/NaN#Operations_generatin­g_NaN
38. Taylorova řada
    https://cs.wikipedia.org/wi­ki/Taylorova_%C5%99ada