Skip to content

Instantly share code, notes, and snippets.

@christophercrouzet

christophercrouzet/distance.h

Last active April 16, 2024 13:02
Show Gist options
  • Save christophercrouzet/464977b435c958765c10 to your computer and use it in GitHub Desktop.
Save christophercrouzet/464977b435c958765c10 to your computer and use it in GitHub Desktop.
Download ZIP
Compute the distance in ulps between floating-point numbers in C++11.
Raw
distance.h
#ifndef DISTANCE_H
#define DISTANCE_H
#include <algorithm>
#include <cassert>
#include <cmath>
#include <limits>
#include "floatingpointtraits.h"
namespace internal
{
template<typename T>
T
epsilonForExponent(int exponent)
{
using Traits = traits::FloatingPoint<T>;
return exponent < Traits::smallestNormalNumberExponent ?
// The epsilon is constant amongst the subnormal numbers.
Traits::subnormalNumbersEpsilon :
// The amount of numbers per exponent is constant amongst the
// normal numbers and is equal to `2^bitsForSignificand`, with
// `bitsForSignificand` being the number of bits used to define
// the significand of a floating-point number.
//
// Knowing the amount of normal numbers that can be representated for
// a given exponent `x`, the epsilon is equal to the absolute
// difference `2^(x + 1) - 2^x` divided by the amount of numbers.
//
// After simplification, `(2^(x + 1) - 2^x) / 2^bitsForSignificand`
// becomes `2^(x - bitsForSignificand)`.
std::exp2(T(exponent - std::numeric_limits<T>::digits + 1));
}
template<typename T>
T
distanceTowardsExponentUpperBound(T value, int exponent)
{
return (std::exp2(T(exponent + 1)) - value)
/ epsilonForExponent<T>(exponent);
}
template<typename T>
T
distanceTowardsExponentLowerBound(T value, int exponent)
{
return (value - std::exp2(T(exponent)))
/ epsilonForExponent<T>(exponent);
}
template<typename T>
T
distance(T from, T to)
{
using Traits = traits::FloatingPoint<T>;
assert(std::isfinite(from));
assert(std::isfinite(to));
assert(from != to);
assert(from < to);
assert(from > 0);
assert(to > 0);
assert(std::signbit(from) == std::signbit(to));
int exponentFrom = std::ilogb(from);
int exponentTo = std::ilogb(to);
T out = 0;
if (exponentFrom == exponentTo) {
out = (to - from) / epsilonForExponent<T>(exponentFrom);
}
else {
// 2^x 2^(x+1) 2^(x+2) 2^(x+3)
// | | | | | | |
// from to
//
// \___________/
// ^
// distance between
// `2^(exponentFrom + 1)`
// and `2^exponentTo`
// Compute the distance between `2^(exponentFrom + 1)`
// and `2^exponentTo`.
assert(exponentFrom < exponentTo);
if (exponentTo - exponentFrom > 1) {
// Do every in-between exponents represent normal numbers?
if (exponentFrom + 1 >= Traits::smallestNormalNumberExponent) {
out += Traits::amountOfNormalNumbersPerExponent
* (exponentTo - exponentFrom - 1);
}
// Do every in-between exponents represent subnormal numbers?
else if (exponentTo <= Traits::smallestNormalNumberExponent) {
int exponent = exponentFrom + 1;
while (exponent < exponentTo) {
out += std::exp2(
T(exponent - Traits::smallestSubnormalNumberExponent)
);
++exponent;
}
}
else {
int exponent = exponentFrom + 1;
// Distance within subnormal numbers.
while (exponent < Traits::smallestNormalNumberExponent) {
out += std::exp2(
T(exponent - Traits::smallestSubnormalNumberExponent)
);
++exponent;
}
// Distance within normal numbers.
out += Traits::amountOfNormalNumbersPerExponent
* (exponentTo - exponent);
}
}
// 2^x 2^(x+1) 2^(x+2) 2^(x+3)
// | | | | | | |
// from to
//
// \_/ \__________/
// ^ ^
// distance towards distance towards
// exponent exponent
// upper bound lower bound
// Add both values at once to reduce the precision error
// for large numbers.
out += (
distanceTowardsExponentUpperBound(from, exponentFrom)
+ distanceTowardsExponentLowerBound(to, exponentTo)
);
}
assert(out >= 0);
assert(out == std::trunc(out));
return out;
}
} // namespace internal.
template<typename T>
T
distance(T from, T to)
{
// The distance in ulps between two numbers x and y is equal to
// the number of times that `std::nextafter()` needs to be called
// in order to go from x to y.
// The result is to become inaccurate if the difference between the
// two numbers is of at least `2^std::numeric_limits<T>::digits`.
if (std::isnan(from) || std::isnan(to)) {
return std::numeric_limits<T>::quiet_NaN();
}
if (from == to) {
return 0;
}
if (from > to) {
std::swap(from, to);
}
if (std::isinf(from)) {
// `from` equals to -infinity since `from < to`.
return 1 + distance(-std::numeric_limits<T>::max(), to);
}
if (std::isinf(to)) {
// `to` equals to +infinity since `from < to`.
return 1 + distance(from, std::numeric_limits<T>::max());
}
if (std::signbit(from) && std::signbit(to)) {
// `from` and `to` are negative.
T positiveFrom = -from;
from = -to;
to = positiveFrom;
}
static constexpr T smallestSubnormal = std::numeric_limits<T>::denorm_min();
// If either number to compare is null, the exponent returned by
// `std::logb()` will be infinity. Hence we need to compute the exponent
// of the nearest value from zero of the same sign then add one to
// the result to make up for the distance from zero to the nearest value.
if (from == 0) {
// `to` has to be positive since `from < to` and `from == 0`.
if (to == smallestSubnormal) {
return 1;
}
return 1 + internal::distance(smallestSubnormal, to);
}
if (to == 0) {
// `from` has to be negative since `from < to` and `to == 0`.
if (from == -smallestSubnormal) {
return 1;
}
return 1 + internal::distance(smallestSubnormal, -from);
}
// When the signs of the values to compare differ, we need to compare
// the distance between each number with the nearest value from zero
// of the same sign, then sum the results with an extra increment of
// two to make up for the distances from zero to the nearest values.
if (std::signbit(from) != std::signbit(to)) {
// `from` has to be negative and `to` positive since `from < to`
// and their sign differ.
return 2
+ (from == -smallestSubnormal ? 0 :
internal::distance(smallestSubnormal, -from))
+ (to == smallestSubnormal ? 0 :
internal::distance(smallestSubnormal, to));
}
return internal::distance(from, to);
}
#endif // DISTANCE_H
Raw
floatingpointtraits.h
#ifndef FLOATINGPOINTTRAITS_H
#define FLOATINGPOINTTRAITS_H
#include <limits>
namespace traits
{
namespace internal
{
template<typename T_Scalar>
constexpr T_Scalar
powerOfTwo(T_Scalar exponent)
{
return exponent == 0 ? 1 :
exponent > 0 ? (2 * powerOfTwo(exponent - 1)) :
(1 / (2 * powerOfTwo(-exponent - 1)));
}
} // namespace internal.
template<typename T>
struct FloatingPoint
{
// The floating-point numbers are treated here as being defined in the form
// `value = significand * 2^exponent` with the significand value ranging
// between 1 and `std::numeric_limits<T>::radix` (usually 2).
//
// As a result, the exponents for normalized floating-point numbers are
// defined between `std::numeric_limits<T>::min_exponent - 1` and
// `std::numeric_limits<T>::max_exponent - 1`.
static constexpr int smallestSubnormalNumberExponent =
std::numeric_limits<T>::min_exponent - std::numeric_limits<T>::digits;
static constexpr int smallestNormalNumberExponent =
std::numeric_limits<T>::min_exponent - 1;
static constexpr int largestNumberExponent =
std::numeric_limits<T>::max_exponent - 1;
// Within the subnormal numbers range, incrementing the exponent
// doubles the amount of numbers that can be represented within
// two successive power of two but the epsilon value
// between two adjacent numbers remains constant.
static constexpr T subnormalNumbersEpsilon =
std::numeric_limits<T>::denorm_min();
// Within the normal numbers range, incrementing the exponent
// halves the precision of the epsilon value between two adjacent numbers
// but the amount of numbers that can be represented
// within two successive power of two remains constant.
static constexpr T amountOfNormalNumbersPerExponent =
internal::powerOfTwo(T(std::numeric_limits<T>::digits - 1));
};
} // namespace traits.
#endif // FLOATINGPOINTTRAITS_H
Raw
test.cpp
#include <algorithm>
#include <cmath>
#include <cstddef>
#include <iostream>
#include <limits>
#include <vector>
#include "distance.h"
// TestDataHelper struct definition.
//
// O-(''Q)
// -----------------------------------------------------------------------------
template<typename T>
struct TestDataHelper;
template<>
struct TestDataHelper<float>
{
static std::vector<float>
randomValues()
{
return {
// Small non-subnormal positive values.
+1.39707453e-38f, +2.21643714e-37f, +1.53309747e-36f, +5.44752738e-36f,
+7.06794505e-35f, +1.84537148e-33f, +2.08331712e-31f, +1.02366176e-30f,
+2.29569601e-29f, +5.65505465e-28f, +1.77006721e-26f, +1.89976750e-25f,
+1.22189279e-24f, +5.86009623e-24f, +1.30081472e-21f, +3.40062714e-21f,
+3.83077371e-19f, +1.36280487e-18f, +1.04683161e-16f, +1.55927257e-15f,
+3.97193697e-15f, +4.05875637e-13f, +4.18267244e-12f, +1.74538890e-11f,
+1.83292637e-09f, +9.54729718e-09f, +4.42461356e-07f, +7.11808434e-06f,
+5.86588685e-05f, +4.04127757e-04f, +5.06849261e-04f, +4.82005589e-02f,
// Small non-subnormal negative values.
-2.14937100e-38f, -2.71179122e-37f, -9.57290701e-37f, -6.31573370e-35f,
-1.47271295e-33f, -3.28248569e-33f, -2.17465867e-31f, -4.65401272e-30f,
-6.75446923e-29f, -2.70445390e-28f, -6.56840642e-27f, -3.70014616e-25f,
-4.95167103e-25f, -1.26538959e-24f, -1.05587212e-21f, -4.36966519e-21f,
-1.84314271e-19f, -4.65708227e-18f, -3.73290762e-17f, -1.00099356e-16f,
-1.74559379e-14f, -1.22694416e-13f, -4.03344762e-12f, -8.98345634e-11f,
-1.48624490e-09f, -1.78331998e-08f, -7.09012298e-08f, -9.00824162e-08f,
-1.46330603e-05f, -2.85773694e-05f, -1.58181638e-02f, -7.58057311e-02f,
// Non-subormal positive values.
+4.88976270e-01f, +1.19546652e+01f, +1.89442337e+02f, +8.07529480e+02f,
+5.27478125e+04f, +1.03697994e+06f, +1.01791550e+07f, +8.48552560e+07f,
+1.94855168e+09f, +4.80912835e+10f, +1.86110296e+11f, +9.40841003e+12f,
+2.46404077e+13f, +6.87113700e+14f, +5.04988940e+16f, +4.11258993e+17f,
+5.15260451e+18f, +2.49693391e+20f, +4.47107576e+21f, +6.11249438e+22f,
+7.44703633e+23f, +3.87430826e+24f, +1.32378945e+26f, +2.85050748e+27f,
+2.39292746e+28f, +7.68203072e+29f, +1.08868013e+31f, +3.10118796e+32f,
+4.07830897e+33f, +2.71865918e+34f, +8.56937906e+35f, +3.12659492e+36f,
// Non-subnormal negative values.
-3.39485437e-01f, -1.56564827e+01f, -2.16126160e+02f, -3.25303986e+02f,
-5.32189492e+04f, -4.52800750e+05f, -1.61457820e+07f, -5.02500000e+07f,
-1.03071405e+09f, -2.36702618e+10f, -4.21891473e+11f, -2.82603238e+12f,
-3.93928486e+13f, -1.50511237e+15f, -1.39496994e+16f, -6.45580589e+17f,
-2.24618928e+18f, -1.69562444e+20f, -1.47329562e+21f, -2.74004879e+22f,
-6.26370868e+23f, -8.40602217e+24f, -2.16801609e+25f, -2.51693132e+27f,
-6.08741020e+27f, -3.05793630e+28f, -3.38248861e+30f, -8.34515207e+31f,
-2.16951298e+33f, -4.12623712e+34f, -8.95593406e+35f, -7.85781113e+36f
};
}
};
template<>
struct TestDataHelper<double>
{
static std::vector<double>
randomValues()
{
return {
// Small non-subnormal positive values.
+3.0604972324194311e-307, +6.7491369930635272e-298,
+1.2846638073575816e-288, +5.0247553419568776e-279,
+4.0144063431822133e-269, +4.7459035055565366e-259,
+2.0401907596490763e-249, +3.8025168520181211e-240,
+3.6456045197196816e-230, +9.9241634710788449e-221,
+6.8893689662980400e-211, +2.9837649077685835e-201,
+5.2167466997989244e-192, +4.7992688077867061e-182,
+5.3891808439698168e-173, +7.7960799668029759e-163,
+1.0052256981883652e-153, +1.2498559560086430e-143,
+3.5067421791398320e-134, +1.2282650562660921e-124,
+1.5362642475071183e-114, +3.2875973732195868e-105,
+2.1746266038685319e-095, +8.9533266943042601e-086,
+3.5714233379701883e-077, +1.3852159133877683e-066,
+2.0516826485713844e-057, +2.7477890999323729e-047,
+1.0679200883247513e-037, +2.8345002387170384e-028,
+1.9521560393356937e-018, +3.1777674228585397e-009,
// Small non-subnormal negative values.
-1.7383605493941595e-307, -9.6956074769945193e-298,
-4.9770852119977892e-289, -1.4669491027954632e-278,
-8.6915775146338386e-269, -1.1106655247588516e-259,
-7.3168183393325645e-250, -5.8254071371361915e-240,
-1.1338286119438010e-230, -3.0240742312836975e-221,
-1.3741160322166637e-211, -1.4525561283802520e-201,
-1.0384409636816511e-193, -5.7109105329562271e-182,
-1.5953227983048856e-172, -7.8909694679908292e-163,
-3.4873707865917991e-153, -1.5588070924250605e-143,
-1.5736907380797844e-134, -1.9038270734370708e-124,
-1.1563293298544775e-114, -4.4382880589977151e-105,
-2.3367479279479162e-095, -1.0532705449454580e-085,
-1.4665270599979929e-076, -8.3948246520387701e-067,
-6.8060547453816635e-057, -1.0392267729858452e-048,
-1.1687404942662141e-038, -4.0779083987473568e-028,
-1.6219008512955870e-018, -3.0783231043952328e-009,
// Non-subormal positive values.
+6.7312446103573420e-001, +3.8858757947132740e+009,
+1.1455964479950627e+019, +2.6592437869848886e+028,
+2.4985860075237618e+038, +1.6109623409158992e+047,
+4.2369807664889535e+057, +1.5637549761414510e+067,
+7.4983403013310807e+076, +2.5109324402578082e+085,
+5.9868103724342808e+095, +8.9987155805842543e+105,
+2.8203223088335084e+115, +3.1377134127018829e+124,
+3.5156026976269647e+134, +1.4779277452464080e+144,
+9.4415568981768957e+152, +1.6073698558667417e+163,
+1.4004080019487282e+173, +9.3136105930023645e+182,
+1.0118520940275815e+192, +5.5297416988082782e+201,
+4.7584413460502903e+211, +2.9796800105717338e+221,
+1.3362577067396765e+231, +3.2818907351504651e+240,
+2.1071732315449056e+250, +5.5629934491690242e+258,
+2.0822344175713011e+269, +7.6580761754503568e+278,
+4.7898170200223161e+288, +2.3128750888569858e+297,
// Non-subnormal negative values.
-5.8704355499724636e-001, -1.3500635827495389e+007,
-4.6455188847831931e+018, -4.5309514619112853e+028,
-7.8194319617532854e+037, -7.7119759867394170e+047,
-1.7045641726416442e+057, -1.5371363224427223e+067,
-6.0620078249192394e+076, -3.0144771120153138e+086,
-1.1470814505749805e+096, -8.3868019711026089e+105,
-3.7375577621073641e+114, -1.9425765657233467e+124,
-1.4793740239794770e+134, -2.7443574941121530e+144,
-5.7205791700491208e+153, -3.6879548517134091e+163,
-1.8016128006371260e+173, -3.8470118656198290e+182,
-2.3126990253624127e+191, -2.6678985854623137e+201,
-4.6002717010861424e+211, -7.9764470305111089e+220,
-1.4290597280721725e+230, -3.7292810617649615e+240,
-2.0089528530568632e+250, -1.6183160488902539e+259,
-4.1547168751415793e+269, -2.2123849505286914e+278,
-7.9538043916486098e+288, -3.1152582179840712e+298
};
}
};
template<>
struct TestDataHelper<long double>
{
static std::vector<long double>
randomValues()
{
return {
// Small non-subnormal positive values.
+3.58282670787618450878420909559674e-4932L,
+1.62883629625361471571712463377950e-4777L,
+1.69441374344082021886576182279463e-4623L,
+2.69916949090091731381431631791078e-4469L,
+2.36536866224627153547619921004631e-4315L,
+2.89697256717443857344884862677359e-4161L,
+3.77460085768415251075870818123390e-4007L,
+1.77989848895537471705686032347030e-3853L,
+1.16413968270339948915642498161788e-3698L,
+1.46942053260383916091252013764818e-3544L,
+8.85970262002374548436429553034007e-3391L,
+1.04826442877821464413699353151200e-3236L,
+3.74937276875201167011053225397026e-3082L,
+2.79109522698032122597649947585959e-2928L,
+2.80091695784469832083775134207639e-2774L,
+8.93348883926216941492930024802830e-2620L,
+8.12350667123689300115156782030686e-2466L,
+5.01446717346937767719739318200038e-2312L,
+3.57812396924381463969152079583015e-2158L,
+3.01768684799339966738642785957235e-2004L,
+9.07361796275422773687026208216255e-1850L,
+3.00306635466183849490092161878697e-1695L,
+1.88146409338356252313196228793816e-1542L,
+1.05972294475657602584954056187726e-1386L,
+1.12843384874826996859495088588033e-1232L,
+7.91818278375874936729442396782907e-1079L,
+1.14339354812653795328828966787417e-0924L,
+1.43404875444924498025974694010899e-0770L,
+2.22622926552606276354044656392614e-0616L,
+3.40691175671305633931996026204289e-0462L,
+1.02287657890247897410837092017875e-0308L,
+9.39482208358113538806750841134851e-0154L,
// Small non-subnormal negative values.
-7.95854562181148464398082777459519e-4932L,
-6.22473208329651542592376552525345e-4778L,
-1.59828680785798362362525663639723e-4623L,
-2.94609364860936828666180088020532e-4469L,
-5.93458172190120968067616843705561e-4316L,
-5.54810901370806640508099414028459e-4163L,
-2.27601643120846017040650409867248e-4007L,
-6.11898699636314122912007441835105e-3853L,
-9.44419908064015116446221936787634e-3699L,
-1.44671903820800985503153702874234e-3544L,
-2.45145547766773694263896644179753e-3390L,
-8.39309591132445607397453815237565e-3237L,
-2.45502717936137010666000092017051e-3082L,
-5.72424601551190375477980222469233e-2928L,
-1.19943121840992513223520046278285e-2774L,
-7.07148308639403273872628905995017e-2620L,
-1.25971781536015552125985668481582e-2465L,
-5.59814968167226860779377330471251e-2312L,
-2.36615172488543519406672259815460e-2158L,
-2.84753324367803088833302796266704e-2003L,
-8.38466745243769523440223947712476e-1850L,
-4.33452691530985023364379201860311e-1695L,
-2.30841390709073700804913811922985e-1541L,
-8.69324741960651061253783525167992e-1387L,
-3.46248724517850957034519136088992e-1233L,
-1.85570222524090605615601857065524e-1079L,
-6.67983343779271540922515704553626e-0925L,
-3.50065635294276473508139906308249e-0770L,
-8.36959296195393587017024879318322e-0618L,
-5.68622678677639244952791108185755e-0462L,
-5.92277589890811315016823838706879e-0308L,
-4.88236194831246362038973367483548e-0154L,
// Non-subormal positive values.
+7.75451667378974193582143181746602e-0001L,
+8.18747764860630142734663937327801e+0153L,
+6.75650897956847698003406335895887e+0307L,
+2.06645993387854745597867156866492e+0462L,
+8.51649039451405783712258449808567e+0615L,
+1.83207742524216690445026159846819e+0770L,
+2.74634419669010670153328353012089e+0924L,
+6.86956651599906249309415771283356e+1078L,
+4.57905508348972992219608199958868e+1232L,
+9.63439909963123682720894258650396e+1386L,
+9.79311112661649927556238875593887e+1540L,
+2.46889066901589191310543163264119e+1695L,
+7.30546769457199332385742494967952e+1848L,
+4.14389882294952567831279721597136e+2003L,
+4.80710819974574571209032428122381e+2156L,
+2.67985278591945207304181353277517e+2311L,
+4.41517511906402766020183775613779e+2465L,
+9.52924674401985469851847208443548e+2619L,
+1.44336835408402149899154876933275e+2774L,
+5.61942338569628350672613997318828e+2927L,
+3.01557925967236314546330369390672e+3080L,
+3.36585320075309151721639085708904e+3236L,
+5.63498890448744403067210692383959e+3390L,
+4.68767643022242570206021013958514e+3544L,
+9.12469486888258807869869841577778e+3698L,
+1.30579257836396511174164009821036e+3853L,
+6.88238225645436543058264485284186e+4006L,
+3.69187910982698170390472104264172e+4160L,
+1.18744288766699830930908737176785e+4315L,
+3.85319067318926462378055714222182e+4469L,
+6.28556133801454349813291497176754e+4623L,
+4.83035690090808708121180695952860e+4777L,
// Non-subnormal negative values.
-1.36292775046919707553332717919758e-0001L,
-4.80954610372164205569586907066979e+0152L,
-9.79128101795403237698192511941189e+0307L,
-2.39529386231975978599491649552795e+0462L,
-1.26495223194752910363594559918618e+0616L,
-2.53493611684741600958793090441123e+0770L,
-1.42207802702287601232538815503980e+0924L,
-4.12699866174710633532633847847134e+1078L,
-5.62281420395745621397323955134295e+1232L,
-6.27938213667111819363873693067516e+1386L,
-9.62783763918038306662273487588752e+1540L,
-2.23427273726655093305887833282419e+1695L,
-4.73502996009396339783632659654225e+1847L,
-4.38439981203564321165180594555575e+2003L,
-3.81086623572480393767852254645586e+2157L,
-5.95365622488650089116015163124686e+2311L,
-3.04928657162836443324297752604405e+2464L,
-3.17144802208917896133100818060613e+2619L,
-2.17754574253693438291920938967984e+2773L,
-9.38966030651618989702737182380971e+2927L,
-9.16691562712588123972775853227657e+3081L,
-3.50013751606097730963867881121116e+3236L,
-3.93577009995343245197904699407424e+3390L,
-1.78787756130996009297617245911249e+3544L,
-6.80121998422414203638526785130998e+3698L,
-1.03075862267472006807661601397973e+3853L,
-1.20952197453021433913654536242583e+4007L,
-2.06622873381076943007676980837342e+4161L,
-1.14087912585522638587751490485598e+4314L,
-1.83315320024850014230760531718997e+4469L,
-5.11276493498405241258159676881630e+4623L,
-3.82961368007045063583055724804425e+4777L
};
}
};
// TestData struct definition.
//
// O-(''Q)
// -----------------------------------------------------------------------------
template<typename T>
struct TestData
{
static std::vector<T>
handPickedValues()
{
return {
// Zeroes.
+0.0,
-0.0,
// Powers of two.
+0.25, +0.5, +1, +2, +4, +16, +256, +65536,
-0.25, -0.5, -1, -2, -4, -16, -256, -65536,
// Powers of ten.
+10, +100, +1000, +10000, +100000, +1000000,
-10, -100, -1000, -10000, -100000, -1000000,
// Normal numbers.
+0.1, +0.7, +3, +7, +9, +22, +328, +83285,
-0.1, -0.7, -3, -7, -9, -22, -328, -83285,
// Subnormal numbers.
+std::numeric_limits<T>::denorm_min(),
-std::numeric_limits<T>::denorm_min(),
+std::numeric_limits<T>::min() / 2,
-std::numeric_limits<T>::min() / 2,
// Subnormal to normal boundary numbers.
+std::numeric_limits<T>::min(),
-std::numeric_limits<T>::min()
};
}
static std::vector<T>
randomValues()
{
return TestDataHelper<T>::randomValues();
}
static std::vector<T>
values()
{
std::vector<T> handPickedValues = TestData::handPickedValues();
std::vector<T> randomValues = TestData::randomValues();
std::vector<T> out;
out.reserve(handPickedValues.size() + randomValues.size());
out.insert(out.begin(), handPickedValues.begin(), handPickedValues.end());
out.insert(out.end(), randomValues.begin(), randomValues.end());
return out;
}
};
// Helper functions.
//
// O-(''Q)
// -----------------------------------------------------------------------------
template<typename T>
T
_next(T value)
{
return std::nextafter(value, std::numeric_limits<T>::infinity());
}
template<typename T>
T
_previous(T value)
{
return std::nextafter(value, -std::numeric_limits<T>::infinity());
}
template<typename T>
T
_increment(T value, unsigned int shift)
{
while (shift-- > 0) {
value = std::nextafter(value, std::numeric_limits<T>::infinity());
}
return value;
}
template<typename T>
T
_decrement(T value, unsigned int shift)
{
while (shift-- > 0) {
value = std::nextafter(value, -std::numeric_limits<T>::infinity());
}
return value;
}
// Test suite.
//
// O-(''Q)
// -----------------------------------------------------------------------------
template<typename T>
void distance()
{
std::vector<T> shifts = {
0, 1, 2, 3, 4, 5, 7, 9, 11, 16, 23, 54, 89, 123, 249, 1137
};
for (auto value : TestData<T>::values()) {
assert(distance(value, value) == 0);
assert(distance(value, _next(value)) == 1);
assert(distance(value, _previous(value)) == 1);
assert(distance(_next(value), value) == 1);
assert(distance(_previous(value), value) == 1);
assert(distance(_next(value), _next(value)) == 0);
assert(distance(_previous(value), _previous(value)) == 0);
assert(distance(_next(value), _previous(value)) == 2);
assert(distance(_previous(value), _next(value)) == 2);
for (auto shift : shifts) {
assert(distance(value, _increment(value, shift)) == shift);
assert(distance(value, _decrement(value, shift)) == shift);
assert(distance(_increment(value, shift), value) == shift);
assert(distance(_decrement(value, shift), value) == shift);
assert(distance(_increment(value, shift), _increment(value, shift)) == 0);
assert(distance(_decrement(value, shift), _decrement(value, shift)) == 0);
assert(distance(_increment(value, shift), _decrement(value, shift)) == shift * 2);
assert(distance(_decrement(value, shift), _increment(value, shift)) == shift * 2);
}
}
int largestBiasedExponent =
2 * std::numeric_limits<T>::max_exponent - 1;
T amountOfSubnormalNumbers =
std::exp2(T(std::numeric_limits<T>::digits - 1));
T amountOfNormalNumbers =
(largestBiasedExponent - 1)
* std::exp2(T(std::numeric_limits<T>::digits));
T amountOfNumbers =
largestBiasedExponent
* std::exp2(T(std::numeric_limits<T>::digits));
// Some distances below are so large that the comparison tests
// are highly inaccurate.
assert(std::isnan(distance(-std::numeric_limits<T>::quiet_NaN(), -std::numeric_limits<T>::quiet_NaN())));
assert(distance(-std::numeric_limits<T>::infinity(), -std::numeric_limits<T>::infinity()) == 0);
assert(distance(-std::numeric_limits<T>::infinity(), T(0)) == amountOfNumbers / 2);
assert(distance(-std::numeric_limits<T>::infinity(), +std::numeric_limits<T>::infinity()) == amountOfNumbers);
assert(distance(-std::numeric_limits<T>::max(), -std::numeric_limits<T>::infinity()) == 1);
assert(distance(-std::numeric_limits<T>::max(), -std::numeric_limits<T>::min()) == amountOfNormalNumbers / 2);
assert(distance(-std::numeric_limits<T>::max(), T(0)) == (amountOfNumbers - 1) / 2);
assert(distance(-std::numeric_limits<T>::max(), +std::numeric_limits<T>::max()) == amountOfNumbers - 2);
assert(distance(-std::numeric_limits<T>::min(), T(0)) == amountOfSubnormalNumbers);
assert(distance(-std::numeric_limits<T>::min(), +std::numeric_limits<T>::min()) == amountOfSubnormalNumbers * 2);
assert(std::isnan(distance(T(0), +std::numeric_limits<T>::quiet_NaN())));
assert(distance(+std::numeric_limits<T>::min(), -std::numeric_limits<T>::min()) == amountOfSubnormalNumbers * 2);
assert(distance(+std::numeric_limits<T>::min(), T(0)) == amountOfSubnormalNumbers);
assert(distance(+std::numeric_limits<T>::max(), -std::numeric_limits<T>::max()) == amountOfNumbers - 2);
assert(distance(+std::numeric_limits<T>::max(), T(0)) == (amountOfNumbers - 1) / 2);
assert(distance(+std::numeric_limits<T>::max(), +std::numeric_limits<T>::min()) == amountOfNormalNumbers / 2);
assert(distance(+std::numeric_limits<T>::max(), +std::numeric_limits<T>::infinity()) == 1);
assert(distance(+std::numeric_limits<T>::infinity(), -std::numeric_limits<T>::infinity()) == amountOfNumbers);
assert(distance(+std::numeric_limits<T>::infinity(), T(0)) == amountOfNumbers / 2);
assert(distance(+std::numeric_limits<T>::infinity(), +std::numeric_limits<T>::infinity()) == 0);
assert(std::isnan(distance(+std::numeric_limits<T>::quiet_NaN(), +std::numeric_limits<T>::quiet_NaN())));
}
// Main procedure.
//
// O-(''Q)
// -----------------------------------------------------------------------------
int main(int argc, char **argv)
{
distance<float>();
distance<double>();
distance<long double>();
std::cout << "All good!" << std::endl;
return 0;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment