Întreprindem aici o sinteză ideatică, lejeră și incompletă, vizând rosturile (oricând, se poate pleca de la [2] pentru referințe istorice precise și pentru diverse aspecte legate de aritmetica în „virgulă flotantă”). Precizăm totuși că ne-am apucat de aceasta nu de dragul artei, ci în contextul lucrului din [1], întâmpinând anumite dileme asupra corectitudinii rezultatelor obținute în diverse limbaje, operând cu numere foarte mari.

Geneza scurtată a lucrurilor

Cam știm de numerele reale: manualele școlare expiră pe rând mulțimile de numere $\mathbb{N}$, $\mathbb{Z}$, $\mathbb{Q}$ și $\mathbb{R}$ și $\mathbb{C}$ (împreună cu structurile matematice subiacente); de fapt, termenul a fost introdus pe când calculele concrete începuseră să evidențieze și valori „nereale” (tipic — valoarea imposibilă "$\sqrt{-1}$" care apăruse pe neașteptate, în calculul rădăcinilor ecuației de gradul trei).
În practică (unde măsurăm, estimăm și înregistrăm distanțe, durate etc.), exprimăm și lucrăm cu aproximări, convenind tacit asupra unei baze comune de numerație — uzual 10, dar pentru sistemele electronice de calcul 2 (sau poate, o anumită putere de 2).

Orice număr real $x$ are o parte întreagă $E_x=\lfloor{x}\rfloor$, care în baza fixată $\beta$ se reprezintă cu un număr finit de cifre, $E_x=\pm\,\overline{a_1a_2\ldots a_k}=\pm\sum_{i=1}^{k}a_i\beta^i\in\mathbb{Z}$ și o parte fracționară $M_x=\{x\}$ care se reprezintă cu un număr eventual nesfârșit de cifre, $M_x=\overline{a_{-1}a_{-2}a_{-3}\ldots}$ cu valoarea $M_x=\sum_{i=1}^\infty a_{-i}\beta^{-i}\in\boldsymbol{[}0,1\boldsymbol{)}$ — unde $a_j$ sunt cifre din baza $\beta$ (întregii din intervalul $[0, \beta-1]$).
$E_x$ caracterizează „mărimea”, adică intervalul de întregi consecutivi în care se află $x$, iar $M_x$ precizează apoi, în acel interval, valoarea respectivă (care de fapt… este greu de „precizat”, mai ales dacă $x\in\mathbb{R}-\mathbb{Q}$).

Avem $x=\lfloor x\rfloor + \{x\}$, iar scrierea obișnuită este $x=E_x\boldsymbol{.}M_x$ — separând convențional prin "." (sau "virgulă") cifrele celor două părți (și subânțelegând baza $\beta$).
Putem „simplifica” această notație, ca în acest exemplu ad hoc: $579\boldsymbol{.}13=5\times \boldsymbol{10^2}+7\times 10^1+$ $9\times 10^0 + 1\times 10^{-1}+3\times 10^{-2}=\boldsymbol{10^2}\left(5 + 7\times 10^{-1}+9\times 10^{-2} + 1\times 10^{-3}+3\times 10^{-4}\right)=5\boldsymbol{.}7913\times \boldsymbol{10^{2}}$; factorizând prin cea mai mare putere de $\beta$ existentă în scrierea obișnuită "$E_x\boldsymbol{.}M_x$" — altfel spus, mutând (sau "flotând") virgula spre stânga (sau spre dreapta, după caz), până ce în fața ei rămâne o singură cifră, dar nenulă — vedem că orice număr real $x\boldsymbol{\ne 0}$ se reprezintă în mod unic (până la o anumită excepție, binecunoscută) sub forma $x=\pm \boldsymbol{M\times\beta^E}$, unde $M$ (numit semnificant, sau chiar „mantisă”) este astfel încât $\boldsymbol{1\le M < \beta}$ iar $E$ (numit expozant, sau exponentul lui $x$) este un număr întreg. (un exemplu mai potrivit era $0.005793\approx 5.79\times 10^{-3}$).
Dacă începând de la un rang încolo, cifrele din $M_x$ sunt toate egale cu $\beta-1$, atunci avem totuși două reprezentări pentru $x$; de exemplu trebuie să acceptăm că $x=8\boldsymbol{.}62999\ldots$ $\boldsymbol{=}8\boldsymbol{.}63$ (într-adevăr: $8\boldsymbol{.}62+\dfrac{9}{10^3}\sum_{j=0}^{\infty}10^{-j}=8\boldsymbol{.}62+\dfrac{9}{10^3}\dfrac{1}{1-10^{-1}}=8\boldsymbol{.}62+10^{-2}=8\boldsymbol{.}63$).

Având în vedere limitările inerente calculatorului (uman, sau… mai ales artificial), parametrii $M$ și $E$ trebuie reduși sau adaptați la câte un anumit număr de cifre — ajungând astfel la numerele „flotante” cu precizia dată $\boldsymbol{p}$ (sau „numerele-mașină”): cele pentru care semnificantul $M$ are (maximum) $p$ cifre (în baza $\beta$), iar expozantul $E$ este cuprins între două limite convenite (sau impuse) în prealabil.
Într-un astfel de format $\boldsymbol{(\beta, p)}$, un număr nenul oarecare (dacă nu depășește cadrul fixat pentru expozanți) va putea fi reprezentat (aproximativ) printr-unul dintre cele două numere flotante care încadrează valoarea reală respectivă; desigur… pot apărea fel de fel de îndoieli (și riscuri, pentru practică: "It makes me nervous to fly on airplanes since I know they are designed using floating-point arithmetic.") și devin necesare anumite convenții și reguli de rotunjire.

Mantisele posibile, cu $p$ cifre

Este ușor de generat mantisele posibile în precizia $p$, "$a_0\boldsymbol{.}a_1a_2\ldots a_{p-1}$", cu $a_0\ne 0$:

gen_fractions <- function(B, P) {  # baza și precizia totală
    dig <- 0:(B-1)  # vectorul cifrelor în baza B (aici, B <= 10)
    M <- expand.grid(rep(list(dig), times = P-1)) |> 
         as.matrix() |> 
         apply(1, paste, collapse="") |> 
         sort()  
    sapply(1:(B-1), function(d) paste0(d, ".", M)) |> 
    as.vector()  # vectorul mantiselor (de tip "char"), ordonat crescător
}

Precizăm că expand.grid() furnizează un „tabel” în care fiecare linie conține câte o combinație de valori ale vectorilor indicați — aici, P-1 vectori, egali toți cu dig; rezultă toate șirurile de câte P-1 cifre și prin paste(), le-am adăugat cifra nenulă dinaintea punctului.

Pentru exemplificare, ne amintim regula de calcul „cu două zecimale exacte”, folosită pentru stabilirea mediilor școlare anuale — și considerăm sistemul $\boldsymbol{(\beta=10,\, p=3)}$ care ar corespunde acestei reglementări. Avem 900 de mantise, de la $1.00$ la $9.99$:

> M <- gen_fractions(10,3) |> print(quote=FALSE)
  [1] 1.00 1.01 1.02 1.03 1.04 1.05 1.06 1.07 1.08 1.09 1.10 1.11 1.12 1.13 1.14
 [16] 1.15 1.16 1.17 1.18 ...
 ... ...
 [91] 1.90 1.91 1.92 1.93 1.94 1.95 1.96 1.97 1.98 1.99 2.00 2.01 2.02 2.03 2.04
... ...
[886] 9.85 9.86 9.87 9.88 9.89 9.90 9.91 9.92 9.93 9.94 9.95 9.96 9.97 9.98 9.99

Calculată exact, media școlară poate fi $8.457$ și avem de ales între cele două valori din formatul $(10, p=3)$ care o încadrează: $8.45$ sau $8.46$. Dacă media exactă ar fi $9.666...$, am avea de ales între $9.66$ și $9.67$; iar dacă media exactă ar fi $7.4999...$, avem de ales între $7.49$ și $7.50$. Alegerea aproximării dintre cele două posibile pentru precizia $p$ respectivă, depinde de reglementările existente („cu două zecimale exacte” sugerează $8.45$, $9.66$, respectiv $7.49$; dar dacă ne luăm după un calculator obișnuit, valorile alese vor fi $8.46$, $9.67$, respectiv $7.50$).
Putem considera și alți exponenți decât zero (cel specific pentru mediile școlare, aflate în intervalul $[1,\,10)$); de exemplu, $15.349$ și $0.07894$ s-ar reprezenta (cu $p=3$ cifre semnificative) prin $1.53\times 10^1$ (sau prin $1.54\times 10^1$), respectiv $7.89\times 10^{-2}$ (sau $7.90\times 10^{-2}$).

Să observăm că în cazul B=2, condiția $a_0\ne 0$ obligă $a_0\boldsymbol{=1}$ — astfel că în acest caz este suficient să vizăm numai „fracția” de după virgulă, pe care în treacăt o putem interpreta și ca valoare de număr întreg cu $p-1$ cifre.

Formate binare standard

Bineînțeles că pentru calculul mediilor școlare chiar nu era necesară, inventarea numerelor „flotante”… Cazul cel mai important este $\beta=\boldsymbol{2}$, asumat de mașinile electronice de calcul; acestea operează pe câmpuri binare de o anumită lungime — cel mai obișnuit 32, 64 sau chiar 128 de biți — structurate (ca „numere flotante”, sau „numere-mașină”) pe trei subcâmpuri: primul bit $s$ înregistrează semnul numărului, pe următorii câțiva biți se reprezintă expozantul $E$, iar pe biții rămași se reprezintă semnificantul $f$:

Bitul de semn disociază între numerele pozitive (cu $s=0$) și cele negative ($s=1$); compararea a două numere astfel reprezentate revine la a compara lexicografic cele două câmpuri de câte $32\,/\,64\,/\,128 \text{ biți}$ (ținănd seama totuși, de semn: două numere negative sunt în ordine inversă aceleia dintre valorile absolute ale lor).

Exponentul maxim ar fi $2^n$, unde $n=\boldsymbol{8}\,/\,\boldsymbol{11}\,/\,15$; dar avem și exponenți negativi, deci primul bit din câmpul $E$ ar trebui rezervat pentru a înregistra semnul… Pentru a evita separarea acestui bit de semn (și a-l economisi), s-a convenit ca în câmpul $E$ să se înregistreze deplasamentul exponentului față de $2^{n-1}$ (mijlocul intervalului $[0,2^n$]). Mai precis, în loc de un expozant (admisibil) $h$, în câmpul $E$ se memorează valoarea $h\boldsymbol{+}(\boldsymbol{2^{n-1}-1})\in\boldsymbol{[1, 2^n-2]}$; valorile devenite „extreme”, $E=0$ și $E=2^n-1$ vor putea indica anumite situații speciale (de exemplu, depășirea domeniului specific formatului respectiv).

Numărul zero este reprezentat prin $E=0$ și $f=0$, deci prin două valori („zero cu semn”): fie $\boldsymbol{+}0$, fie $\boldsymbol{-}0$ (după cum bitul de semn $s$ este $0$ sau $1$). Ca semnificație, $+0$ reprezintă un număr real pozitiv care este mai mic decât cel mai mic număr pozitiv reprezentabil în formatul respectiv (analog pentru $-0$).
Dacă $f=0$ și toți biții din câmpul $E$ sunt egali cu $1$, atunci avem valorile $+\infty$ (semnificând un număr real mai mare decât valoarea maximă reprezentabilă în formatul respectiv) și $-\infty$; pentru împărțire este asigurată proprietatea $1/\pm 0=\pm\infty$.
Să observăm că ne-a rămas cam un singur caz: toți biții din $E$ sunt egali cu $1$, iar în câmpul $f$ există măcar un bit $1$; se reprezintă astfel valorile NaN ("Not a Number"), semnificând nedeterminări precum $+\infty-\infty$, sau valori nedefinite (precum $\sqrt{-1}$).

Precizia asigurată de aceste formate este $p=\boldsymbol{24}\,/\,\boldsymbol{53}\,/\,113$ biți — fiindcă în câmpul $f$ nu se înregistrează și primul bit ($a_0$, dinaintea virgulei), dat fiind că acesta (trebuind să fie o cifră nenulă, în baza respectivă) este obligatoriu $\boldsymbol{1}$ (și fiind cunoscut, consemnarea lui devine inutilă).
Dar pot fi considerate și numere subnormale, în care $f$ conține măcar un bit $1$ dar este prefixat (în fața virgulei) nu de un bit $1$ cum este cazul normal (când $E\ne 0$), ci de un bit $\boldsymbol{0}$; mai precis, „numerele subnormale” au $E=0$ și $f\ne 0$ — doar că în acest caz, acel bit $1$ care în mod tacit prefixează partea $f$ (în cazul numerelor normale) este transferat tacit părții $E$, care ajunge astfel să fie interpretat ca fiind $E=1$ (valoarea minimă de exponent admisibil), încât exponentul real $h$ rezultat astfel (din ecuația $E=h+(2^{n-1}-1)$) este $2-2^{n-1}$.

Deci numerele normale au valorile $(-1)^s\cdot 2^\boldsymbol{h}\cdot \boldsymbol{1.}f$ unde $h = E-(2^{n-1}-1)$ cu $E$ întreg și $1\le E \le 2^n -1$, iar $f$ conține $(p-1)$ biți. Cu $s=0$, $E=1$ și $f=0$ obținem cel mai mic număr normal pozitiv: $\boldsymbol{2^{-126}}\approx 1.175\times 10^{-38}$ pentru formatul pe 32 de biți și $\boldsymbol{2^{-1022}}\approx$ $2.225\times 10^{-308}$ pentru formatul "binary64".
Pentru cel mai mare număr normal pozitiv $h=(2^n-2)-(2^{n-1}-1)=2^{n-1}-1$, iar $f$ are toți biții egali cu $1$ (adică $f=1 + 2^{-1} + 2^{-2} + \cdots + 2^{-(p-1)}=2(1-2^{-p})$) — rezultă valorile $\boldsymbol{2^{128}\cdot (1-2^{-24})}\approx 3.402823\times 10^{38}$, respectiv $\boldsymbol{2^{1024}\cdot (1-2^{-53})}\approx 1.797693\times 10^{308}$.
Valorile respective pot fi determinate folosind de exemplu consola din $\mathbf{R}$:

> n <- 11  # pentru formatul "double precision" (pe 64 biți)
> h <- 2^(n-1)
> (2^h)*(1-2^(-53))  # corect era: (2.0^h)*(1-2^(-53))
[1] Inf  # factorul întreg 2^1024 este prea mare!  
> (2^(h-1))*(2-2^(-52))  # decrementăm exponentul și dublăm mantisa  
[1] 1.797693e+308  # cel mai mare număr normal în formatul "binary64"  

Am folosit un artificiu de rescriere tipic, pentru a evita calculul unei puteri prea mari de număr întreg — ceea ce nu era necesar, dacă în loc de "2^h" foloseam "2.0^h", forțând double, din integer… însă am vrut să ocolim aritmetica specifică numerelor flotante (angajând numai aritmetica de bază, cu numere întregi: $2^h$ decurge ca înmulțire cu baza, adică deplasare spre stânga).

Numerele subnormale au valorile $(-1)^s\cdot 2^\boldsymbol{h}\cdot \boldsymbol{0.}f$, unde $h=2-2^{n-1}$, iar $f$ conține măcar un bit $1$. Setând în $f$ numai ultimul bit (cel de rang $p-1$), avem cel mai mic număr pozitiv subnormal: $2^{-126}\cdot 2^{-23}=2^{-149}\approx 1.401\times 10^{-45}$ respectiv, $2^{-1022}\cdot 2^{-52}=2^{-1074}\approx 5\times 10^{-324}$.
Iar cel mai mare număr subnormal se obține setând toți cei $(p-1)$ biți din $f$, deci are valoarea $\boldsymbol{2^h\cdot (1-2^{p-1})}$, adică $\approx 1.175494\times 10^{-38}$, respectiv $\approx 2.225074\times 10^{-308}$.

Asupra distribuției numerelor-mașină

Pentru simplificarea exprimării, considerăm numai numerele pozitive; fie $\mathbb{F}^+$ mulțimea numerelor pozitive reprezentabile normal (dar în unele contexte, subânțelegem în $\mathbb{F}^+$ și pe cele subnormale) în formatul respectiv (binary32, sau binary64).

Mantisa are $p$ biți (24, respectiv 53), dar primul bit este fixat (pe 1 pentru numerele normale, sau pe 0 pentru cele subnormale); mulțimea mantiselor (văzute ca șabloane binare de lungime dată) este în mod implicit, ordonată lexicografic (și dacă le interpretăm ca reprezentări binare de numere întregi, atunci fiecare se obține adunând 1 celei precedente).
Deci pentru fiecare exponent $h$ fixat, în $\mathbb{F}^+$ avem câte exact $2^{p-1}$ numere cuprinse în intervalul $[2^h, 2^{h+1})$, egal distanțate între ele; raportând lungimea intervalului la numărul de valori, rezultă că distanța este egală cu $2^h/2^{p-1}=2^h\times \boldsymbol{\varepsilon}$ (și crește exponențial odată cu $h$), unde factorul constant $\boldsymbol{\varepsilon=2^{-(p-1)}}$ (ponderea ultimului bit al mantisei) este numit de obicei "machine epsilon".

Cazul tipic ar fi $h=0$: numerele din $\mathbb{F^+}$ aflate în intervalul $[1,2)$ — fiind distanțate între ele cu $\varepsilon$ — sunt termenii progresiei aritmetice $\varphi_k=\boldsymbol{1+k\varepsilon},\,0\le k < 2^{p-1}-1$; dacă $x$ este un număr real oarecare din $[1,2)$, atunci $x$ se află între doi termeni consecutivi din șirul $(\varphi_k)$ și distanța lui $x$ la cel mai apropiat dintre aceștia este de cel mult $\varepsilon/2$ (altfel spus, eroarea aproximării lui $x$ printr-unul dintre cei doi termeni din $\mathbb{F}^+$ este de cel mult $\varepsilon/2$).
Apoi, crescând exponentul $h$ — respectiv, coborându-l — numerele din $\mathbb{F}^+$ aflate în intervalul $[2,2^2)$ — respectiv, în $[2^{-1}, 1)$ — sunt exact dublele $2(\varphi_k)$ — respectiv, jumătățile $2^{-1}(\varphi_k)$ — celor precedente; dublarea, respectiv înjumătățirea, continuă apoi din aproape în aproape, pentru ceilalți exponenți $h$.

Fiindcă incrementarea lui $h$ dublează distanța precedentă, urmează că în vecinătatea lui $2^h$ apare o anumită discontinuitate; de exemplu, numerele reale din intervalul $(2-\varepsilon, 2)$ s-ar reprezenta prin 2 cu o eroare mai mică decât $\varepsilon/2$, în timp ce pentru cele de la dreapta lui 2 (din intervalul $(2, 2+\varepsilon)$) eroarea ar fi majorată de $(\boldsymbol{2}\varepsilon)/2=\varepsilon$.

Dacă luăm în considerație și numerele subnormale, adică pe cele cu valoarea absolută mai mică decât cel mai mic număr din $\mathbb{F}^+$ — care este $2^{E_{min}-1}$ cu $E_{min}=2-2^{n-1}$ — atunci să observăm că (dat fiind că pentru acestea avem un același expozant, anume chiar $E_{min}$) distanța dintre ele este aceeași cu cea dintre numerele normale aflate în intervalul $[2^{E_{min}-1}, 2^{E_{min}})$, deci este egală cu $2^{E_{min}-p}$.

Dacă $x$ este un număr real pozitiv oarecare, atunci sau $x\in\mathbb{F}^+$ (adică este reprezentabil exact în formatul respectiv), sau $x$ depășește pe cel mai mare număr din $\mathbb{F}^+$ (caz în care este reprezentat prin $INF=+\infty$), sau este situat între două numere consecutive din $\mathbb{F}^+$ — caz în care este necesară o regulă de rotunjire, pentru a decide care dintre acestea să fie considerat ca reprezentant al lui $x$ în formatul respectiv.
Regula de rotunjire cea mai obișnuită este "round-to-nearest-even" (desemnată prin RN(x)), prin care $x$ este reprezentat prin acel număr din $\mathbb{F}^+$ care îi este cel mai apropiat; ar apărea o dilemă dacă $x$ este la mijloc — iar în acest caz se ține seama de faptul că două numere consecutive din $\mathbb{F}^+$ diferă neapărat la ultimul bit al mantisei și de obicei, se alege ca reprezentat al lui $x$ pe cel care are ultimul bit egal cu 0.

Am văzut mai sus că eroarea absolută depinde în vecinătatea unei puteri de 2, de faptul că $x$ este la stânga sau la dreapta acelei puteri; în schimb, eroarea relativă este „constantă”: $\dfrac{|x-\mathrm{RN}(x)|}{x}\le \varepsilon /2$. Într-adevăr, fie $\varphi_j$ și $\varphi_{j+1}$ cei doi termeni consecutivi din progresia aritmetică $\left(2^h(1 + k\varepsilon)\right)_{0\le k < 2^{p-1}-1}$ formată de numerele din $\mathbb{F}^+$, pentru care $\varphi_j \le x \le \varphi_{j+1}$; distanța la cel mai apropiat capăt este $|x-\mathrm{RN}(x)| \boldsymbol{\le} (\varphi_{j+1}-\varphi_j)/2=2^h\varepsilon /2$ și fiindcă, evident, $2^h \le \varphi_j = 2^h(1+j\varepsilon) \le x$, rezultă $|x-\mathrm{RN}(x)|\le x\,\varepsilon/2$.
Echivalent, $x(1-\varepsilon/2)\le \mathrm{RN}(x)\le x(1+\varepsilon/2)$ de unde avem și exprimarea $\mathrm{RN}(x)=x(1+\delta)$ pentru $|\delta| \le \varepsilon/2$. Iar implementările existente pentru operații, asigură că pentru orice numere-mașină $x$ și $y$ și oricare operație $\circ$ dintre cele patru operații aritmetice de bază (plus extragerea rădăcinii pătrate), avem: $\mathrm{RN}(x\circ y)=(x\circ y)(1+\delta)$, pentru un $\delta\le\varepsilon/2$ (adică, rezultatul rotunjit diferă de cel exact cel mult prin factorul $1+\varepsilon/2$).

Experimente de reprezentare

Se pot găsi în multe locuri, diverse exemple de conversie în „virgulă mobilă” a unor numere reale; aici exploatăm consecvent faptul evidențiat mai sus: numerele din $\mathbb{F}^+$ care au un același exponent formează o progresie aritmetică.

Un prim episod

Generăm aleatoriu un exponent (pe 8 biți) și o mantisă (pe 23 de biți) și întrebăm invers: care număr este reprezentat prin numărul-mașină (în formatul binary32) obținut astfel?

gen_binary <- function(n)
    sample(x=c("0","1"), size=n, replace=TRUE) |> paste0(collapse="")

xb <- paste(gen_binary(1), gen_binary(8), gen_binary(23)) |> print()
[1] "0 11100101 01100011010011001110111"  # spațiat numai pentru lizibilitate

Șablonul binar $xb$ reprezintă un număr $x\in\mathbb{F}^+_{32}$; expozantul "11100101" are valoarea $2^7+2^6+2^5+5=229$, deci exponentul real este $229-(2^7-1)=\boldsymbol{102}$, astfel că $x$ este unul dintre termenii progresiei aritmetice $\left(\boldsymbol{2^{102}(1 + k\varepsilon)}\right)_{0\le k < 2^{23}-1}$, unde $\varepsilon=2^{-23}$; valoarea parametrului $k$ este dată de numărul întreg reprezentat pe ultimii 23 de biți din $xb$, anume (folosind eventual o aplicație "Calculator științific") $k=$"01100011010011001110111"$=3253879$.
Rezultă că $x=2^{102}(1+3253879\times 2^{-23})$, care este număr întreg, $x=\boldsymbol{2^{79}\times 11642487}$ (toate numerele din $\mathbb{{F}^+_{32}}$ care sunt mai mari ca $2^{23}$ sunt întregi, din ce în ce mai dispersați pe măsură ce creștem exponentul).
Putem afișa valoarea exactă a lui $x$ folosind Python (unde putem lucra cu întregi lungi):

>>> 2**79 * 11642487
7037451569413832588168691449856  # aproximativ 7.03745157×10^30

Dar afișarea în formatul "floating-point" obișnuit a lui $x$ (implicând și conversie la baza 10) nu poate fi una exactă (în fond, precizia de 24 biți asigură numai 6-7 cifre zecimale exacte):

> x <- 2^102 * (1 + 3253879/2^23)  # în consola R
> x
[1] 7.037452e+30  # primele 6 cifre zecimale exacte, dintre cele 31 (rotunjind la a 7-a)
> options(digits=16)  # pentru afișare cu 15 cifre exacte (rotunjind la a 16-a)
> x
[1] 7.037451569413833e+30  # Obs.: intern, R folosește modelul "binary64"

Putem vedea folosind $\mathbf{C}$ — unde avem și tipul float, corespunzător formatului "binary32" — că totuși $x$ este păstrat exact (cu toate cifrele), confirmând că $x\in\mathbb{F}^+_{32}$:

#include <stdio.h>
int main() {
    union {
        int as_int;
        float as_float;
    } num;
    num.as_int = 0b01110010101100011010011001110111;  /* 0x72B1A677 */
    printf("%d\t%f\n", num.as_int, num.as_float);
}
/* se afișează:
1924245111	7037451569413832588168691449856.000000 */

Șablonul binar $xb$ a fost interpretat o dată ca fiind de tip int (rezultând valoarea $1+2+2^2+2^4+2^5+\cdots+2^{29}+2^{30}=1924245111$) și apoi, ca float (rezultând exact $x$).

Flotantul care urmează după $x$ în $\mathbb{F}^+_{32}$ este $y=2^{102}\left(1 + (\boldsymbol{k+1})\varepsilon\right)=2^{79}\times 11642488=$ 7037452173876742395483278802944 $\approx$ 7.037452173876742e+30; toate numerele reale $z$ cuprinse între $x$ și $y$ — aflate ambele în $\mathbb{F}^+_{32}$, la distanță mare $\approx 6\times 10^{23}$ unul de altul — vor fi aproximate prin $\mathrm{RN}(z)$ fie cu $x$, fie cu $y$ (și distingerea între ele a acestor valori $z$ impune mărirea temporară a preciziei inițiale $p$ — adică, lucrul „în precizie multiplă”).

Un al doilea episod

Putem folosi și „invers”, programul $\mathbf{C}$ redat mai sus — pornim de la partea num.as_float:

    num.as_float = 1.0 / 10;  /* un număr "real" oarecare... 
    printf("%d\t%f\n", num.as_int, num.as_float);
/* se afișează:
1036831949	0.100000 */

Drept rezultat al împărțirii la 10 a numărului cu tipul implicit float 1.0, s-a afișat (cam cum ne așteptăm) 0.10000… Afișarea este orientată spre cel „cuminte” — însă valoarea afișată nu este totuna cu aceea memorată.

Este ușor de văzut că $0.1_{(10)}=\sum_{j=1}^\infty(2^{-4j}+2^{-4j-1})=0.000\overline{1100}_{(2)}$ unde secvența de 4 biți supraliniată se repetă în continuare la nesfârșit; prin urmare, inversul lui $10$ nu face parte din $\mathbb{F}^+_{32}$ (și în memorie, el este aproximat printr-unul dintre cele două numere vecine din $\mathbb{F}^+_{32}$ între care este încadrat).

Dacă $h$ este exponentul celor două numere din $\mathbb{F}^+_{32}$ între care se află $10^{-1}$, atunci pentru un anumit $k$, $0\le k < 2^{p-1}-1$ (unde $p$ este precizia corespunzătoare formatului respectiv, deci aici $p=24$) trebuie să avem:

de unde găsim $k=\boldsymbol{\left\lfloor\varepsilon^{-1}(2^{-h}\cdot 10^{-1}-1)\right\rfloor}$.

Putem găsi exponentul $h$ plecând de la valoarea afișată mai sus pentru num.as_int: 1036831949=0x3DCCCCCD=0b00111101110011001100110011001101; partea de exponent este 0b01111011=123, deci exponentul adevărat este $h=123-127=-4$ (bineînțeles… puteam afla exponentul și direct: avem $2^3 < 10 < 2^4$ și inversând, regăsim $h=-4$).

Rezultă $k=\left\lfloor 2^{23}\dfrac{3}{5}\right\rfloor=5033164$; prin urmare, încadrarea lui $10^{-1}$ între numere consecutive din $\mathbb{F}^+_{32}$ este: $$2^{-4}(1 + 5033164\cdot\varepsilon) < \dfrac{1}{10} < 2^{-4}(1 + 5033165\cdot\varepsilon)$$

ceea ce (încercând să folosim numai aritmetica întregilor) revine la:

Rămâne să comparăm numerele întregi rămase la capete, cu valoarea $2^{27}/10$ (care se obține banal, punând virgulă în fața ultimei cifre); desigur, apelăm la interpretorul Python:

>>> 2**27 / 10
13421772.8
>>> 5033164 + 2**23
13421772  # diferența este -0.8
>>> 5033165 + 2**23
13421773  # diferența este +0.2

Rezultă că $10^{-1}$ este mai apropiat de capătul din dreapta, deci

În unele contexte, în loc de $\mathrm{RN}()$ am putea prefera regula "round-down", prin care s-ar produce capătul din stânga:

Repetând calculele cu $p=53$ în loc de $p=24$ (caz în care va rezulta o fracție de întregi diferită de cea de mai sus) putem găsi și reprezentarea în $\mathbb{F}^+_{64}$ a lui $10^{-1}$, confruntând apoi cu rezultatul de mai sus (constatând că mărind precizia, crește și apropierea de valoarea reală).