Voici du code qui gère tous les nombres finis. Il attend l'arithmétique IEEE 754. J'ai remplacé ma version précédente par un code plus simple et plus propre. Au lieu de deux implémentations du calcul de distance, celui-ci a deux implémentations de conversion d'un nombre en virgule flottante en son encodage (une en copiant les bits, une en le manipulant mathématiquement). Ensuite, le calcul de la distance est assez simple (les valeurs négatives doivent être ajustées, puis la distance est simplement une soustraction).
#include
#include
#include
#include
#include
typedef double Float; // Le type en virgule flottante à utiliser.
typedef std::uint64_t UInt; // Entier non signé de même taille que Float.
/* Définir une valeur avec seulement le bit de poids fort d'un UInt défini. C'est aussi l'encodage du -0 en virgule flottante.
*/
static constexpr UInt HighBit
= std::numeric_limits::max() ^ std::numeric_limits::max() >> 1;
// Retourne l'encodage d'un nombre en virgule flottante en copiant ses bits.
static UInt EncodingBits(Float x)
{
UInt result;
std::memcpy(&result, &x, sizeof result);
return result;
}
// Retourne l'encodage d'un nombre en virgule flottante en utilisant des mathématiques.
static UInt EncodingMath(Float x)
{
static constexpr int SignificandBits = std::numeric_limits::digits;
static constexpr int MinimumExponent = std::numeric_limits::min_exponent;
// Encoder le bit de poids fort.
UInt result = std::signbit(x) ? HighBit : 0;
// Si la valeur est zéro, les bits restants sont zéros, nous avons donc fini.
if (x == 0) return result;
/* La bibliothèque C fournit une routine peu connue pour diviser un nombre en virgule flottante en une mantisse et un exposant.
Notez que cela produit une mantisse normalisée, pas l'encodage réel de la mantisse. Notamment,
cela amène les mantisses des sous-normaux à au moins 1/2. Nous ajusterons cela ci-dessous. De plus, cette routine normalise à [1/2, 1),
alors que l'IEEE 754 est généralement exprimée avec [1, 2), mais cela ne nous dérange pas.
*/
int xe;
Float xf = std::frexp(std::fabs(x), &xe);
// Tester si le nombre est subnormal.
if (xe < MinimumExponent)
{
/* Pour une valeur subnormale, l'encodage de l'exposant est zéro, donc nous devons
uniquement insérer les bits de la mantisse. Cela met à l'échelle la mantisse
pour que son bit de poids faible soit mis à l'endroit de la position 1 puis il l'insère
dans l'encodage.
*/
result |= (UInt) std::ldexp(xf, xe - MinimumExponent + SignificandBits);
}
else
{
/* Pour une valeur normale, la mantisse est encodée sans son bit de poids fort.
Donc, nous soustrayons .5 pour supprimer ce bit puis mettons à l'échelle la
mantisse pour que son bit de poids faible soit mis à l'endroit de la position 1.
*/
result |= (UInt) std::ldexp(xf - .5, SignificandBits);
/* L'exposant est encodé avec un biais de (dans la terminologie C++)
MinimumExponent - 1. Nous le soustrayons donc pour obtenir l'encodage de l'exposant
puis le décalons à la position du champ de l'exposant. Ensuite, nous l'insérons dans l'encodage.
*/
result |= ((UInt) xe - MinimumExponent + 1) << (SignificandBits-1);
}
return result;
}
/* Retourne l'encodage d'un nombre en virgule flottante. Pour l'illustration,
on obtient l'encodage avec deux méthodes différentes et on compare les résultats.
*/
static UInt Encoding(Float x)
{
UInt xb = EncodingBits(x);
UInt xm = EncodingMath(x);
if (xb != xm)
{
std::cerr << "Erreur interne lors de l'encodage de " << x << ".\n";
std::cerr << "\tEncodingBits dit " << xb << ".\n";
std::cerr << "\tEncodingMath dit " << xm << ".\n";
std::exit(EXIT_FAILURE);
}
return xb;
}
/* Retourne la distance de a à b sous forme du nombre de valeurs représentables en
Float d'une à l'autre. b doit être supérieur ou égal à a. 0 est
compté une seule fois.
*/
static UInt Distance(Float a, Float b)
{
UInt ae = Encoding(a);
UInt be = Encoding(b);
/* Pour les valeurs représentées de +0 à l'infini, les points flottants binaires IEEE 754
sont en ordre croissant et sont consécutifs. Donc on peut simplement soustraire deux encodages pour obtenir le nombre de valeurs représentables
entre eux (en incluant un point de départ mais pas l'autre).
Malheureusement, les nombres négatifs ne sont pas consécutifs et ont
une direction différente. Pour faire face à cela, si le nombre est négatif, nous transformons
son encodage en soustrayant de l'encodage du -0. Cela nous donne une
séquence consécutive d'encodages allant du plus grand nombre négatif fini de grandeur au plus grand nombre fini, en ordre croissant
sauf pour l'enroulement à la valeur UInt maximale.
Notez que cela mappe également l'encodage du -0 sur 0 (l'encodage du +0),
donc les deux zéros deviennent un seul point, donc ils sont comptés une seule fois.
*/
if (HighBit & ae) ae = HighBit - ae;
if (HighBit & be) be = HighBit - be;
// Retourner la distance entre les deux encodages transformés.
return be - ae;
}
static void Try(Float a, Float b)
{
std::cout << "[" << a << ", " << b << "] contient "
<< Distance(a,b) + 1 << " valeurs représentables.\n";
}
int main(void)
{
if (sizeof(Float) != sizeof(UInt))
{
std::cerr << "Erreur, UInt doit être un entier non signé de la même taille que Float.\n";
std::exit(EXIT_FAILURE);
}
/* Préparer quelques valeurs de test : plus petite valeur positive (subnormale), plus grande valeur
subnormale, plus petite valeur normale.
*/
Float S1 = std::numeric_limits::denorm_min();
Float N1 = std::numeric_limits::min();
Float S2 = N1 - S1;
// Test 0 <= a <= b.
Try( 0, 0);
Try( 0, S1);
Try( 0, S2);
Try( 0, N1);
Try( 0, 1./3);
Try(S1, S1);
Try(S1, S2);
Try(S1, N1);
Try(S1, 1./3);
Try(S2, S2);
Try(S2, N1);
Try(S2, 1./3);
Try(N1, N1);
Try(N1, 1./3);
// Test a <= b <= 0.
Try(-0., -0.);
Try(-S1, -0.);
Try(-S2, -0.);
Try(-N1, -0.);
Try(-1./3, -0.);
Try(-S1, -S1);
Try(-S2, -S1);
Try(-N1, -S1);
Try(-1./3, -S1);
Try(-S2, -S2);
Try(-N1, -S2);
Try(-1./3, -S2);
Try(-N1, -N1);
Try(-1./3, -N1);
// Test a <= 0 <= b.
Try(-0., +0.);
Try(-0., S1);
Try(-0., S2);
Try(-0., N1);
Try(-0., 1./3);
Try(-S1, +0.);
Try(-S1, S1);
Try(-S1, S2);
Try(-S1, N1);
Try(-S1, 1./3);
Try(-S2, +0.);
Try(-S2, S1);
Try(-S2, S2);
Try(-S2, N1);
Try(-S2, 1./3);
Try(-N1, +0.);
Try(-N1, S1);
Try(-N1, S2);
Try(-N1, N1);
Try(-1./3, 1./3);
Try(-1./3, +0.);
Try(-1./3, S1);
Try(-1./3, S2);
Try(-1./3, N1);
Try(-1./3, 1./3);
return 0;
}
