4.3 Lebegőpontos reprezentáció és numerikus hibák

Egy monitorozó rendszerben a késleltetés, a CPU-terhelés és az időbélyeg gyakran lebegőpontos érték. Ezek azonban nem valós számok korlátlan pontosságú tárolói: véges számú reprezentálható érték közül választanak. Emiatt két különböző időpont azonos tárolt értékké válhat, egy összegzés eredménye függhet a műveleti sorrendtől, és egy optimalizálási kapcsoló megváltoztathatja a numerikus szerződést.

A Promóciók, signedness és túlcsordulás részhez hasonlóan itt is külön kell választanunk a típus, a konverzió és a művelet szabályait. Az alap ISO C17; az IEEE 754 szerinti, kettes alapú működést külön implementációs feltételként kezeljük.

4.3.1 Típusok és ellenőrizhető numerikus profil

Az ISO C17 három alapvető valós lebegőpontos típusa a float, a double és a long double. Értékkészletük egymásba ágyazott: a float értékei a double, a double értékei a long double értékkészletének részét képezik. A long double ettől még nem feltétlenül pontosabb a double-nél.

A sizeof önmagában nem mondja meg a pontosságot: tárolási méretet mér C bájtban, amely kitöltést is tartalmazhat. A numerikus tulajdonságokat a <float.h> írja le.

A DBL_ makróknak megfelelő FLT_ és LDBL_ változataik is vannak. A subnormál támogatásról a *_HAS_SUBNORM makrók adnak információt; értékük akár meghatározhatatlan támogatást is jelezhet.

Konkrét laborprofil, nem ISO-követelmény: Linux, x86-64 System V ABI, LP64, CHAR_BIT == 8, alapértelmezett ABI-szerinti long double.

Ez nem univerzális „64 bites gépes” táblázat. Más ABI vagy fordítókapcsoló más long double-t eredményezhet.

A literál típusa is számít: 0.1 típusa double, 0.1f típusa float, 0.1L típusa long double. A korábban float-ra kerekített érték későbbi double-re alakítása nem állítja vissza az elveszett információt.

4.3.2 Reprezentáció és időadatok felbontása

Kettes alapú formátumban egy nem nulla normál szám alakja:

$$x=(-1)^s\left(1+\frac{k}{2^{p-1}}\right)2^e,$$

ahol $s$ az előjelbit, $p$ a szignifikand pontossága, $k$ egész, $0\le k<2^{p-1}$, $e$ pedig a megengedett tartományba eső kitevő.

Azonos kitevő mellett a szomszédos értékek távolsága:

$$\Delta=2^{e-(p-1)}.$$

Bizonyításvázlat. Két szomszédos reprezentációban $k$ eggyel változik. A szignifikand különbsége $2^{-(p-1)}$; ezt $2^e$-vel szorozva kapjuk a távolságot.

Binary64 esetén $p=53$. Egy körülbelül $1{,}7\cdot10^9$ másodperces Unix-időbélyegre:

$$2^{30}\le t<2^{31},$$

ezért

$$\Delta=2^{30-52}=2^{-22}\text{ s} \approx238{,}4\text{ ns}.$$

A nagy abszolút időbélyeg tehát már nem képes minden nanoszekundumos eltérést megkülönböztetni. Ha előbb két időbélyeget double-re alakítunk, majd kivonjuk őket, az elveszett részletek nem nyerhetők vissza.

Gyakorlati megoldás: az időpontokat megfelelő egész reprezentációban tartjuk, tartományellenőrzött kivonással képezzük a különbséget, és csak az időtartamot alakítjuk másodpercben megadott double-re. A kisebb nagyságrendhez finomabb abszolút felbontás tartozik.

4.3.2.1 Subnormál értékek és különleges állapotok

A subnormál számok a legkisebb normál szám és a nulla közötti tartomány egy részét fedik le, csökkenő relatív pontossággal. Binary64 esetén:

• a legkisebb pozitív normál érték $2^{-1022}$;
• a legkisebb pozitív subnormál érték $2^{-1074}$;
• a subnormál értékek lépésköze állandó.

Ez a fokozatos alulcsordulás (gradual underflow). Nem minden végrehajtási profil őrzi meg: egyes hardveres üzemmódok a subnormál eredményeket nullára cserélik.

Az IEEE 754 szerinti megvalósítás további értékeket is biztosít:

• Infinity: pozitív vagy negatív végtelen; nem azonos a legnagyobb véges számmal.
• NaN, Not a Number: például érvénytelen numerikus művelet eredményének jelölése. NaN esetén az x == x hamis, az x != x igaz; vizsgálatára mégis az isnan való.
• Signed zero: a pozitív és negatív nulla numerikusan egyenlő, de a signbit megkülönbözteti őket.

Az ISO C17 önmagában nem ígéri minden ilyen érték és minden IEEE-kivételszabály támogatását. Különösen nem szabad abból kiindulni, hogy a nullával osztás minden C implementáción egyszerűen végtelent ad.

Mérési pipeline-ban általában külön adatminőségi szabály szükséges: elfogadunk-e végtelent, vagy minden nem véges értéket hibának tekintünk? A NaN nem helyettesíti automatikusan a hiányzó adat explicit jelölését.

4.3.3 Kerekítés, hiba és összehasonlítás

A $0{,}1=1/10$ kettes számrendszerben végtelen törtként írható fel. Véges bináris reprezentáció ezért csak közelíti. Ugyanez igaz sok művelet pontos matematikai eredményére.

Legyen $x$ a pontos, $\widehat{x}$ a tárolt érték. Az abszolút és relatív hiba:

$$E_{\mathrm{abs}}=|\widehat{x}-x|, \qquad E_{\mathrm{rel}}=\frac{|\widehat{x}-x|}{|x|}, \quad x\ne0.$$

Legközelebbi értékre kerekítésnél, normál tartományban, túlcsordulás és alulcsordulás nélkül egy binary64 művelet szokásos hibamodellje:

$$\operatorname{fl}(a\circ b)=(a\circ b)(1+\delta), \qquad |\delta|\le2^{-53}.$$

Ez egyetlen megfelelően kerekített művelet hibakorlátja, nem egy teljes algoritmusé. A DBL_EPSILON ebben a profilban $2^{-52}$: sem a legkisebb pozitív érték, sem univerzális összehasonlítási küszöb.

Közel azonos mennyiségek kivonásakor a bemenetek hibája az eredményhez képest nagy lehet. Maga a kivonás akár pontos is lehet, miközben a különbség relatív hibája jelentős: a probléma a már pontatlan bemenetekből indul.

Numerikus összehasonlításra gyakori szerződés:

$$|a-b|\le \max\left(\tau_{\mathrm{abs}}, \tau_{\mathrm{rel}}\max(|a|,|b|)\right).$$

Az abszolút tolerancia a mennyiség mértékegységében értendő; a relatív tolerancia dimenzió nélküli. Nulla közelében az abszolút tag szükséges. A toleranciát a szenzor felbontása, az algoritmus hibája és az alkalmazás követelménye alapján választjuk, nem egy véletlenszerű „kis epsilonból”.

4.3.3.1 Kidolgozott C17-példa

A következő függvény nemnegatív, véges késleltetésekre készült. Hibás érték vagy tolerancia esetén hamisat ad. A $0\le\tau_{\mathrm{rel}}\le1$ feltétel és a nemnegatív bemenetek megakadályozzák a különbség, illetve a relatív küszöb túlcsordulását. Nem ad át tulajdonjogot és nem tart fenn kölcsönzött hozzáférést.

#include <float.h>
#include <math.h>
#include <stdbool.h>
#include <stdio.h>

static bool latency_close(double a, double b,
                          double atol, double rtol)
{
    if (!isfinite(a) || !isfinite(b) ||
        !isfinite(atol) || !isfinite(rtol) ||
        a < 0.0 || b < 0.0 || atol < 0.0 ||
        rtol < 0.0 || rtol > 1.0) {
        return false;
    }

    return fabs(a - b) <= fmax(atol, rtol * fmax(a, b));
}

int main(void)
{
    double sum = 0.1 + 0.2;
    double reference = 0.3;

    printf("radix=%d, double precision=%d, eval=%d\n",
           FLT_RADIX, DBL_MANT_DIG, FLT_EVAL_METHOD);
    printf("sizes: %zu %zu %zu C_byte\n",
           sizeof(float), sizeof(double), sizeof(long double));
    printf("sum=%.*g, hex=%a\n",
           DBL_DECIMAL_DIG, sum, sum);
    printf("exact=%d, tolerant=%d\n",
           sum == reference,
           latency_close(sum, reference, 1e-12, 1e-9));
    return 0;
}

Javasolt, rögzített kiinduló fordítópárok: GCC 14.2.0 és Clang 19.1.7, glibc 2.40 mellett, a fenti x86-64 célprofilon. A parancsok külön laborhoz tartoznak; nem jelentenek elvégzett fordítást vagy tesztet.

gcc -std=c17 -O2 -Wall -Wextra -Wpedantic \
    -march=x86-64 -fno-fast-math -ffp-contract=off \
    numeric.c -lm -o numeric-gcc

clang -std=c17 -O2 -Wall -Wextra -Wpedantic \
    -march=x86-64 -fno-fast-math -ffp-contract=off \
    numeric.c -lm -o numeric-clang

Binary64 és legközelebbi értékre kerekítés mellett a várt decimális összeg 0.30000000000000004; a pontos összehasonlítás hamis, a választott toleranciás összehasonlítás igaz. Ezek profilhoz kötött várt eredmények.

4.3.4 Könyvtár, környezet és optimalizálás

A <math.h> műveletei között szerepel a sqrt, log, fabs, nextafter és fma. A megfelelő f és l végződésű változatok általában float, illetve long double értékekkel dolgoznak. A fpclassify, isfinite, isnan és signbit az értékek osztályozását segítik.

Hibakezelésnél a math_errhandling jelzi az errno-alapú és/vagy lebegőpontos kivételjelzést. Nem helyes minden függvényhívás után kizárólag az errno-t figyelni.

A <fenv.h> a lebegőpontos környezetet teszi elérhetővé:

• fegetround, fesetround: kerekítési mód lekérdezése és módosítása;
• feclearexcept, fetestexcept: támogatott kivételflagek törlése és vizsgálata.

A lebegőpontos kivételflag nem feltétlenül jelent processzmegszakítást vagy signalt. A támogatott módokat, makrókat és a függvények visszatérési értékét ellenőrizni kell.

ISO C17 mechanizmus: a #pragma STDC FENV_ACCESS ON jelzi, hogy a program megfigyeli vagy módosítja a környezetet. Fordítókorlát: támogatása és optimalizálással való együttműködése verziófüggő; a -frounding-math önmagában nem bizonyít teljes támogatást.

A műveletösszevonás (contraction) szintén fontos. A külön szorzás és összeadás két kerekítést végezhet:

$$\operatorname{round}(\operatorname{round}(ab)+c),$$

míg az összevont művelet egyet:

$$\operatorname{round}(ab+c).$$

Az eredmények eltérhetnek. A #pragma STDC FP_CONTRACT OFF ennek szabványos szabályozási eszköze; tényleges támogatását ellenőrizni kell. A laborparancs ezért külön -ffp-contract=off kapcsolót használ. Az explicit fma(a,b,c) továbbra is egyszeri kerekítésű műveletet kér.

A GCC/Clang -ffast-math kapcsolója nem pusztán gyorsabb végrehajtást választ. Verzió- és célfüggő engedményeket kapcsolhat be: átrendezheti a műveleteket, feltételezheti a NaN és infinity hiányát, figyelmen kívül hagyhatja az előjeles nullát vagy a dinamikus kerekítést. Egyes profilok a subnormál kezelést is módosítják. Emiatt akár egy isfinite-alapú adatellenőrzés szerződése is sérülhet.

A reprodukálhatósági jegyzőkönyv rögzítse a fordító és libc verzióját, cél-ABI-t, kapcsolókat, kerekítési módot, contraction-szabályt, subnormál kezelést és releváns float.h értékeket. Ez szükséges, de a matematikai könyvtárak eltérő közelítései miatt nem mindig elegendő a platformok közötti bitazonossághoz.

4.3.5 Hexadecimális literálok és kiírás

A hexadecimális lebegőpontos literálban a p utáni kitevő kettes alapú:

$$\texttt{0x1.8p+1} =\left(1+\frac8{16}\right)2^1=3.$$

A 0x1.999999999999ap-4 a binary64-ben szokásos, $0{,}1$-hez legközelebbi értéket írja le pontosan. Ez alkalmas regressziós tesztek bemeneteinek rögzítésére.

A %a hexadecimális kiírás kettes alapú profilban alapértelmezett pontossággal pontosan megkülönböztethető alakot ad. Decimális visszaolvasáshoz használjuk például:

/* Részlet; C17, <float.h> és <stdio.h> szükséges. */
printf("%.*g\n", DBL_DECIMAL_DIG, value);

A value itt double. long double esetén LDBL_DECIMAL_DIG és %.*Lg, illetve %La szükséges. A printf variadikus hívásában a float automatikusan double-re promótálódik.

A reprodukálható szöveghez a locale-t, a visszaolvasási eljárást és a kerekítési környezetet is rögzíteni kell. A fenti teljes program nem módosít locale-t, így a kezdeti C locale-t használja.

4.3.6 Gyakori hibák

• „A double pontosan tárolja a decimális adatot.” Csak bizonyos értékeket; pontos elszámoláshoz gyakran skálázott egész reprezentáció kell.
• „Minden == hibás.” Pontosan reprezentálható vagy diszkrét állapotértékeknél helyes lehet. A toleranciás közelség viszont nem tranzitív, ezért nem általános egyenlőségi reláció.
• „Az epsilon a mérési tolerancia.” Az epsilon reprezentációs tulajdonság; a tolerancia alkalmazási szerződés.
• „A több tizedesjegyes kiírás pontosabbá teszi az adatot.” Csak a már tárolt értéket mutatja részletesebben.
• „A fast-math eredmény ugyanaz, csak gyorsabb.” A megengedett numerikus viselkedés is változik.

4.3.7 Feladatok

1. Magyarázó: mekkora a binary64 értékek távolsága közvetlenül $2^{32}$ másodperc fölött? Megőrizhető minden 100 ns-os eltérés?

2. Kódolási: egészítsd ki a példát tesztekkel. Ellenőrizd a kétféle nullát, egy negatív késleltetés elutasítását, az abszolút tolerancia szerepét és egy túl nagy relatív tolerancia elutasítását.

3. Diagnosztikai: legyen a = 1e16, b = -1e16, c = 1.0. Binary64, legközelebbi kerekítés és külön műveletek mellett miért különbözhet (a + b) + c és a + (b + c)? Mit következtethetsz egy fast-math fordításból?

4.3.8 Megoldások

4.3.8.1 Időbélyeg felbontása

Itt $e=32$, ezért:

$$\Delta=2^{32-52}=2^{-20}\text{ s} \approx953{,}674\text{ ns}.$$

Nem őrizhető meg minden 100 ns-os eltérés. Ellenőrzésként a nextafter(t, DBL_MAX) - t különbség vizsgálható t = 0x1p32 mellett, a rögzített profilban.

4.3.8.2 A tolerancia tesztelése

Az eredeti C17-programhoz adható részlet:

#include <assert.h>

static void test_latency_close(void)
{
    assert(latency_close(0.0, -0.0, 0.0, 0.0));
    assert(!latency_close(-0.1, 0.1, 1.0, 0.0));
    assert(latency_close(0.0, 5e-10, 1e-9, 0.0));
    assert(!latency_close(0.0, 5e-10, 0.0, 1e-9));
    assert(!latency_close(1.0, 1.0, 0.0, 2.0));
}

Hívd meg a függvényt a main elején, és ne definiáld az NDEBUG makrót. Az ellenőrzési feltétel: minden állítás teljesül. A harmadik és negyedik teszt különíti el az abszolút és relatív küszöb szerepét.

4.3.8.3 Műveleti sorrend

Az a + b pontosan nulla, így az első kifejezés eredménye $1$. A $10^{16}$ körüli binary64 lépésköz $2$; a b + c félúton lévő pontos eredménye pároshoz kerekítéskor visszakerekedik b-re. A második kifejezés így $0$.

Fast-math mellett az átrendezés megengedett lehet, ezért a zárójelezésből önmagában nem következik ez az eredménypár. Azonos kimenet sem bizonyítja a kapcsoló ártalmatlanságát: csak az adott fordítás megfigyelése.

4.3.9 Összefoglaló

A lebegőpontos pontosság a nagyságrendtől, reprezentációtól és műveleti sorrendtől függ. Megbízható mérési kódhoz tartomány- és toleranciaszerződés, különlegesérték-kezelés, valamint rögzített fordítási profil szükséges. A pontos kiírás láthatóvá teszi a tárolt értéket; az elveszett információt azonban sem több számjegy, sem későbbi típusbővítés nem állítja vissza.